Improve performance of add/sub when exponent of two bigdecimals have huge difference (#478)

* Improve performance of add/sub when exponent of two bigdecimals have huge difference

BigDecimal('1e+10000000000').add(1, 10) allocates GB of memory.
BigDecimal('1e+10000000000').add(BigDecimal('0.1e+9999999990'), 10) will give the same rounding result with less memory allocation and computation time.

* Remove workaround of add/sub performance bug

Author: tompng

ext/bigdecimal/bigdecimal.c

@@ -986,7 +986,7 @@ BigDecimal_mode(int argc, VALUE *argv, VALUE self)
 static size_t
 GetAddSubPrec(Real *a, Real *b)
 {
-    if (!VpIsDef(a) || !VpIsDef(b)) return (size_t)-1L;
+    if (VpIsZero(a) || VpIsZero(b)) return Max(a->Prec, b->Prec);
     ssize_t min_a = a->exponent - a->Prec;
     ssize_t min_b = b->exponent - b->Prec;
     return Max(a->exponent, b->exponent) - Min(min_a, min_b);
@@ -1290,13 +1290,32 @@ BigDecimal_addsub_with_coerce(VALUE self, VALUE r, size_t prec, int operation)
     if (VpIsNaN(a.real)) return CheckGetValue(a);
     if (VpIsNaN(b.real)) return CheckGetValue(b);
 
-    mx = GetAddSubPrec(a.real, b.real);
-    if (mx == (size_t)-1L) {
-        /* a or b is inf */
+    if (VpIsInf(a.real) || VpIsInf(b.real)) {
         c = NewZeroWrap(1, BASE_FIG);
         VpAddSub(c.real, a.real, b.real, operation);
     }
     else {
+
+        // Optimization when exponent difference is large
+        // (1.234e+1000).add(5.678e-1000, 10) == (1.234e+1000).add(0.1e+990, 10) in every rounding mode
+        if (prec && !VpIsZero(a.real) && !VpIsZero(b.real)) {
+            size_t precRoom = roomof(prec, BASE_FIG);
+            if (a.real->exponent - (ssize_t)Max(a.real->Prec, precRoom) - 1 > b.real->exponent) {
+                BDVALUE b2 = NewZeroWrap(1, BASE_FIG);
+                VpSetOne(b2.real)
+                VpSetSign(b2.real, b.real->sign);
+                b2.real->exponent = a.real->exponent - (ssize_t)Max(a.real->Prec, precRoom) - 1;
+                b = b2;
+            } else if (b.real->exponent - (ssize_t)Max(b.real->Prec, precRoom) - 1 > a.real->exponent) {
+                BDVALUE a2 = NewZeroWrap(1, BASE_FIG);
+                VpSetOne(a2.real)
+                VpSetSign(a2.real, a.real->sign);
+                a2.real->exponent = b.real->exponent - (ssize_t)Max(b.real->Prec, precRoom) - 1;
+                a = a2;
+            }
+        }
+
+        mx = GetAddSubPrec(a.real, b.real);
         c = NewZeroWrap(1, (mx + 1) * BASE_FIG);
         size_t pl = VpGetPrecLimit();
         if (prec) VpSetPrecLimit(prec);

lib/bigdecimal/math.rb

@@ -574,14 +574,7 @@ def expm1(x, prec)
 
     return BigDecimal(-1) if exp_prec <= 0
 
-    exp = exp(x, exp_prec)
-
-    if exp.exponent > prec + BigDecimal.double_fig
-      # Workaroudn for https://github.com/ruby/bigdecimal/issues/464
-      exp
-    else
-      exp.sub(1, prec)
-    end
+    exp(x, exp_prec).sub(1, prec)
   end
 
   #   erf(decimal, numeric) -> BigDecimal
@@ -608,12 +601,7 @@ def erf(x, prec)
       log10_erfc = -xf ** 2 / Math.log(10) - Math.log10(xf * Math::PI ** 0.5)
       erfc_prec = [prec + log10_erfc.ceil, 1].max
       erfc = _erfc_asymptotic(x, erfc_prec)
-      if erfc
-        # Workaround for https://github.com/ruby/bigdecimal/issues/464
-        return BigDecimal(1) if erfc.exponent < -prec - BigDecimal.double_fig
-
-        return BigDecimal(1).sub(erfc, prec)
-      end
+      return BigDecimal(1).sub(erfc, prec) if erfc
     end
 
     prec2 = prec + BigDecimal.double_fig

test/bigdecimal/test_bigdecimal.rb

@@ -523,6 +523,52 @@ def test_add_sub_underflow
     assert_negative_zero((-x) - (-y))
   end
 
+  def test_add_sub_huge_with_zero
+    huge = BigDecimal("0.1e#{EXPONENT_MAX}")
+    assert_equal(huge, huge + 0)
+    assert_equal(huge, huge - 0)
+    assert_equal(huge, huge.add(BigDecimal(0), 10))
+    assert_equal(huge, huge.sub(BigDecimal(0), 10))
+  end
+
+  def test_add_sub_huge_exponent_difference
+    huge = BigDecimal('1e+1000000000000')
+    zero = BigDecimal(0)
+    small = BigDecimal('1e-1000000000000')
+    assert_equal(huge, huge.add(small, 10))
+    assert_equal(huge, huge.sub(small, 10))
+    assert_equal(huge, small.add(huge, 10))
+    assert_equal(-huge, small.sub(huge, 10))
+    assert_equal(huge, zero.add(huge, 10))
+    assert_equal(huge, huge.add(zero, 10))
+    assert_equal(small, zero.add(small, 10))
+    assert_equal(small, small.add(zero, 10))
+
+    BigDecimal.mode(BigDecimal::ROUND_MODE, BigDecimal::ROUND_DOWN)
+    assert_equal(huge, huge.add(small, 10))
+    assert_equal(huge * BigDecimal('0.9999999999'), huge.sub(small, 10))
+    assert_equal(huge, huge.add(zero, 10))
+    assert_equal(huge, huge.sub(zero, 10))
+
+    BigDecimal.mode(BigDecimal::ROUND_MODE, BigDecimal::ROUND_UP)
+    assert_equal(huge * BigDecimal('1.000000001'), huge.add(small, 10))
+    assert_equal(huge, huge.sub(small, 10))
+    assert_equal(huge, huge.add(zero, 10))
+    assert_equal(huge, huge.sub(zero, 10))
+
+    [BigDecimal::ROUND_UP, BigDecimal::ROUND_DOWN, BigDecimal::ROUND_HALF_UP].each do |mode|
+      BigDecimal.mode(BigDecimal::ROUND_MODE, mode)
+      xs = [BigDecimal(1), BigDecimal(2)]
+      ys = [BigDecimal('0.5e-89'), BigDecimal('0.5e-90'), BigDecimal('0.5e-91')]
+      xs.product(ys).each do |x, y|
+        assert_equal((x + y).mult(1, 91), x.add(y, 91))
+        assert_equal((y + x).mult(1, 91), y.add(x, 91))
+        assert_equal((x - y).mult(1, 91), x.sub(y, 91))
+        assert_equal((y - x).mult(1, 91), y.sub(x, 91))
+      end
+    end
+  end
+
   def test_mult_div_overflow_underflow_sign
     BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, false)
     BigDecimal.mode(BigDecimal::EXCEPTION_UNDERFLOW, false)