BigMath.exp overflow/underflow check (#523)

* BigMath.exp overflow/underflow check

`BigMath.exp(BigDecimal('1e100000000'), 20)` was slow.
Skip calculation when x.exponent >= 21 which always result in overflow or underflow.

* Add underflow check to BigMath.erfc

Author: tompng

lib/bigdecimal.rb

@@ -57,6 +57,13 @@ def self.infinity_computation_result # :nodoc:
       BigDecimal::INFINITY
     end
 
+    def self.underflow_computation_result # :nodoc:
+      if BigDecimal.mode(BigDecimal::EXCEPTION_ALL).anybits?(BigDecimal::EXCEPTION_UNDERFLOW)
+        raise FloatDomainError, 'Exponent underflow'
+      end
+      BigDecimal(0)
+    end
+
     def self.nan_computation_result # :nodoc:
       if BigDecimal.mode(BigDecimal::EXCEPTION_ALL).anybits?(BigDecimal::EXCEPTION_NaN)
         raise FloatDomainError, "Computation results to 'NaN'"
@@ -350,7 +357,17 @@ def exp(x, prec)
     prec = BigDecimal::Internal.coerce_validate_prec(prec, :exp)
     x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :exp)
     return BigDecimal::Internal.nan_computation_result if x.nan?
-    return x.positive? ? BigDecimal::Internal.infinity_computation_result : BigDecimal(0) if x.infinite?
+    if x.infinite? || x.exponent >= 21 # exp(10**20) and exp(-10**20) overflows/underflows 64-bit exponent
+      if x.positive?
+        return BigDecimal::Internal.infinity_computation_result
+      elsif x.infinite?
+        # exp(-Infinity) is +0 by definition, this is not an underflow.
+        return BigDecimal(0)
+      else
+        return BigDecimal::Internal.underflow_computation_result
+      end
+    end
+
     return BigDecimal(1) if x.zero?
 
     # exp(x * 10**cnt) = exp(x)**(10**cnt)

lib/bigdecimal/math.rb

@@ -637,7 +637,7 @@ def erfc(x, prec)
     return BigDecimal::Internal.nan_computation_result if x.nan?
     return BigDecimal(1 - x.infinite?) if x.infinite?
     return BigDecimal(1).sub(erf(x, prec + BigDecimal::Internal::EXTRA_PREC), prec) if x < 0.5
-    return BigDecimal(0) if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)
+    return BigDecimal::Internal.underflow_computation_result if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow)
 
     if x >= 8
       y = _erfc_asymptotic(x, prec)

test/bigdecimal/helper.rb

@@ -71,7 +71,7 @@ def assert_infinite_calculation(positive:)
       BigDecimal.mode(BigDecimal::EXCEPTION_INFINITY, false)
       positive ? assert_positive_infinite(yield) : assert_negative_infinite(yield)
       BigDecimal.mode(BigDecimal::EXCEPTION_INFINITY, true)
-      assert_raise_with_message(FloatDomainError, /Infinity/) { yield }
+      assert_raise_with_message(FloatDomainError, /infinity|overflow/i) { yield }
     end
   end
 
@@ -83,6 +83,19 @@ def assert_negative_infinite_calculation(&block)
     assert_infinite_calculation(positive: false, &block)
   end
 
+  def assert_underflow_calculation(accept_overflow: false)
+    BigDecimal.save_exception_mode do
+      BigDecimal.mode(BigDecimal::EXCEPTION_UNDERFLOW, false)
+      BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, false)
+      assert_equal(BigDecimal(0), yield)
+      BigDecimal.mode(BigDecimal::EXCEPTION_UNDERFLOW, true)
+      BigDecimal.mode(BigDecimal::EXCEPTION_OVERFLOW, true)
+      # Accept internal overflow (e.g. overflow calculating denominator part)
+      pattern = accept_overflow ? /underflow|overflow/i : /underflow/i
+      assert_raise_with_message(FloatDomainError, pattern) { yield }
+    end
+  end
+
   def assert_nan_calculation(&block)
     BigDecimal.save_exception_mode do
       BigDecimal.mode(BigDecimal::EXCEPTION_NaN, false)

test/bigdecimal/test_bigdecimal.rb

@@ -2296,13 +2296,24 @@ def test_exp_with_negative
   end
 
   def test_exp_with_negative_infinite
+    # exp(-infinity) is exactly zero. This is not an underflow.
     assert_equal(0, BigMath.exp(NEGATIVE_INFINITY, 20))
   end
 
   def test_exp_with_positive_infinite
     assert_positive_infinite_calculation { BigMath.exp(BigDecimal::INFINITY, 20) }
   end
 
+  def test_exp_with_overflow_underflow
+    assert_underflow_calculation { BigMath.exp(-1e+100, 20) }
+    assert_underflow_calculation { BigMath.exp(-0.9e+20, 20) }
+    assert_positive_infinite_calculation { BigMath.exp(1e+100, 20) }
+    assert_positive_infinite_calculation { BigMath.exp(0.9e+20, 20) }
+    huge = BigDecimal("0.1e#{EXPONENT_MAX / 100}")
+    assert_positive_infinite_calculation { BigMath.exp(huge, 20) }
+    assert_underflow_calculation { BigMath.exp(-huge, 20) }
+  end
+
   def test_exp_with_nan
     assert_nan_calculation { BigMath.exp(BigDecimal::NAN, 20) }
   end

test/bigdecimal/test_bigmath.rb

@@ -521,7 +521,8 @@ def test_erfc
     end
     assert_equal(0, BigMath.erfc(PINF, N))
     assert_equal(2, BigMath.erfc(MINF, N))
-    assert_equal(0, BigMath.erfc(BigDecimal('1e400'), 10))
+    assert_underflow_calculation(accept_overflow: true) { BigMath.erfc(4999999999, 10) }
+    assert_underflow_calculation { BigMath.erfc(BigDecimal('1e400'), 10) }
     assert_equal(2, BigMath.erfc(BigDecimal('-1e400'), 10))
     assert_equal(1, BigMath.erfc(BigDecimal('1e-400'), N))