optimize full-digit

Author: tompng

lib/bigdecimal/math.rb

@@ -743,32 +743,100 @@ def gamma(x, prec)
     (x + (a - 1)).power(x - 0.5, prec2).mult(BigMath.exp(1 - x, prec2), prec2).mult(sum, prec)
   end
 
+  # Calculates prod{x-k} and its coefficients for given ks, xn and prec with baby-step giant-step method.
+  # xn is an array of precalculated powers of x: [1, x, x**2, x**3, ...]
+  private_class_method def _x_minus_k_prod_coef(ks, xn, prec)
+    coef = [1]
+    ks.each do |k|
+      coef_next = [0] * (coef.size + 1)
+      coef.each_with_index do |c, i|
+        coef_next[i] -= k * c
+        coef_next[i + 1] += c
+      end
+      coef = coef_next
+    end
+
+    prd = coef.each_with_index.map do |c, i|
+      xn[i].mult(c, prec)
+    end.reduce do |sum, value|
+      sum.add(value, prec)
+    end
+    [prd, coef]
+  end
+
   def gamma_lagrange(x, prec)
+    # Calculate approximage gamma by Lagrange interpolation of b**x/x! at x=b-l, b-l+1, ..., b+l
+
     x = BigDecimal(x) - 1
     l = prec
     shift = x < 2 * prec ? 2 * prec - x.floor : 0
     x += shift
-    # Lagrange interpolate of b**x/x! at x=b-l, b-l+1, ..., b+l
     b = x.round
-    prods = ((b-l)..(b+l)).map {|i| x - i } + shift.times.map {|i| x - i }
-    prods = prods.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a } while prods.size != 1
-    prod = prods.first
 
     c0s = [*(1..b-l), *(1..2*l)]
     c0s = c0s.each_slice(2).map {|a, b| b ? BigDecimal(a).mult(b, prec) : BigDecimal(a) } while c0s.size != 1
     c0 = c0s.first
 
+    # Optimize this calculation for full-digit-x case and small-digit-x case
+    # sum = BigDecimal(0)
+    # prod = [*(b-l)..(b+l), *0...shift).map {|i| x - i }.reduce { _1.mult(_2, prec) }
+    # c = BigDecimal(1)
+    # ((b-l)..(b+l)).each do |i|
+    #   if i != b - l
+    #     c = c.mult(-b * (b + l - i + 1), prec).div((i - b + l) * i, prec)
+    #   end
+    #   sum = sum.add(c.div(x - i, prec), prec)
+    # end
     if x.n_significant_digits > prec / 10
-      sum = BigDecimal(0)
+      # Reduce multiplications by baby-step giant-step method
+
+      batch_size = prec.bit_length
+      # When expanding prod{x-k}, the coefficient of x**n might be huge.
+      # Increase internal calculation precision to avoid loss of precision due to cancellation.
+      internal_xn_prec = prec + (Math.log10(b + l) * batch_size).ceil
+      xn = [BigDecimal(1)]
+      xn << xn.last.mult(x, internal_xn_prec) while xn.size <= batch_size
+
       c = BigDecimal(1)
-      ((b-l)..(b+l)).each do |i|
-        if i != b - l
-          c = c.mult(-b * (b + l - i + 1), prec).div((i - b + l) * i, prec)
+      sum = BigDecimal(0)
+      prod = BigDecimal(1)
+
+      ((b-l)..(b+l)).to_a.each_slice(batch_size) do |batch_ks|
+        # Calculate prod{x-k} in this batch
+        batch_prod, prod_coef = _x_minus_k_prod_coef(batch_ks, xn, internal_xn_prec)
+
+        # Calculate coefficients of batch_prod/(x-k)
+        batch_coef = [0] * batch_ks.size
+        c_scale = 1r
+        batch_ks.each do |k|
+          c_scale = c_scale * (-b * (b + l - k + 1)) / ((k - b + l) * k) if k != b - l
+          rem = 0
+          (batch_ks.size - 1).downto(0) do |i|
+            quo = prod_coef[i + 1] + rem
+            rem = quo * k
+            batch_coef[i] += c_scale * quo
+          end
         end
-        sum = sum.add(c.div(x - i, prec), prec)
+
+        batch_sum = BigDecimal(0)
+        batch_coef.each_with_index do |coef, i|
+          batch_sum = batch_sum.add(xn[i].mult(coef.numerator, internal_xn_prec).div(coef.denominator, internal_xn_prec), internal_xn_prec)
+        end
+        sum = sum.add(batch_sum.mult(c, prec).div(batch_prod, prec), prec)
+        c = c.mult(c_scale.numerator, prec).div(c_scale.denominator, prec)
+        prod = prod.mult(batch_prod, prec)
+      end
+
+      # Perform shift.times {|i| prod = prod.mult(x - i, prec) } with batch processing
+      shift.times.to_a.each_slice(batch_size) do |batch_ks|
+        shift_prod, _prod_coef = _x_minus_k_prod_coef(batch_ks, xn, internal_xn_prec)
+        prod = prod.mult(shift_prod, prec)
       end
     else
       # Binary splitting
+      prods = ((b-l)..(b+l)).map {|i| x - i } + shift.times.map {|i| x - i }
+      prods = prods.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a } while prods.size != 1
+      prod = prods.first
       fractions = (b - l + 1..b + l).map do |i|
         denominator = (x - i).mult((i - b + l) * i, prec)
         numerator = (x - i + 1).mult(-b * (b + l - i + 1), prec)