Switch Lagrange interpolation and Spouge approximation

Author: tompng

lib/bigdecimal/math.rb

@@ -737,6 +737,8 @@ def gamma(x, prec)
       return pi.div(gamma(1 - x, prec2).mult(sin, prec2), prec)
     elsif x.frac.zero? && x < 1000 * prec
       return _gamma_positive_integer(x, prec2).mult(1, prec)
+    elsif x < (prec / prec.bit_length)**2
+      return _gamma_lagrange(x, prec2).mult(1, prec)
     end
 
     a, sum = _gamma_spouge_sum_part(x, prec2)
@@ -764,8 +766,9 @@ def gamma(x, prec)
     [prd, coef]
   end
 
-  def gamma_lagrange(x, prec)
+  private_class_method def _gamma_lagrange(x, prec)
     # Calculate approximage gamma by Lagrange interpolation of b**x/x! at x=b-l, b-l+1, ..., b+l
+    # Estimated time compexity: O(prec**2), precondition: 0 < x < const * (prec / prec.bit_length)**2
 
     x = BigDecimal(x) - 1
     l = prec
@@ -904,8 +907,13 @@ def lgamma(x, prec)
       end
 
       prec2 += [-diff1_exponent, -diff2_exponent, 0].max
-      a, sum = _gamma_spouge_sum_part(x, prec2)
-      log_gamma = BigMath.log(sum, prec2).add((x - 0.5).mult(BigMath.log(x.add(a - 1, prec2), prec2), prec2) + 1 - x, prec)
+
+      if x < (prec / prec.bit_length)**2
+        log_gamma = BigMath.log(_gamma_lagrange(x, prec2), prec)
+      else
+        a, sum = _gamma_spouge_sum_part(x, prec2)
+        log_gamma = BigMath.log(sum, prec2).add((x - 0.5).mult(BigMath.log(x.add(a - 1, prec2), prec2), prec2) + 1 - x, prec)
+      end
       [log_gamma, 1]
     end
   end

test/bigdecimal/test_bigmath.rb

@@ -567,11 +567,11 @@ def test_gamma
       BigDecimal('0.28242294079603478742934215780245355184774949260912e456569'),
       BigMath.gamma(100000, 50)
     )
-    precisions = [50, 100, 150]
-    assert_converge_in_precision(precisions) {|n| gamma(BigDecimal("0.3"), n) }
-    assert_converge_in_precision(precisions) {|n| gamma(BigDecimal("-1.9" + "9" * 30), n) }
-    assert_converge_in_precision(precisions) {|n| gamma(BigDecimal("1234.56789"), n) }
-    assert_converge_in_precision(precisions) {|n| gamma(BigDecimal("-987.654321"), n) }
+    precisions = [100, 200, 300, 400]
+    assert_converge_in_precision() {|n| gamma(BigDecimal("0.3"), n) }
+    assert_converge_in_precision() {|n| gamma(BigDecimal("-1.9" + "9" * 30), n) }
+    assert_converge_in_precision() {|n| gamma(BigDecimal("1234.56789"), n) }
+    assert_converge_in_precision() {|n| gamma(BigDecimal("-987.654321"), n) }
   end
 
   def test_lgamma
@@ -587,7 +587,7 @@ def test_lgamma
       assert_equal(sign, bigsign)
     end
     assert_equal([BigMath.log(PI(120).sqrt(120), 100), 1], lgamma(BigDecimal("0.5"), 100))
-    precisions = [50, 100, 150]
+    precisions = [50, 100, 150, 200]
     assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("0." + "9" * 80), n).first }
     assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("1." + "0" * 80 + "1"), n).first }
     assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("1." + "9" * 80), n).first }