Implement BigMath.log1p and BigMath.expm1

Author: tompng

lib/bigdecimal/math.rb

@@ -22,6 +22,8 @@
 #   atanh(x, prec)
 #   log2 (x, prec)
 #   log10(x, prec)
+#   log1p(x, prec)
+#   expm1(x, prec)
 #   PI  (prec)
 #   E   (prec) == exp(1.0,prec)
 #
@@ -524,6 +526,42 @@ def log10(x, prec)
     v.mult(1, prec)
   end
 
+  # call-seq:
+  #   BigMath.log1p(decimal, numeric)    -> BigDecimal
+  #
+  # Computes log(1 + decimal) to the specified number of digits of precision, +numeric+.
+  def log1p(x, prec)
+    BigDecimal::Internal.validate_prec(prec, :log1p)
+    x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log1p)
+    raise Math::DomainError, 'Out of domain argument for log1p' if x < -1
+
+    return BigMath.log(x + 1, prec)
+  end
+
+  # call-seq:
+  #   BigMath.expm1(decimal, numeric)    -> BigDecimal
+  #
+  # Computes exp(decimal) - 1 to the specified number of digits of precision, +numeric+.
+  #
+  def expm1(x, prec)
+    BigDecimal::Internal.validate_prec(prec, :expm1)
+    x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :expm1)
+    return BigDecimal(-1) if x.infinite? == -1
+
+    exp_prec = prec
+    if x < -1
+      # log10(exp(x)) = x * log10(e)
+      lg_e = 0.4342944819032518
+      exp_prec = prec + (lg_e * x).ceil + 2
+    elsif x < 1
+      exp_prec = prec - x.exponent + 2
+    else
+      exp_prec = prec
+    end
+    exp_prec > 0 ? BigMath.exp(x, exp_prec).sub(1, prec) : BigDecimal(-1)
+  end
+
+
   # call-seq:
   #   PI(numeric) -> BigDecimal
   #

test/bigdecimal/test_bigmath.rb

@@ -363,4 +363,24 @@ def test_log10
       end
     end
   end
+
+  def test_log1p
+    assert_raise(Math::DomainError) { log1p(MINF, N) }
+    assert_raise(Math::DomainError) { log1p(BigDecimal("-1.01"), N) }
+    assert_in_epsilon(Math.log(0.01), log1p(BigDecimal("-0.99"), N))
+    assert_positive_infinite_calculation { log1p(PINF, N) }
+    assert_in_exact_precision(BigMath.log(1 + BigDecimal("1e-20"), 100), log1p(BigDecimal("1e-20"), 100), 100)
+  end
+
+  def test_expm1
+    assert_equal(-1, expm1(MINF, N))
+    assert_positive_infinite_calculation { expm1(PINF, N) }
+    assert_equal(-1, expm1(BigDecimal("-400"), 100))
+    assert_equal(-1, expm1(BigDecimal("-231"), 100))
+    assert_not_equal(-1, expm1(BigDecimal("-229"), 100))
+    assert_in_exact_precision(BigMath.exp(-220, 100) - 1, expm1(BigDecimal("-220"), 100), 100)
+    assert_in_exact_precision(BigMath.exp(-3, 100) - 1, expm1(BigDecimal("-3"), 100), 100)
+    assert_in_exact_precision(BigMath.exp(BigDecimal("1.23e-10"), 120) - 1, expm1(BigDecimal("1.23e-10"), 100), 100)
+    assert_in_exact_precision(BigMath.exp(123, 120) - 1, expm1(BigDecimal("123"), 100), 100)
+  end
 end