From 82c96140225aa05ba105fda0923e617ed2e355cd Mon Sep 17 00:00:00 2001 From: tompng Date: Thu, 9 Apr 2026 03:50:35 +0900 Subject: [PATCH 01/11] Lagrange interpolation based gamma/lgamma calculation. Calculates `gamma(x)` by Lagrange interpolation of `b^x/x!` where `b` is `x.round`. Implements Binary Splitting Method version for small-digit number and Baby-Step Giant-Step version for full-digit number. Fallback to Stirling's asymptotic expansion if `x` is extremely large. --- bigdecimal.gemspec | 1 + lib/bigdecimal/math.rb | 120 +------- lib/bigdecimal/math/gamma.rb | 495 ++++++++++++++++++++++++++++++++ test/bigdecimal/test_bigmath.rb | 55 ++-- 4 files changed, 535 insertions(+), 136 deletions(-) create mode 100644 lib/bigdecimal/math/gamma.rb diff --git a/bigdecimal.gemspec b/bigdecimal.gemspec index 38a0f6d6..b5c3c255 100644 --- a/bigdecimal.gemspec +++ b/bigdecimal.gemspec @@ -30,6 +30,7 @@ Gem::Specification.new do |s| lib/bigdecimal/ludcmp.rb lib/bigdecimal/math.rb lib/bigdecimal/math/erf.rb + lib/bigdecimal/math/gamma.rb lib/bigdecimal/newton.rb lib/bigdecimal/util.rb sample/linear.rb diff --git a/lib/bigdecimal/math.rb b/lib/bigdecimal/math.rb index f5754766..af7fb2cf 100644 --- a/lib/bigdecimal/math.rb +++ b/lib/bigdecimal/math.rb @@ -626,22 +626,8 @@ def erfc(x, prec) # #=> "0.17724538509055160272981674833411e1" # def gamma(x, prec) - prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma) - x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma) - prec2 = prec + BigDecimal::Internal::EXTRA_PREC - if x < 0.5 - raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero? - - # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z) - pi = PI(prec2) - sin = _sinpix(x, pi, prec2) - 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) - end - - a, sum = _gamma_spouge_sum_part(x, prec2) - (x + (a - 1)).power(x - 0.5, prec2).mult(BigMath.exp(1 - x, prec2), prec2).mult(sum, prec) + require 'bigdecimal/math/gamma' + Gamma.gamma(x, prec) end # call-seq: @@ -654,106 +640,8 @@ def gamma(x, prec) # #=> [0.57236494292470008707171367567653e0, 1] # def lgamma(x, prec) - prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma) - x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma) - prec2 = prec + BigDecimal::Internal::EXTRA_PREC - if x < 0.5 - return [BigDecimal::INFINITY, 1] if x.frac.zero? - - loop do - # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z) - pi = PI(prec2) - sin = _sinpix(x, pi, prec2) - log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec) - return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC - - # Retry with higher precision if loss of significance is too large - prec2 = prec2 * 3 / 2 - end - elsif x.frac.zero? && x < 1000 * prec - log_gamma = BigMath.log(_gamma_positive_integer(x, prec2), prec) - [log_gamma, 1] - else - # if x is close to 1 or 2, increase precision to reduce loss of significance - diff1_exponent = (x - 1).exponent - diff2_exponent = (x - 2).exponent - extremely_near_one = diff1_exponent < -prec2 - extremely_near_two = diff2_exponent < -prec2 - - if extremely_near_one || extremely_near_two - # If x is extreamely close to base = 1 or 2, linear interpolation is accurate enough. - # Taylor expansion at x = base is: (x - base) * digamma(base) + (x - base) ** 2 * trigamma(base) / 2 + ... - # And we can ignore (x - base) ** 2 and higher order terms. - base = extremely_near_one ? 1 : 2 - d = BigDecimal(1)._decimal_shift(1 - prec2) - log_gamma_d, sign = lgamma(base + d, prec2) - return [log_gamma_d.mult(x - base, prec2).div(d, prec), sign] - 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) - [log_gamma, 1] - end - end - - # Returns sum part: sqrt(2*pi) and c[k]/(x+k) terms of Spouge's approximation - private_class_method def _gamma_spouge_sum_part(x, prec) # :nodoc: - x -= 1 - # Spouge's approximation - # x! = (x + a)**(x + 0.5) * exp(-x - a) * (sqrt(2 * pi) + (1..a - 1).sum{|k| c[k] / (x + k) } + epsilon) - # where c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(a - k) / (k - 1)! - # and epsilon is bounded by a**(-0.5) * (2 * pi) ** (-a - 0.5) - - # Estimate required a for given precision - a = (prec / Math.log10(2 * Math::PI)).ceil - - # Calculate exponent of c[k] in low precision to estimate required precision - low_prec = 16 - log10f = Math.log(10) - x_low_prec = x.mult(1, low_prec) - loggamma_k = 0 - ck_exponents = (1..a-1).map do |k| - loggamma_k += Math.log10(k - 1) if k > 1 - -loggamma_k - k / log10f + (k - 0.5) * Math.log10(a - k) - BigDecimal::Internal.float_log(x_low_prec.add(k, low_prec)) / log10f - end - - # Estimate exponent of sum by Stirling's approximation - approx_sum_exponent = x < 1 ? -Math.log10(a) / 2 : Math.log10(2 * Math::PI) / 2 + x_low_prec.add(0.5, low_prec) * Math.log10(x_low_prec / x_low_prec.add(a, low_prec)) - - # Determine required precision of c[k] - prec2 = [ck_exponents.max.ceil - approx_sum_exponent.floor, 0].max + prec - - einv = BigMath.exp(-1, prec2) - sum = (PI(prec) * 2).sqrt(prec).mult(BigMath.exp(-a, prec), prec) - y = BigDecimal(1) - (1..a - 1).each do |k| - # c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(-k) / (k-1)! / (x + k) - y = y.div(1 - k, prec2) if k > 1 - y = y.mult(einv, prec2) - z = y.mult(BigDecimal((a - k) ** k), prec2).div(BigDecimal(a - k).sqrt(prec2).mult(x.add(k, prec2), prec2), prec2) - # sum += c[k] / (x + k) - sum = sum.add(z, prec2) - end - [a, sum] - end - - private_class_method def _gamma_positive_integer(x, prec) # :nodoc: - return x if x == 1 - numbers = (1..x - 1).map {|i| BigDecimal(i) } - while numbers.size > 1 - numbers = numbers.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a } - end - numbers.first - end - - # Returns sin(pi * x), for gamma reflection formula calculation - private_class_method def _sinpix(x, pi, prec) # :nodoc: - x = x % 2 - sign = x > 1 ? -1 : 1 - x %= 1 - x = 1 - x if x > 0.5 # to avoid sin(pi*x) loss of precision for x close to 1 - sign * sin(x.mult(pi, prec), prec) + require 'bigdecimal/math/gamma' + Gamma.lgamma(x, prec) end # call-seq: diff --git a/lib/bigdecimal/math/gamma.rb b/lib/bigdecimal/math/gamma.rb new file mode 100644 index 00000000..f2de143d --- /dev/null +++ b/lib/bigdecimal/math/gamma.rb @@ -0,0 +1,495 @@ +# frozen_string_literal: true +require 'bigdecimal/math' + +module BigMath + + # Calculates gamma/lgamma + # Algorithm overview: + # + # Lagrange interpolation of f(x) = b**x / x! at integer nodes x_i = b-l, ..., b+l. + # BSM(Binary Splitting Method) version for small digit numbers, O(PREC*log(PREC)^3) + # BSGS(Baby-Step Giant-Step) version for full digit numbers, O(PREC^2*log(log(PREC))) + # Both orders of magnitude faster than Spouge's approximation which is O(PREC^2*log(PREC)) + # (Complexities assume quasi-linear multiplication, counting large-by-small products + # as (n/m) * M(m) = n * log(m) bit ops. BigDecimal multiplies the small coefficients + # by schoolbook instead: an extra log factor asymptotically, but faster at any feasible PREC.) + # Requires fast calculation of factorial(nearly_x_integer). + # + # Factorial Doubling for fast calculation of large factorials: + # Using Legendre duplication formula, we can calculate factorial(2n) from factorial(n) and factorial(n + 0.5). + # Calculating factorial(n + 0.5) is done by the BSM version of Lagrange interpolation in quasi-linear time. + # This will drastically reduce the cost of calculating large factorials. + # O(PREC*log(PREC)^3*log(factorial_argument)) + # + # Stirling's approximation with Bernoulli numbers + # Only used when x is extremely large. + + module Gamma # :nodoc: + + # Calculates gamma function with given precision. + def self.gamma(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma) + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + + if x < 0.5 + raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero? + + # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z) + pi = BigMath::PI(prec2) + sin = sinpix(x, pi, prec2) + pi.div(gamma(1 - x, prec2).mult(sin, prec2), prec) + else + # Digits of x beyond the working precision cannot affect the result. + # Rounding must happen before the integer test: an x indistinguishable from + # an integer must take the exact integer path, because gamma_lagrange + # requires a non-integer x (an integer x makes a node distance exactly zero). + x = x.mult(1, prec2 + x.exponent + 10) + if x.frac.zero? + integer_factorial(x.to_i - 1, prec2).mult(1, prec) + else + base, large_factorial_arg, small_factorial_arg, exp2 = gamma_lagrange(x, prec2) + ans = base.mult(integer_factorial(small_factorial_arg, prec2), prec2) + ans = ans.mult(BigDecimal(2).power(exp2, prec2), prec2) unless exp2.zero? + ans.mult(integer_factorial(large_factorial_arg, prec2), prec) + end + end + end + + # Calculates log gamma and its sign with given precision. + def self.lgamma(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma) + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + if x < 0.5 + return [BigDecimal::INFINITY, 1] if x.frac.zero? + + loop do + # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z) + pi = BigMath::PI(prec2) + sin = sinpix(x, pi, prec2) + log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec) + return [log_gamma, sin > 0 ? 1 : -1] if log_gamma != 0 && prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC + + # Retry with higher precision if loss of significance is too large + prec2 = prec2 * 3 / 2 + end + else + # if x is close to 1 or 2, increase precision to reduce loss of significance + diff1_exponent = x < 3 ? (x - 1).exponent : 0 + diff2_exponent = x < 3 ? (x - 2).exponent : 0 + extremely_near_one = diff1_exponent < -prec2 + extremely_near_two = diff2_exponent < -prec2 + + if extremely_near_one || extremely_near_two + # If x is extremely close to base = 1 or 2, linear interpolation is accurate enough. + # Taylor expansion at x = base is: (x - base) * digamma(base) + (x - base) ** 2 * trigamma(base) / 2 + ... + # And we can ignore (x - base) ** 2 and higher order terms. + base = extremely_near_one ? 1 : 2 + d = BigDecimal(1)._decimal_shift(1 - prec2) + log_gamma_d, sign = lgamma(base + d, prec2) + return [log_gamma_d.mult(x - base, prec2).div(d, prec), sign] + end + + prec2 += [-diff1_exponent, -diff2_exponent, 0].max + + # Same rounding as in Gamma.gamma. Must come after the near 1 and 2 handling + # above, which needs the exact distance from x to 1 or 2. + x = x.mult(1, prec2 + x.exponent + 10) + + # When x is extremely large, the cost of Bernoulli number generation for Stirling's + # asymptotic expansion is smaller than the cost of multiple steps of doubling method. + # The condition is based on heuristic cost estimation and empirical tuning. + if x > prec2 && x.exponent > Integer.sqrt(prec2) / 6 + [lgamma_stirling(x, prec2).mult(1, prec), 1] + elsif x.frac.zero? + [integer_factorial_log(x.to_i - 1, prec2).mult(1, prec), 1] + else + base, large_factorial_arg, small_factorial_arg, exp2 = gamma_lagrange(x, prec2) + lgamma = BigMath.log(base, prec2) + lgamma = lgamma.add(BigMath.log(2, prec2) * exp2, prec2) unless exp2.zero? + lgamma = lgamma.add(integer_factorial_log(small_factorial_arg, prec2), prec2) + lgamma = lgamma.add(integer_factorial_log(large_factorial_arg, prec2), prec) + [lgamma, 1] + end + end + 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, ...] + def self.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 + + prod = coef.each_with_index.map do |c, i| + xn[i].mult(c, prec) + end.reduce do |sum, value| + sum.add(value, prec) + end + [prod, coef] + end + + # Bit length to keep in bit-dropping integer products and fraction merges. + # Only about prec * log2(10) bits are needed. 10 / 3 slightly exceeds log2(10), + # and 64 extra bits absorb the ~1 bit lost per truncation over the tree depth. + def self.drop_cap_bits(prec) + (prec * 10 + 192) / 3 + end + + # Calculate numbers.reduce(:*) of integers by Binary Splitting Method. + # Returns [mantissa, exp2] representing mantissa * 2**exp2. + # With cap, lower bits of intermediate products are dropped to keep each + # multiplication cost bounded by the target precision. Without cap, the product is exact. + def self.int_bsm_prod(numbers, cap = nil) + numbers = numbers.to_a + exp2 = 0 + while numbers.size > 1 + numbers = numbers.each_slice(2).map do |a, b| + next a unless b + v = a * b + if cap && (s = v.bit_length - cap) > 0 + exp2 += s + v >>= s + end + v + end + end + [numbers.first || 1, exp2] + end + + # Calculate factorial for integer n + def self.integer_factorial(n, prec) + power_part, exp2, exp_sqrtpi = integer_factorial_parameter(n, prec) + ans = BigDecimal(2).power(exp2, prec) + power_part.each_with_index do |base, index| + ans = ans.mult(base.power(1 << index, prec), prec) + end + if exp_sqrtpi != 0 + pi = BigMath::PI(prec) + # exp_sqrtpi is 2**k - 1 (odd): the doubling recursion squares the child's + # sqrt(pi) exponent and adds one, so only the last level's sqrt survives + pipow = pi.power(exp_sqrtpi / 2, prec).mult(pi.sqrt(prec), prec) + ans = ans.div(pipow, prec) + end + ans + end + + # Calculate log factorial for integer n + def self.integer_factorial_log(n, prec) + power_part, exp2, exp_sqrtpi = integer_factorial_parameter(n, prec) + ans = exp2.zero? ? BigDecimal(0) : BigMath.log(2, prec) * exp2 + power_part.each_with_index do |base, index| + ans = ans.add(BigMath.log(base, prec) * (1 << index), prec) + end + if exp_sqrtpi != 0 + pi = BigMath::PI(prec) + ans = ans.sub(BigMath.log(pi, prec) * (BigDecimal(exp_sqrtpi) / 2), prec) + end + ans + end + + # Calculates parameters for integer factorial calculation. + # Returns [base_power_part, exp2, exp_sqrtpi] that can produce factorial(n) as: + # factorial(n) = prod { base_power_part[i]**(1 << i) } * 2**exp2 / sqrt(pi)**exp_sqrtpi + # These parameters are used to avoid overflow when calculating log factorial and lgamma for large n. + def self.integer_factorial_parameter(n, prec) + base_power_part, factorial_power_part, exp2, exp_sqrtpi = integer_factorial_recursive(n, prec) + fact_x = 1 + fact_y = BigDecimal(1) + # factorial_power_part is non-decreasing (deeper recursion levels have smaller b, + # and gamma_lagrange_l grows as b shrinks), so fact_y can be extended incrementally. + factorial_power_part.each_with_index do |factorial_arg, index| + # Exact product (no bit drop): these ranges total only O(prec * log(prec)) digits, + # and a dropped 2**s here would be raised to 2**index, exceeding the representable + # exponent range while base_power_part[index] underflows by the same amount. + mantissa, = int_bsm_prod(fact_x + 1..factorial_arg) + fact_y = fact_y.mult(mantissa, prec) + fact_x = factorial_arg + base_power_part[index] = base_power_part[index].mult(fact_y, prec) + end + [base_power_part, exp2, exp_sqrtpi] + end + + # Returns [base_power_part, factorial_power_part, exp2, exp_sqrtpi] that can produce factorial(n) as: + # factorial(n) = prod { base_power_part[i]**(1 << i) } * prod { factorial(factorial_power_part[i])**(1 << i) } * 2**exp2 / sqrt(pi)**(exp_sqrtpi) + # If n is large, this method recursively calculates factorial for smaller n by Legendre duplication formula. + def self.integer_factorial_recursive(n, prec) + if n < 4 * prec + mantissa, exp2 = int_bsm_prod(1..n, drop_cap_bits(prec)) + return [[BigDecimal(mantissa)], [], exp2, 0] + end + + # Use Legendre duplication formula to calculate double factorials: + # factorial(n) = factorial(n/2.0) * factorial((n-1)/2.0) * 2**n / sqrt(pi) + # gamma_lagrange((n + 1) / 2 + 0.5, prec) computes the half-integer factorial + # (whichever of the two factors above is a half-integer). + half_arg = BigDecimal((n + 1) / 2) + BigDecimal('0.5') + base, large_factorial_arg, small_factorial_arg, lagrange_exp2 = gamma_lagrange(half_arg, prec) + + range_mantissa, = int_bsm_prod(large_factorial_arg + 1..n / 2) # exact: total size is O(prec) digits + base = base.mult(range_mantissa, prec) + base_power_part, factorial_power_part, exp2, exp_sqrtpi = integer_factorial_recursive(large_factorial_arg, prec) + [ + [base] + base_power_part, + [small_factorial_arg] + factorial_power_part, + exp2 * 2 + n + lagrange_exp2, + exp_sqrtpi * 2 + 1 + ] + end + + # Estimate the required number of interpolation points `l` to achieve `prec` digits. + # + # Assuming the nodes stay strictly positive, the function b^x/x! approximates a + # Gaussian curve e^(-y^2 / 2b) around its peak (x=b). + # The Taylor coefficient of degree 2l is roughly c_2l = 1 / (l! * (2b)^l). + # Multiplying this by the distance product of 2l+1 nodes (approx (l/e)^(2l)), + # the overall truncation error E is bounded by: E ~ (l / 2eb)^l. + # + # Setting E <= 10^-prec gives the implicit equation: + # l * log10(2 * e * b / l) = prec => l = prec / log10(2 * e * b / l) + def self.gamma_lagrange_l(b, prec) + # Initial guess of l. When b >= 2 * prec - 1 (guaranteed by the shift in gamma_lagrange), + # this is safely larger than the actual l. + l = prec + + # Solves the implicit equation via fixed-point iteration. + # Due to the slow growth of the logarithm, 2 iterations are practically sufficient. + 2.times { l = prec / Math.log10(2 * Math::E * b / l) } + l.ceil + 10 # Adds safety margin + end + + # Calculate approximate gamma by Lagrange interpolation of f(x) = b**x / x! + # Nodes are placed at x_i = b-l, b-l+1, ..., b+l. + # b: x.round, l: number of nodes on one side (total nodes = 2*l+1) + # + # Mathematically, we use the barycentric interpolation form: + # f(x) \approx \omega(x) \sum_{i} \frac{w_i f(x_i)}{x - x_i} + # Therefore, \Gamma(x+1) = x! = b**x / f(x) + # + # Time complexity: + # - O(PREC*log(PREC)^3) for small-digit x (Binary Splitting) + # - O(PREC^2) for full-digit x (Baby-step Giant-step) + # + # Returns [base, large_factorial_arg, small_factorial_arg, exp2] that can produce gamma(x) as: + # gamma(x) = base * 2**exp2 * factorial(large_factorial_arg) * factorial(small_factorial_arg) + def self.gamma_lagrange(x, prec) + # Shift x to establish a safe center (b) for the barycentric interpolation. + # + # We must keep all interpolation nodes strictly positive (b - l > 0). Approaching + # x = 0 breaks the Gaussian approximation used to estimate `l` and provides no + # useful information for the interpolation. + # + # While b =~ 1.36 * prec is the strict theoretical minimum to stay positive, we + # heuristically use b = 2 * prec. This moves the nodes safely away from x = 0, + # stabilizes the curve, and empirically yields the optimal total computation cost. + # (See `gamma_lagrange_l` for the mathematical derivation of the approximation). + shift = x < 2 * prec ? 2 * prec - x.floor : 0 + x += shift + + x = BigDecimal(x) - 1 + b = x.round + l = gamma_lagrange_l(b, prec) + exp2 = 0 + + # --- Reference: Naive interpolation logic --- + # 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) # represents w_i * f(x_i) (normalized) + # (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 + # -------------------------------------------- + + # Choose between BSM and BSGS based on total bit cost: + # BSM: (l * n_sig / prec) full-digit multiplications, each costing prec * log(prec) + # bit ops, total l * n_sig * log(prec). + # BSGS: l * prec bit ops with batch_size = log2(prec) (see below). + # Cross-over: n_sig * log(prec) > prec. + if x.n_significant_digits * prec.bit_length > prec + # Reduce full-precision multiplications/divisions using a Batched Evaluation + # inspired by the Baby-Step Giant-Step (BSGS) method. + + # Normal BSGS uses batch_size = sqrt(l), but here the integer coefficients of the + # expanded prod { x - k } over a batch grow like (b + l)**batch_size, so a smaller + # batch keeps both the coefficient size and the per-batch evaluation cost low. + # batch_size = log2(prec) brings the total BSGS bit cost down to O(l * prec * log(log(prec))). + batch_size = prec.bit_length + + # When expanding prod { x - k }, the coefficient of x**n might be huge. + # Increase internal calculation precision to avoid catastrophic cancellation. + # When x is within 10**-q of a node, batch_prod cancels by q more digits; + # without the extra digits the computed batch_prod can even collapse to + # exactly zero, and the division below would raise or produce NaN. + nearest_node_distance = x - x.round(0, BigDecimal::ROUND_HALF_UP) + near_node_digits = [1 - nearest_node_distance.exponent, 0].max + internal_xn_prec = prec + (Math.log10(b + l) * batch_size).ceil + near_node_digits + xn = [BigDecimal(1)] + xn << xn.last.mult(x, internal_xn_prec) while xn.size <= batch_size + + c = BigDecimal(1) + 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) using Synthetic Division (Ruffini's rule) + 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 + + 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 + # batch_prod loses relative accuracy when x is extremely close to a node in this + # batch. This is harmless: the same computed value is divided into sum here and + # multiplied into prod below, so the error cancels in the final prod * sum. + 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 Method (BSM) for short-digit inputs. + # Scaling x by 10**frac_digits makes every node term x - i an exact integer, + # so the whole tree runs on Integer arithmetic. The scale cancels inside the + # fraction merges; only prod and sum need explicit rescaling. + frac_digits = [x.n_significant_digits - x.exponent, 0].max + s10 = 10**frac_digits + xs = x._decimal_shift(frac_digits).to_i + cap = drop_cap_bits(prec) + + prod_factors = (b - l..b + l).map {|i| xs - i * s10 } + shift.times.map {|i| xs - i * s10 } + mantissa, dropped_exp2 = int_bsm_prod(prod_factors, cap) + prod = BigDecimal(mantissa)._decimal_shift(-frac_digits * prod_factors.size) + # prod is missing the dropped 2**dropped_exp2 factor and gamma is proportional + # to 1 / prod. Returning the compensation as exp2 lets the factorial doubling + # fold it into its own 2**exp2 channel instead of paying a power here. + exp2 = -dropped_exp2 + + # State represents a partial evaluation of the series as: [sum_num, mult_num, den] + # Conceptually, each state translates to the following mathematical expression: + # (sum_num / den) + (mult_num / den) * (rest_of_the_series) + # + # The initial state [denominator, numerator, denominator] simply represents: + # (denominator / denominator) + (numerator / denominator) * rest + # = 1 + (numerator / denominator) * rest + # + fractions = (b - l + 1..b + l).map do |i| + denominator = (xs - i * s10) * ((i - b + l) * i) + numerator = (xs - (i - 1) * s10) * (-b * (b + l - i + 1)) + [denominator, numerator, denominator] + end + + while fractions.size > 1 + fractions = fractions.each_slice(2).map do |a, c| + c ||= [1, 0, 1] + # Merge operation for BSM: + # a[0]/a[2] + a[1]/a[2] * (c[0]/c[2] + c[1]/c[2] * rest) + # = (a[0]*c[2] + a[1]*c[0]) / (a[2]*c[2]) + (a[1]*c[1]) / (a[2]*c[2]) * rest + v0 = a[0] * c[2] + a[1] * c[0] + v1 = a[1] * c[1] + v2 = a[2] * c[2] + # Drop lower bits to avoid the integers growing too large; see drop_cap_bits. + # All three components share the shift, so the represented ratios are unchanged. + s = v2.bit_length - cap + if s > 0 + v0 >>= s + v1 >>= s + v2 >>= s + end + [v0, v1, v2] + end + end + fraction = fractions.first + sum = BigDecimal((fraction[0] + fraction[1]) * s10).div(fraction[2] * (xs - (b - l) * s10), prec) + end + + # Reconstruct Gamma(x_original) by reversing the scaling and applying shift formula + e = x - (b - l) + if e.frac == 0.5 + # The factorial doubling path always has x = integer + 0.5. For this exponent shape, + # integer power and sqrt are much cheaper than the exp/log based fractional power. + power_part = BigDecimal(b).power(e.to_i, prec).mult(BigDecimal(b).sqrt(prec), prec) + else + power_part = BigDecimal(b).power(e, prec) + end + base = power_part.div(prod.mult(sum, prec), prec) + large_factorial_arg = b - l + small_factorial_arg = 2 * l + [base, large_factorial_arg, small_factorial_arg, exp2] + end + + # Calculates bernoulli number. + # bns: calculated bernoulli numbers for memoization + def self.bernoulli(n, bns, prec) + return bns[0] ||= BigDecimal(1) if n == 0 + return bns[1] ||= BigDecimal(-0.5) if n == 1 + return bns[n] ||= BigDecimal(0) if n.odd? + bns[n] ||= ( + comb = 1 + sum = BigDecimal(0) + n.times do |i| + sum = sum.add(comb * bernoulli(i, bns, prec), prec) + comb = comb * (n - i + 1) / (i + 1) + end + sum.div(-n - 1, prec) + ) + end + + # Calculate gamma using Stirling's asymptotic expansion. + # While the condition of this asymptotic expansion is x > prec * log(10) / 2 / pi, + # we'll use this method only when x is extremely large to reduce the cost of Bernoulli number generation. + def self.lgamma_stirling(x, prec) + x = BigDecimal(x) + y = (x * (BigMath.log(x, prec) - 1)).add(BigMath.log(2 * BigMath::PI(prec).div(x, prec), prec) / 2, prec) + bns = [] + xn = x + x2 = x.mult(x, prec) + (1..).each do |k| + xn = xn.mult(x2, prec) if k != 1 + d = bernoulli(2 * k, bns, prec).div(xn, prec).div(2 * k * (2 * k - 1), prec) + y = y.add(d, prec) + break if d.exponent < y.exponent - prec + end + y + end + + # Returns sin(pi * x), for gamma reflection formula calculation + def self.sinpix(x, pi, prec) + x = x % 2 + sign = x > 1 ? -1 : 1 + x %= 1 + x = 1 - x if x > 0.5 # to avoid sin(pi*x) loss of precision for x close to 1 + sign * BigMath.sin(x.mult(pi, prec), prec) + end + end + + private_constant :Gamma +end diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index 05bcbc82..f9891615 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -575,11 +575,20 @@ 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) } + 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) } + assert_converge_in_precision {|n| gamma(BigDecimal("1e8"), n) } + assert_converge_in_precision {|n| gamma(BigDecimal("1e15") + BigDecimal("0.5"), n) } + assert_converge_in_precision {|n| gamma(BigDecimal(1).div(3, n * 2), n) } + assert_converge_in_precision {|n| gamma(10000000 + BigDecimal(1).div(3, n * 2), n) } + # x extremely close to an interpolation node stresses the batch_prod cancellation in BSGS + assert_converge_in_precision {|n| gamma(BigDecimal(5) + BigDecimal("1e-40"), n) } + # x closer to a node than the working precision must not produce a zero node distance + assert_converge_in_precision {|n| gamma(BigDecimal(5) + BigDecimal("1e-300"), n) } + assert_converge_in_precision {|n| gamma(BigDecimal("-3") - BigDecimal("1e-2000"), n) } + assert_equal(BigDecimal(24), gamma(BigDecimal(5) + BigDecimal("1e-2000"), 50)) end def test_lgamma @@ -595,23 +604,29 @@ 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] - 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 } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("2." + "0" * 80 + "1"), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("-1." + "9" * 30), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("-3." + "0" * 30 + "1"), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("10"), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("0.3"), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("-1.9" + "9" * 30), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("987.65421"), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("-1234.56789"), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal("1e+400"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("0." + "9" * 80), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("1." + "0" * 80 + "1"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("1." + "9" * 80), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("2." + "0" * 80 + "1"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("-1." + "9" * 30), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("-3." + "0" * 30 + "1"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("10"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("0.3"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("-1.9" + "9" * 30), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("987.65421"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("-1234.56789"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("1e+18"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal("1e+400"), n).first } # gamma close 1 or -1 cases - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal('-3.143580888349980058694358781820227899566'), n).first } - assert_converge_in_precision(precisions) {|n| lgamma(BigDecimal('-4.991544640560047722345260122806465721667'), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal('-3.143580888349980058694358781820227899566'), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal('-4.991544640560047722345260122806465721667'), n).first } + + # x closer to a node than the working precision must not produce a zero node distance + assert_converge_in_precision {|n| lgamma(BigDecimal(5) + BigDecimal("1e-300"), n).first } + # x closer to 1 or 2 than any tested precision takes the linear interpolation path + assert_converge_in_precision {|n| lgamma(BigDecimal(1) + BigDecimal("1e-2000"), n).first } + assert_converge_in_precision {|n| lgamma(BigDecimal(2) + BigDecimal("1e-2000"), n).first } end def test_frexp From bf1fe2fe118eb7dc75558be11db8b2d12f3ac2ef Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 00:17:26 +0900 Subject: [PATCH 02/11] Experimental multipoint evaluation of Lagrange interpolation for gamma Lift the BSM [sum_num, mult_num, den] triple to polynomials in the batch offset z, build them with a product tree over fixed-point coefficients via Kronecker substitution onto Integer (GMP) multiplication, and evaluate at the arithmetic progression z = 0, m, 2m, ... Polynomial work is O(PREC^1.5 * polylog); evaluation is currently per-point Horner. Matches the BSGS implementation exactly up to 50000 digits in tests. Crossover is around 10000 digits (1.9x faster at 50000 digits). Co-Authored-By: Claude Fable 5 --- gamma_mp_check.rb | 54 ++++++ lib/bigdecimal/math/gamma_multipoint.rb | 248 ++++++++++++++++++++++++ 2 files changed, 302 insertions(+) create mode 100644 gamma_mp_check.rb create mode 100644 lib/bigdecimal/math/gamma_multipoint.rb diff --git a/gamma_mp_check.rb b/gamma_mp_check.rb new file mode 100644 index 00000000..b6414b7d --- /dev/null +++ b/gamma_mp_check.rb @@ -0,0 +1,54 @@ +# Check & benchmark for the experimental multipoint gamma (gamma_multipoint.rb) +# Usage: ruby -Ilib -Itmp/arm64-darwin24/stage/lib gamma_mp_check.rb [mode] +# mode: acc (default) | bench | debug +require 'bigdecimal' +require 'bigdecimal/math' +require 'bigdecimal/math/gamma_multipoint' +require 'benchmark' + +MP = BigMath.const_get(:Gamma)::Multipoint +G = BigMath.const_get(:Gamma) + +def rel_err_exp(a, b, prec) + e = a.sub(b, prec + 50).div(b, 10).abs + e.zero? ? :exact : e.exponent +end + +mode = ARGV[0] || 'acc' + +case mode +when 'debug' + # Tiny case: compare mp against the regular implementation step by step + prec = 50 + x = BigDecimal(2).sqrt(150) + a = MP.gamma(x, prec) + b = BigMath.gamma(x, prec + 20) + puts "mp = #{a.to_s("F")[0, 60]}" + puts "ref = #{b.to_s("F")[0, 60]}" + puts "rel_err_exp = #{rel_err_exp(a, b, prec)}" +when 'acc' + [100, 200, 500, 1000, 2000].each do |prec| + cases = { + "sqrt2" => BigDecimal(2).sqrt(2 * prec + 50), + "1/3" => BigDecimal(1).div(3, 2 * prec + 50), + "near-node 7+eps" => BigDecimal(7) + BigDecimal(1).div(3, prec + 50)._decimal_shift(-(prec / 2)), + "0.6" => BigDecimal("0.6") + BigDecimal(1).div(7, 2 * prec + 50)._decimal_shift(-3), + } + cases.each do |name, x| + t_mp = Benchmark.realtime { @mp = MP.gamma(x, prec) } + ref = BigMath.gamma(x, prec + 50) + e = rel_err_exp(@mp, ref, prec) + ok = e == :exact || e <= -prec + puts format("%s prec=%-5d %-16s rel_err_exp=%-6s mp=%.2fs", ok ? "OK " : "FAIL", prec, name, e, t_mp) + end + end +when 'bench' + [2000, 5000, 10000].each do |prec| + x = BigDecimal(2).sqrt(2 * prec + 50) + t_mp = Benchmark.realtime { @mp = MP.gamma(x, prec) } + t_ref = Benchmark.realtime { @ref = BigMath.gamma(x, prec) } + refhi = BigMath.gamma(x, prec + 50) + puts format("prec=%-6d mp=%.2fs bsgs=%.2fs (%.1fx) mp_err=%s bsgs_err=%s", + prec, t_mp, t_ref, t_ref / t_mp, rel_err_exp(@mp, refhi, prec), rel_err_exp(@ref, refhi, prec)) + end +end diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb new file mode 100644 index 00000000..ca6df3c7 --- /dev/null +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -0,0 +1,248 @@ +# frozen_string_literal: true + +# Experimental multipoint-evaluation version of the Lagrange interpolation +# used by BigMath.gamma, targeting full-digit x. +# +# The BSGS version in gamma.rb costs O(PREC^2 * polylog): every node needs a +# scalar multiplication against a full-precision power of x. This file evaluates +# the same barycentric sum with sqrt-size batches instead: +# - The [sum_num, mult_num, den] triple of the BSM branch is lifted to +# polynomials in the batch offset z. One triple tree describes all batches. +# - Polynomial arithmetic runs on fixed-point coefficients via Kronecker +# substitution onto Integer multiplication, so it needs quasi-linear Integer +# multiplication (GMP-backed Ruby). +# - The polynomials are evaluated at the arithmetic progression z = 0, mb, +# 2*mb, ... (currently by per-point Horner with small multipliers; a fast +# Newton-basis transform can replace it later). +# Polynomial work is O(PREC^1.5 * polylog). +# +# The batch denominator values E(z) used in prod are derived from the same +# computed F2(z) used in sum (E = F2 * (x-A-z) / (B*I) with exact integer B, I), +# so the near-node cancellation between prod and sum stays exact, like the +# batch_prod reuse in the BSGS branch. + +require 'bigdecimal/math/gamma' + +module BigMath + module Gamma + module Multipoint # :nodoc: + + # ---------- Kronecker substitution convolution on Integer ---------- + + def self.pack(coeffs, slot_hex) + coeffs.reverse_each.map {|c| c.to_s(16).rjust(slot_hex, '0') }.join.to_i(16) + end + + def self.pack_signed(coeffs, slot_hex) + v = pack(coeffs.map {|c| c > 0 ? c : 0 }, slot_hex) + v -= pack(coeffs.map {|c| c < 0 ? -c : 0 }, slot_hex) if coeffs.any? {|c| c < 0 } + v + end + + def self.unpack_signed(n, slot_hex, size) + half = 1 << (slot_hex * 4 - 1) + bias = (('8' + '0' * (slot_hex - 1)) * size).to_i(16) + s = (n + bias).to_s(16).rjust(slot_hex * size, '0') + (0...size).map {|i| s[(size - 1 - i) * slot_hex, slot_hex].to_i(16) - half } + end + + # Convolution of signed Integer coefficient arrays. + def self.convolve(a, b) + out_size = a.size + b.size - 1 + if a.size < 16 || b.size < 16 + out = Array.new(out_size, 0) + a.each_with_index {|c, i| b.each_with_index {|d, j| out[i + j] += c * d } } + return out + end + max_bits = 1 + a.each {|c| bits = c.abs.bit_length; max_bits = bits if bits > max_bits } + b.each {|c| bits = c.abs.bit_length; max_bits = bits if bits > max_bits } + w = 2 * max_bits + out_size.bit_length + 2 + slot_hex = (w + 3) / 4 + prod = pack_signed(a, slot_hex) * pack_signed(b, slot_hex) + unpack_signed(prod, slot_hex, out_size) + end + + # ---------- fixed-point polynomials ---------- + # Represented as [coeffs, exp]: sum of coeffs[d] * 2**exp * z**d. + # A single exp per polynomial (fixed-point): small coefficients keep less + # relative precision, which only affects small contributions to the value. + + def self.fp_normalize(coeffs, exp, keep_bits) + max = 0 + coeffs.each {|c| bits = c.abs.bit_length; max = bits if bits > max } + s = max - keep_bits + return [coeffs, exp] if s <= 0 + [coeffs.map {|c| c >> s }, exp + s] + end + + def self.fp_mult(p1, p2, keep_bits) + fp_normalize(convolve(p1[0], p2[0]), p1[1] + p2[1], keep_bits) + end + + def self.fp_add(p1, p2, keep_bits) + c1, e1 = p1 + c2, e2 = p2 + if e1 > e2 + c2 = c2.map {|c| c >> (e1 - e2) } + e = e1 + elsif e2 > e1 + c1 = c1.map {|c| c >> (e2 - e1) } + e = e2 + else + e = e1 + end + out = Array.new(c1.size > c2.size ? c1.size : c2.size, 0) + c1.each_with_index {|c, i| out[i] += c } + c2.each_with_index {|c, i| out[i] += c } + fp_normalize(out, e, keep_bits) + end + + # Merge of [sum_num, mult_num, den] triples, same as the BSM merge in + # gamma.rb but over polynomials. + def self.triple_merge(a, c, keep_bits) + [ + fp_add(fp_mult(a[0], c[2], keep_bits), fp_mult(a[1], c[0], keep_bits), keep_bits), + fp_mult(a[1], c[1], keep_bits), + fp_mult(a[2], c[2], keep_bits) + ] + end + + # Exact Horner evaluation at an integer point. Returns the Integer mantissa; + # the value is mantissa * 2**poly_exp. + def self.fp_eval_int(poly, z) + acc = 0 + poly[0].reverse_each {|c| acc = acc * z + c } + acc + end + + # Guard bits on top of the target precision, absorbing: + # - coefficient spread and value dynamic range across batches (~m * log2(n1)) + # - rounding of ~log2(m) tree levels and of the evaluation + # Deliberately generous; to be tightened after error measurements. + def self.guard_bits(m, n1) + 4 * m * (n1.bit_length + 4) + 256 + end + + # Same contract as Gamma.gamma_lagrange. + # Interpolation nodes are A .. A + n1 - 1 with A = b - l and n1 = m**2 + # (m odd so that the barycentric reconstruction keeps positive sign); + # slightly wider than the symmetric b-l .. b+l, which only adds accuracy. + def self.gamma_lagrange(x, prec) # :nodoc: + shift = x < 2 * prec ? 2 * prec - x.floor : 0 + x += shift + x = BigDecimal(x) - 1 + b = x.round + l = Gamma.gamma_lagrange_l(b, prec) + + m = Integer.sqrt(2 * l) + 1 + m += 1 if m.even? + n1 = m * m + a0 = b - l + + keep = Gamma.drop_cap_bits(prec) + guard_bits(m, n1) + s2 = 1 << keep + + # Fixed-point mantissas (keep fractional bits) of x - a0 and x + fd = [x.n_significant_digits - x.exponent, 0].max + p10 = 10**fd + xa = ((x - a0)._decimal_shift(fd).to_i << keep) / p10 + + # Triple tree over leaves t = z + j (j = 1 .. m-1), as polynomials in z: + # den_t = (x - a0 - t) * (t * (a0 + t)) + # num_t = (x - a0 - t + 1) * (-b * (n1 - t)) + identity = [[[1], 0], [[0], 0], [[1], 0]] + fractions = (1..m - 1).map do |j| + xaj = xa - j * s2 + den = fp_normalize( + [xaj * (j * (a0 + j)), xaj * (a0 + 2 * j) - s2 * (j * (a0 + j)), xaj - s2 * (a0 + 2 * j), -s2], + -keep, keep + ) + xaj1 = xaj + s2 + num = fp_normalize( + [-b * xaj1 * (n1 - j), b * (xaj1 + s2 * (n1 - j)), -b * s2], + -keep, keep + ) + [den, num, den] + end + while fractions.size > 1 + fractions = fractions.each_slice(2).map do |p, q| + q ||= identity + triple_merge(p, q, keep) + end + end + f0, f1, f2 = fractions.first + f01 = fp_add(f0, f1, keep) + e01 = f01[1] + e2 = f2[1] + + sum = BigDecimal(0) + prod = BigDecimal(1) + c_k = BigDecimal(1) + m.times do |k| + z = k * m + v01 = fp_eval_int(f01, z) + v2 = fp_eval_int(f2, z) + xaz = x - (a0 + z) + + term = c_k.mult(BigDecimal(v01).mult(1, prec), prec).div(BigDecimal(v2).mult(1, prec), prec).div(xaz, prec) + sum = sum.add(term, prec) + + # E(z) = prod of (x - a0 - z - j) over the batch, derived from the same + # computed F2 value: E = F2 * (x - a0 - z) / (B * I) with + # B * I = prod of (z + j) * (a0 + z + j) for j = 1 .. m-1. + bik = 1 + (1..m - 1).each {|j| bik *= (z + j) * (a0 + z + j) } + ek = BigDecimal(v2).mult(1, prec).mult(xaz, prec).div(bik, prec) + prod = prod.mult(ek, prec) + + if k < m - 1 + rnum = 1 + rden = 1 + (1..m).each do |j2| + rnum *= n1 - z - j2 + rden *= (z + j2) * (a0 + z + j2) + end + c_k = c_k.mult(rnum * (-b)**m, prec).div(rden, prec) + end + end + sum = sum.mult(BigDecimal(2).power(e01 - e2, prec), prec) unless e01 == e2 + prod = prod.mult(BigDecimal(2).power(m * e2, prec), prec) unless e2.zero? + + # Shift product: batches of (x - i) for i = 0 ... shift, remainder handled directly + if shift > 0 + xi = (x._decimal_shift(fd).to_i << keep) / p10 + leaves = (0...m).map {|j| fp_normalize([xi - j * s2, -s2], -keep, keep) } + while leaves.size > 1 + leaves = leaves.each_slice(2).map {|p, q| q ? fp_mult(p, q, keep) : p } + end + esp = leaves.first + full = shift / m + full.times do |k| + prod = prod.mult(BigDecimal(fp_eval_int(esp, k * m)).mult(1, prec), prec) + end + prod = prod.mult(BigDecimal(2).power(full * esp[1], prec), prec) if full > 0 && !esp[1].zero? + (full * m...shift).each {|i| prod = prod.mult(x - i, prec) } + end + + base = BigDecimal(b).power(x - a0, prec).div(prod.mult(sum, prec), prec) + [base, a0, n1 - 1, 0] + end + + # gamma via the multipoint Lagrange evaluation, for testing. + # Only supports non-integer x >= 0.5 on the Lagrange path; other inputs + # are delegated to the regular implementation. + def self.gamma(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma) + return Gamma.gamma(x, prec) if x < 0.5 || x.frac.zero? + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + base, large_factorial_arg, small_factorial_arg, exp2 = gamma_lagrange(x, prec2) + ans = base.mult(Gamma.integer_factorial(small_factorial_arg, prec2), prec2) + ans = ans.mult(BigDecimal(2).power(exp2, prec2), prec2) unless exp2.zero? + ans.mult(Gamma.integer_factorial(large_factorial_arg, prec2), prec) + end + end + end +end From 97c6c2239f0122aeac7abfab387cc43de33d51db Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 00:51:53 +0900 Subject: [PATCH 03/11] Add remainder-tree multipoint evaluation to the gamma experiment Replace the per-point Horner evaluation (the only PREC^2 term of the pipeline) with classical remainder-tree multipoint evaluation over the falling-factorial subproduct moduli: - power series inverses of the reversed moduli are memoized on the shared subproduct tree, - wide dividends are pre-reduced blockwise with R = t**count mod M_root, whose coefficients stay small, so no division wider than 2*count occurs, - Kronecker slot width now uses the two operands' separate maxima, which also speeds up small-by-large coefficient products elsewhere. Values agree exactly with the Horner path in all tests. Measured crossover against Horner is around m = 700 batches (roughly 250000 digits): below it Horner's machine-word constant wins, so eval_mode defaults to :auto. Co-Authored-By: Claude Fable 5 --- gamma_mp_check.rb | 1 + lib/bigdecimal/math/gamma_multipoint.rb | 157 ++++++++++++++++++++++-- 2 files changed, 147 insertions(+), 11 deletions(-) diff --git a/gamma_mp_check.rb b/gamma_mp_check.rb index b6414b7d..8e6c956d 100644 --- a/gamma_mp_check.rb +++ b/gamma_mp_check.rb @@ -8,6 +8,7 @@ MP = BigMath.const_get(:Gamma)::Multipoint G = BigMath.const_get(:Gamma) +MP.eval_mode = ENV['MP_EVAL'].to_sym if ENV['MP_EVAL'] def rel_err_exp(a, b, prec) e = a.sub(b, prec + 50).div(b, 10).abs diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index ca6df3c7..175f6045 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -27,6 +27,18 @@ module BigMath module Gamma module Multipoint # :nodoc: + # :fast = remainder-tree multipoint evaluation (quasi-linear) + # :horner = per-point Horner (simple; the only PREC^2 term of the pipeline, + # but with a machine-word-size constant) + # :auto = :fast only when the batch count is large enough to win. + # Measured crossover on GMP-backed Ruby is around m = 700 batches, + # i.e. roughly 250000 digits of precision. + FAST_EVAL_MIN_BATCHES = 700 + @eval_mode = :auto + class << self + attr_accessor :eval_mode + end + # ---------- Kronecker substitution convolution on Integer ---------- def self.pack(coeffs, slot_hex) @@ -54,10 +66,11 @@ def self.convolve(a, b) a.each_with_index {|c, i| b.each_with_index {|d, j| out[i + j] += c * d } } return out end - max_bits = 1 - a.each {|c| bits = c.abs.bit_length; max_bits = bits if bits > max_bits } - b.each {|c| bits = c.abs.bit_length; max_bits = bits if bits > max_bits } - w = 2 * max_bits + out_size.bit_length + 2 + max_a = 1 + a.each {|c| bits = c.abs.bit_length; max_a = bits if bits > max_a } + max_b = 1 + b.each {|c| bits = c.abs.bit_length; max_b = bits if bits > max_b } + w = max_a + max_b + out_size.bit_length + 2 slot_hex = (w + 3) / 4 prod = pack_signed(a, slot_hex) * pack_signed(b, slot_hex) unpack_signed(prod, slot_hex, out_size) @@ -116,6 +129,111 @@ def self.fp_eval_int(poly, z) acc end + # ---------- fast evaluation at an arithmetic progression ---------- + # Classical remainder-tree multipoint evaluation. The subproduct moduli for + # consecutive integer points are falling-factorial-type polynomials with + # small exact Integer coefficients (about count * log2(count) bits), which + # keeps the divisions well-scaled. + + def self.fp_neg(p) + [p[0].map {|c| -c }, p[1]] + end + + def self.fp_trunc(p, n) + [p[0][0, n] || [0], p[1]] + end + + def self.fp_mult_trunc(p1, p2, n, keep_bits) + fp_normalize(convolve(p1[0], p2[0])[0, n], p1[1] + p2[1], keep_bits) + end + + # Power series inverse to the given length, by Newton iteration. + # The constant term of f must be exactly 1 (monic reversed modulus). + def self.fp_inv_series(f, terms, keep_bits) + y = [[1], 0] + len = 1 + while len < terms + len = 2 * len < terms ? 2 * len : terms + fy = fp_mult_trunc(fp_trunc(f, len), y, len, keep_bits) + y = fp_mult_trunc(y, fp_add([[2], 0], fp_neg(fy), keep_bits), len, keep_bits) + end + y + end + + # Remainder of fp polynomial r modulo a monic exact-Integer polynomial + # m_int (little-endian coefficient array), via reversal and a precomputed + # power series inverse of the reversed modulus. + def self.fp_rem(r, m_int, inv, keep_bits) + dm = m_int.size - 1 + return r if r[0].size <= dm + ql = r[0].size - dm + qrev = fp_mult_trunc([r[0].reverse, r[1]], fp_trunc(inv, ql), ql, keep_bits) + qm = fp_mult([qrev[0].reverse, qrev[1]], [m_int, 0], keep_bits) + fp_trunc(fp_add(r, fp_neg(qm), keep_bits), dm) + end + + # Tree of exact moduli prod{ t - k } over k = lo ... hi. + # Leaf nodes are [modulus]; internal nodes are [modulus, left, right, nil, nil], + # where the two trailing slots memoize the reversed-modulus inverses of the + # children (shared by all evaluations against the same point set). + def self.subproduct_tree(lo, hi) + return [[-lo, 1]] if hi - lo == 1 + mid = (lo + hi) / 2 + left = subproduct_tree(lo, mid) + right = subproduct_tree(mid, hi) + [convolve(left[0], right[0]), left, right, nil, nil] + end + + def self.eval_descend(r, node, keep_bits, out) + if node.size == 1 + out << [r[0][0] || 0, r[1]] + return + end + left = node[1] + right = node[2] + # A dividend has degree < deg(node modulus), so the inverse length needed + # for division by one child is at most the degree of the other child. + node[3] ||= fp_inv_series([left[0].reverse, 0], right[0].size - 1, keep_bits) + node[4] ||= fp_inv_series([right[0].reverse, 0], left[0].size - 1, keep_bits) + eval_descend(fp_rem(r, left[0], node[3], keep_bits), left, keep_bits, out) + eval_descend(fp_rem(r, right[0], node[4], keep_bits), right, keep_bits, out) + end + + # Values of poly at z = 0, stride, 2*stride, ..., (count-1)*stride. + # Returns [mantissas, exp] with a shared exp. + def self.fp_eval_points(poly, stride, count, keep_bits, tree = nil) + sp = 1 + coeffs = poly[0].map {|c| v = c * sp; sp *= stride; v } + scaled = fp_normalize(coeffs, poly[1], keep_bits) + return [[scaled[0][0] || 0], scaled[1]] if count == 1 + + tree ||= subproduct_tree(0, count) + m_root = tree[0] + + r = scaled + if scaled[0].size > count + # Reduce blockwise: h = h_0 + h_1 * R + h_2 * R**2 + ... (mod M_root) + # with R = t**count mod M_root. R has small coefficients (values of + # t**count at the points are at most count**count), so the only wide + # division is the final reduction of a degree < 2*count polynomial. + root_inv = (tree[5] ||= fp_inv_series([m_root.reverse, 0], count, keep_bits)) + rpow = fp_rem([Array.new(count, 0) + [1], 0], m_root, root_inv, keep_bits) + blocks = scaled[0].each_slice(count).map {|blk| [blk, scaled[1]] } + acc = blocks[0] + rp = rpow + (1...blocks.size).each do |i| + acc = fp_add(acc, fp_mult(blocks[i], rp, keep_bits), keep_bits) + rp = fp_rem(fp_mult(rp, rpow, keep_bits), m_root, root_inv, keep_bits) if i + 1 < blocks.size + end + r = fp_rem(acc, m_root, root_inv, keep_bits) + end + + out = [] + eval_descend(r, tree, keep_bits, out) + emax = out.map {|_, e| e }.max + [out.map {|v, e| e == emax ? v : v >> (emax - e) }, emax] + end + # Guard bits on top of the target precision, absorbing: # - coefficient spread and value dynamic range across batches (~m * log2(n1)) # - rounding of ~log2(m) tree levels and of the evaluation @@ -173,16 +291,25 @@ def self.gamma_lagrange(x, prec) # :nodoc: end f0, f1, f2 = fractions.first f01 = fp_add(f0, f1, keep) - e01 = f01[1] - e2 = f2[1] + fast = eval_mode == :fast || (eval_mode == :auto && m > FAST_EVAL_MIN_BATCHES) + if fast + tree = subproduct_tree(0, m) + v01s, e01 = fp_eval_points(f01, m, m, keep, tree) + v2s, e2 = fp_eval_points(f2, m, m, keep, tree) + else + v01s = Array.new(m) {|k| fp_eval_int(f01, k * m) } + v2s = Array.new(m) {|k| fp_eval_int(f2, k * m) } + e01 = f01[1] + e2 = f2[1] + end sum = BigDecimal(0) prod = BigDecimal(1) c_k = BigDecimal(1) m.times do |k| z = k * m - v01 = fp_eval_int(f01, z) - v2 = fp_eval_int(f2, z) + v01 = v01s[k] + v2 = v2s[k] xaz = x - (a0 + z) term = c_k.mult(BigDecimal(v01).mult(1, prec), prec).div(BigDecimal(v2).mult(1, prec), prec).div(xaz, prec) @@ -218,10 +345,18 @@ def self.gamma_lagrange(x, prec) # :nodoc: end esp = leaves.first full = shift / m - full.times do |k| - prod = prod.mult(BigDecimal(fp_eval_int(esp, k * m)).mult(1, prec), prec) + if full > 0 + if fast + evs, ev = fp_eval_points(esp, m, full, keep) + else + evs = Array.new(full) {|k| fp_eval_int(esp, k * m) } + ev = esp[1] + end + full.times do |k| + prod = prod.mult(BigDecimal(evs[k]).mult(1, prec), prec) + end + prod = prod.mult(BigDecimal(2).power(full * ev, prec), prec) unless ev.zero? end - prod = prod.mult(BigDecimal(2).power(full * esp[1], prec), prec) if full > 0 && !esp[1].zero? (full * m...shift).each {|i| prod = prod.mult(x - i, prec) } end From 4cefc0676261bf9c1d25bb004177f7f59ea9d5b3 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 01:16:39 +0900 Subject: [PATCH 04/11] Reduce Kronecker pack/unpack constants and tighten guard bits Pack folds negative coefficients borrow-style into the next slot, so one hex-join replaces the positive/negative double pack and giant subtraction. Unpack recovers signed slots by borrow propagation instead of adding a giant per-slot bias constant. Guard bits: measured loss with guard = 0 is 2.9 - 3.3 * m * n1.bit_length bits over prec = 300..10000, identical for both eval modes and for near-node x (the value dynamic range across batches dominates all other roundings). Set guard = 4 * m * n1.bit_length + 256, a ~20% margin, and record the measurement in the comment. gamma(sqrt2): 10000 digits 7.4s -> 6.1s, 50000 digits 104.9s -> 91.9s (horner mode); fast mode 123.9s -> 107.6s at 50000 digits. Results still agree exactly with the BSGS implementation. Co-Authored-By: Claude Fable 5 --- lib/bigdecimal/math/gamma_multipoint.rb | 60 ++++++++++++++++++------- 1 file changed, 44 insertions(+), 16 deletions(-) diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index 175f6045..ffe6e996 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -41,21 +41,48 @@ class << self # ---------- Kronecker substitution convolution on Integer ---------- - def self.pack(coeffs, slot_hex) - coeffs.reverse_each.map {|c| c.to_s(16).rjust(slot_hex, '0') }.join.to_i(16) - end - + # Packs signed coefficients as consecutive slot_hex*4-bit slots. + # Negative coefficients are folded borrow-style into the next slot, so a + # single hex-join suffices (no negative-part pack and giant subtraction). def self.pack_signed(coeffs, slot_hex) - v = pack(coeffs.map {|c| c > 0 ? c : 0 }, slot_hex) - v -= pack(coeffs.map {|c| c < 0 ? -c : 0 }, slot_hex) if coeffs.any? {|c| c < 0 } - v + w = slot_hex * 4 + full = 1 << w + borrow = 0 + strs = coeffs.map do |c| + v = c + borrow + if v < 0 + borrow = -1 + v += full + else + borrow = 0 + end + v.to_s(16).rjust(slot_hex, '0') + end + n = strs.reverse!.join.to_i(16) + borrow.zero? ? n : n - (1 << (w * coeffs.size)) end + # Splits n back into signed slot values with borrow propagation + # (slots >= 2**(w-1) are negative), avoiding a giant bias addition. def self.unpack_signed(n, slot_hex, size) - half = 1 << (slot_hex * 4 - 1) - bias = (('8' + '0' * (slot_hex - 1)) * size).to_i(16) - s = (n + bias).to_s(16).rjust(slot_hex * size, '0') - (0...size).map {|i| s[(size - 1 - i) * slot_hex, slot_hex].to_i(16) - half } + w = slot_hex * 4 + half = 1 << (w - 1) + full = 1 << w + neg = n.negative? + s = (neg ? -n : n).to_s(16).rjust(slot_hex * size, '0') + carry = 0 + out = Array.new(size) do |i| + v = s[(size - 1 - i) * slot_hex, slot_hex].to_i(16) + carry + if v >= half + carry = 1 + v - full + else + carry = 0 + v + end + end + out.map! {|v| -v } if neg + out end # Convolution of signed Integer coefficient arrays. @@ -234,12 +261,13 @@ def self.fp_eval_points(poly, stride, count, keep_bits, tree = nil) [out.map {|v, e| e == emax ? v : v >> (emax - e) }, emax] end - # Guard bits on top of the target precision, absorbing: - # - coefficient spread and value dynamic range across batches (~m * log2(n1)) - # - rounding of ~log2(m) tree levels and of the evaluation - # Deliberately generous; to be tightened after error measurements. + # Guard bits on top of the target precision. + # Measured loss with guard = 0 is 2.9 - 3.3 * m * n1.bit_length bits over + # prec = 300 .. 10000, identical for both eval modes and for near-node x: + # the value dynamic range across batches dominates every other rounding. + # 4 * m * n1.bit_length keeps a ~20% multiplicative margin over that. def self.guard_bits(m, n1) - 4 * m * (n1.bit_length + 4) + 256 + 4 * m * n1.bit_length + 256 end # Same contract as Gamma.gamma_lagrange. From 179e875f5672bc2a65a87ddf18d7ee8a1fadc2ea Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 01:46:34 +0900 Subject: [PATCH 05/11] Replace the triple tree with a barycentric pair tree The batch leaf factors decompose as den_j = L_j * B_j * I_j and num_j = -b * L_(j-1) * G_j, where only L_j = (x - a0 - j) - z carries full-precision coefficients and B, I, G are small exact-Integer linears. The series numerator then takes the barycentric form F01 = sum_j (Omega / L_j) * w_j, Omega = prod L_i, with w_j collecting only small-coefficient factors. Tree nodes carry [Omega, Phi, BI, GX] with exact-Integer BI/GX side products, so the only wide-by-wide multiplications per merge are Phi_A * Omega_C and Omega_A * Phi_C against the degree-d Omega, instead of four products of degree-3d triples. The j = 0 term is attached at the root, where its Omega * BI part is F2 itself. gamma(sqrt2): 10000 digits 6.1s -> 4.8s, 50000 digits 91.9s -> 69.9s (2.96x over BSGS; crossover is now around 2000 digits). Results still agree exactly with the BSGS implementation in all tests. Co-Authored-By: Claude Fable 5 --- lib/bigdecimal/math/gamma_multipoint.rb | 73 +++++++++++++++---------- 1 file changed, 43 insertions(+), 30 deletions(-) diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index ffe6e996..fcfca504 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -138,14 +138,9 @@ def self.fp_add(p1, p2, keep_bits) fp_normalize(out, e, keep_bits) end - # Merge of [sum_num, mult_num, den] triples, same as the BSM merge in - # gamma.rb but over polynomials. - def self.triple_merge(a, c, keep_bits) - [ - fp_add(fp_mult(a[0], c[2], keep_bits), fp_mult(a[1], c[0], keep_bits), keep_bits), - fp_mult(a[1], c[1], keep_bits), - fp_mult(a[2], c[2], keep_bits) - ] + # Multiplies an fp polynomial by an exact Integer-coefficient polynomial. + def self.fp_mult_intpoly(p, ip, keep_bits) + fp_normalize(convolve(p[0], ip), p[1], keep_bits) end # Exact Horner evaluation at an integer point. Returns the Integer mantissa; @@ -294,31 +289,49 @@ def self.gamma_lagrange(x, prec) # :nodoc: p10 = 10**fd xa = ((x - a0)._decimal_shift(fd).to_i << keep) / p10 - # Triple tree over leaves t = z + j (j = 1 .. m-1), as polynomials in z: - # den_t = (x - a0 - t) * (t * (a0 + t)) - # num_t = (x - a0 - t + 1) * (-b * (n1 - t)) - identity = [[[1], 0], [[0], 0], [[1], 0]] - fractions = (1..m - 1).map do |j| - xaj = xa - j * s2 - den = fp_normalize( - [xaj * (j * (a0 + j)), xaj * (a0 + 2 * j) - s2 * (j * (a0 + j)), xaj - s2 * (a0 + 2 * j), -s2], - -keep, keep - ) - xaj1 = xaj + s2 - num = fp_normalize( - [-b * xaj1 * (n1 - j), b * (xaj1 + s2 * (n1 - j)), -b * s2], - -keep, keep - ) - [den, num, den] + # Barycentric pair tree over node indices j = 0 .. m-1 (j = 0 carries the + # leading term of the series). The leaf factors decompose as + # den_j = L_j * B_j * I_j, num_j = -b * L_(j-1) * G_j + # where only L_j = (x - a0 - j) - z holds full-precision coefficients; + # B_j = z + j, I_j = a0 + j + z, G_j = n1 - j - z are small. The series + # numerator then becomes + # F01 = sum_j (Omega / L_j) * w_j, Omega = prod L_i, + # with w_j collecting only small-coefficient factors. Each node keeps + # [Omega, Phi, BI, GX] (BI = prod B_i * I_i, GX = (-b)**size * prod G_(i+1), + # both exact Integer polynomials) and merges as + # Omega_P = Omega_A * Omega_C + # Phi_P = (Phi_A * Omega_C) * BI_C + (Omega_A * Phi_C) * GX_A + # so the only wide-by-wide multiplications are with Omega (degree d), + # cheaper than merging [sum, mult, den] triples of degree-3d polynomials. + nodes = (1..m - 1).map do |i| + [ + fp_normalize([xa - i * s2, -s2], -keep, keep), + [[1], 0], + [i * (a0 + i), a0 + 2 * i, 1], + [-b * (n1 - i - 1), b] + ] end - while fractions.size > 1 - fractions = fractions.each_slice(2).map do |p, q| - q ||= identity - triple_merge(p, q, keep) + while nodes.size > 1 + nodes = nodes.each_slice(2).map do |na, nc| + next na unless nc + [ + fp_mult(na[0], nc[0], keep), + fp_add( + fp_mult_intpoly(fp_mult(na[1], nc[0], keep), nc[2], keep), + fp_mult_intpoly(fp_mult(na[0], nc[1], keep), na[3], keep), + keep + ), + convolve(na[2], nc[2]), + convolve(na[3], nc[3]) + ] end end - f0, f1, f2 = fractions.first - f01 = fp_add(f0, f1, keep) + sub = nodes.first + f2 = fp_mult_intpoly(sub[0], sub[2], keep) + # Attach the j = 0 term (Phi = 1, GX = -b * (n1 - 1 - z)): + # its Omega_C * BI_C part is exactly F2. + l0 = fp_normalize([xa, -s2], -keep, keep) + f01 = fp_add(f2, fp_mult_intpoly(fp_mult(l0, sub[1], keep), [-b * (n1 - 1), b], keep), keep) fast = eval_mode == :fast || (eval_mode == :auto && m > FAST_EVAL_MIN_BATCHES) if fast tree = subproduct_tree(0, m) From fdd760a74f0a8ba6361b0d72b9be055009fb0175 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 03:04:35 +0900 Subject: [PATCH 06/11] Wire the multipoint evaluation into the gamma dispatch Gamma.gamma_lagrange now routes full-digit x to the multipoint pipeline when Integer multiplication is GMP-backed (Integer::GMP_VERSION) and prec is at least Multipoint.min_prec (3000, above the measured ~2000-digit crossover against BSGS). Multipoint.enabled = false is the kill switch; without GMP the pipeline stays off automatically because Toom-Cook multiplication would make it asymptotically worse than BSGS. This also removes the trap that reflected arguments (x < 0.5) fell back to the O(PREC^2) BSGS: reflection, lgamma and factorial doubling all reach the dispatch through gamma_lagrange. The test-only Multipoint.gamma wrapper is gone; gamma_mp_check.rb toggles Multipoint.enabled instead. BigMath.gamma(sqrt(2)/3): 8000/16000/32000 digits 4.3/16.8/65.7s -> 2.6/8.3/26.4s. BSM and factorial-doubling paths are unaffected. Co-Authored-By: Claude Fable 5 --- bigdecimal.gemspec | 1 + gamma_mp_check.rb | 30 ++++++++++++++++-------- lib/bigdecimal/math/gamma.rb | 11 ++++++++- lib/bigdecimal/math/gamma_multipoint.rb | 31 +++++++++++++------------ test/bigdecimal/test_bigmath.rb | 4 ++++ 5 files changed, 51 insertions(+), 26 deletions(-) diff --git a/bigdecimal.gemspec b/bigdecimal.gemspec index b5c3c255..5cdb58b3 100644 --- a/bigdecimal.gemspec +++ b/bigdecimal.gemspec @@ -31,6 +31,7 @@ Gem::Specification.new do |s| lib/bigdecimal/math.rb lib/bigdecimal/math/erf.rb lib/bigdecimal/math/gamma.rb + lib/bigdecimal/math/gamma_multipoint.rb lib/bigdecimal/newton.rb lib/bigdecimal/util.rb sample/linear.rb diff --git a/gamma_mp_check.rb b/gamma_mp_check.rb index 8e6c956d..720e6b78 100644 --- a/gamma_mp_check.rb +++ b/gamma_mp_check.rb @@ -1,14 +1,23 @@ # Check & benchmark for the experimental multipoint gamma (gamma_multipoint.rb) # Usage: ruby -Ilib -Itmp/arm64-darwin24/stage/lib gamma_mp_check.rb [mode] # mode: acc (default) | bench | debug +# MP_EVAL=fast|horner forces the evaluation mode require 'bigdecimal' require 'bigdecimal/math' -require 'bigdecimal/math/gamma_multipoint' +require 'bigdecimal/math/gamma' require 'benchmark' MP = BigMath.const_get(:Gamma)::Multipoint -G = BigMath.const_get(:Gamma) MP.eval_mode = ENV['MP_EVAL'].to_sym if ENV['MP_EVAL'] +abort 'multipoint is disabled (Integer::GMP_VERSION not found)' unless MP.enabled +MP.min_prec = 1 # exercise the multipoint path at every precision + +def bsgs_gamma(x, prec) + MP.enabled = false + BigMath.gamma(x, prec) +ensure + MP.enabled = true +end def rel_err_exp(a, b, prec) e = a.sub(b, prec + 50).div(b, 10).abs @@ -19,11 +28,10 @@ def rel_err_exp(a, b, prec) case mode when 'debug' - # Tiny case: compare mp against the regular implementation step by step prec = 50 x = BigDecimal(2).sqrt(150) - a = MP.gamma(x, prec) - b = BigMath.gamma(x, prec + 20) + a = BigMath.gamma(x, prec) + b = bsgs_gamma(x, prec + 20) puts "mp = #{a.to_s("F")[0, 60]}" puts "ref = #{b.to_s("F")[0, 60]}" puts "rel_err_exp = #{rel_err_exp(a, b, prec)}" @@ -34,10 +42,12 @@ def rel_err_exp(a, b, prec) "1/3" => BigDecimal(1).div(3, 2 * prec + 50), "near-node 7+eps" => BigDecimal(7) + BigDecimal(1).div(3, prec + 50)._decimal_shift(-(prec / 2)), "0.6" => BigDecimal("0.6") + BigDecimal(1).div(7, 2 * prec + 50)._decimal_shift(-3), + "reflect sqrt2/3" => BigDecimal(2).sqrt(2 * prec + 50).div(3, 2 * prec + 50), + "reflect -sqrt2" => -BigDecimal(2).sqrt(2 * prec + 50), } cases.each do |name, x| - t_mp = Benchmark.realtime { @mp = MP.gamma(x, prec) } - ref = BigMath.gamma(x, prec + 50) + t_mp = Benchmark.realtime { @mp = BigMath.gamma(x, prec) } + ref = bsgs_gamma(x, prec + 50) e = rel_err_exp(@mp, ref, prec) ok = e == :exact || e <= -prec puts format("%s prec=%-5d %-16s rel_err_exp=%-6s mp=%.2fs", ok ? "OK " : "FAIL", prec, name, e, t_mp) @@ -46,9 +56,9 @@ def rel_err_exp(a, b, prec) when 'bench' [2000, 5000, 10000].each do |prec| x = BigDecimal(2).sqrt(2 * prec + 50) - t_mp = Benchmark.realtime { @mp = MP.gamma(x, prec) } - t_ref = Benchmark.realtime { @ref = BigMath.gamma(x, prec) } - refhi = BigMath.gamma(x, prec + 50) + t_mp = Benchmark.realtime { @mp = BigMath.gamma(x, prec) } + t_ref = Benchmark.realtime { @ref = bsgs_gamma(x, prec) } + refhi = bsgs_gamma(x, prec + 50) puts format("prec=%-6d mp=%.2fs bsgs=%.2fs (%.1fx) mp_err=%s bsgs_err=%s", prec, t_mp, t_ref, t_ref / t_mp, rel_err_exp(@mp, refhi, prec), rel_err_exp(@ref, refhi, prec)) end diff --git a/lib/bigdecimal/math/gamma.rb b/lib/bigdecimal/math/gamma.rb index f2de143d..39eee553 100644 --- a/lib/bigdecimal/math/gamma.rb +++ b/lib/bigdecimal/math/gamma.rb @@ -9,7 +9,9 @@ module BigMath # Lagrange interpolation of f(x) = b**x / x! at integer nodes x_i = b-l, ..., b+l. # BSM(Binary Splitting Method) version for small digit numbers, O(PREC*log(PREC)^3) # BSGS(Baby-Step Giant-Step) version for full digit numbers, O(PREC^2*log(log(PREC))) - # Both orders of magnitude faster than Spouge's approximation which is O(PREC^2*log(PREC)) + # Multipoint evaluation version (gamma_multipoint.rb) replaces BSGS for large PREC + # when Integer multiplication is GMP-backed, O(PREC^1.5*log(PREC)^2) + # All orders of magnitude faster than Spouge's approximation which is O(PREC^2*log(PREC)) # (Complexities assume quasi-linear multiplication, counting large-by-small products # as (n/m) * M(m) = n * log(m) bit ops. BigDecimal multiplies the small coefficients # by schoolbook instead: an extra log factor asymptotically, but faster at any feasible PREC.) @@ -280,6 +282,11 @@ def self.gamma_lagrange_l(b, prec) # Returns [base, large_factorial_arg, small_factorial_arg, exp2] that can produce gamma(x) as: # gamma(x) = base * 2**exp2 * factorial(large_factorial_arg) * factorial(small_factorial_arg) def self.gamma_lagrange(x, prec) + # Full-digit x above the crossover precision: use the experimental multipoint + # evaluation (O(PREC^1.5 * polylog) polynomial pipeline) instead of BSGS. + # See gamma_multipoint.rb; requires GMP-backed Integer multiplication. + return Multipoint.gamma_lagrange(x, prec) if Multipoint.use?(x, prec) + # Shift x to establish a safe center (b) for the barycentric interpolation. # # We must keep all interpolation nodes strictly positive (b - l > 0). Approaching @@ -493,3 +500,5 @@ def self.sinpix(x, pi, prec) private_constant :Gamma end + +require 'bigdecimal/math/gamma_multipoint' diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index fcfca504..0c20cca2 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -35,8 +35,23 @@ module Multipoint # :nodoc: # i.e. roughly 250000 digits of precision. FAST_EVAL_MIN_BATCHES = 700 @eval_mode = :auto + + # Dispatch control (used by Gamma.gamma_lagrange). The multipoint path + # requires GMP-backed Integer multiplication: with Toom-Cook the pipeline + # is asymptotically worse than BSGS. min_prec is the measured crossover + # against BSGS (~2000 digits) with margin. enabled = false is the kill switch. + @enabled = !!defined?(Integer::GMP_VERSION) + @min_prec = 3000 + class << self - attr_accessor :eval_mode + attr_accessor :eval_mode, :enabled, :min_prec + end + + # Same full-digit criterion as the BSM/BSGS branch, applied to the shifted x. + def self.use?(x, prec) + return false unless @enabled && prec >= @min_prec + shift = x < 2 * prec ? 2 * prec - x.floor : 0 + (x + shift - 1).n_significant_digits * prec.bit_length > prec end # ---------- Kronecker substitution convolution on Integer ---------- @@ -405,20 +420,6 @@ def self.gamma_lagrange(x, prec) # :nodoc: [base, a0, n1 - 1, 0] end - # gamma via the multipoint Lagrange evaluation, for testing. - # Only supports non-integer x >= 0.5 on the Lagrange path; other inputs - # are delegated to the regular implementation. - def self.gamma(x, prec) - prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma) - x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma) - return Gamma.gamma(x, prec) if x < 0.5 || x.frac.zero? - - prec2 = prec + BigDecimal::Internal::EXTRA_PREC - base, large_factorial_arg, small_factorial_arg, exp2 = gamma_lagrange(x, prec2) - ans = base.mult(Gamma.integer_factorial(small_factorial_arg, prec2), prec2) - ans = ans.mult(BigDecimal(2).power(exp2, prec2), prec2) unless exp2.zero? - ans.mult(Gamma.integer_factorial(large_factorial_arg, prec2), prec) - end end end end diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index f9891615..e7d87549 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -589,6 +589,8 @@ def test_gamma assert_converge_in_precision {|n| gamma(BigDecimal(5) + BigDecimal("1e-300"), n) } assert_converge_in_precision {|n| gamma(BigDecimal("-3") - BigDecimal("1e-2000"), n) } assert_equal(BigDecimal(24), gamma(BigDecimal(5) + BigDecimal("1e-2000"), 50)) + # crosses the multipoint dispatch threshold (BSGS at 1500, multipoint at 3000 when GMP is available) + assert_converge_in_precision([1500, 3000]) {|n| gamma(BigDecimal(1).div(3, n * 2), n) } end def test_lgamma @@ -617,6 +619,8 @@ def test_lgamma assert_converge_in_precision {|n| lgamma(BigDecimal("-1234.56789"), n).first } assert_converge_in_precision {|n| lgamma(BigDecimal("1e+18"), n).first } assert_converge_in_precision {|n| lgamma(BigDecimal("1e+400"), n).first } + # crosses the multipoint dispatch threshold (BSGS at 1500, multipoint at 3000 when GMP is available) + assert_converge_in_precision([1500, 3000]) {|n| lgamma(BigDecimal(1).div(3, n * 2), n).first } # gamma close 1 or -1 cases assert_converge_in_precision {|n| lgamma(BigDecimal('-3.143580888349980058694358781820227899566'), n).first } From c0ba31fc36f316bd9249937bbb98743203b01767 Mon Sep 17 00:00:00 2001 From: tompng Date: Wed, 29 Jul 2026 02:58:52 +0900 Subject: [PATCH 07/11] Add a value-domain engine (BGS shift of evaluation values) Represent the 2x2 batch transition product P_s(z) = prod [[den_t, 0], [num_t, num_t]] by the values of its entries at z = u * s instead of coefficients, and double via P_2s(z) = P_s(z) * P_s(z + s): the tables are extended by shift of evaluation values (Bostan-Gaudry-Schost) - one convolution with exact binomial weights, small-integer reciprocal kernel and exact incremental delta - then combined pointwise. No product tree and no separate evaluation step remain, so the total cost is a geometric sum over doublings: O(PREC^1.5 * log PREC), one log less than the coefficient engine. S = 2**kappa is even, which also removes the odd node count constraint of the coefficient engine. Measured loss with guard = 0 is 1.35 - 1.49 * S * n1.bit_length bits (prec 300..10000, near-node identical): the feared extrapolation amplification does not appear beyond the table dynamic range, so the guard is set to 2 * S * (bit_length + 4) + 256, smaller than the coefficient engine needs. gamma(sqrt2) 50000 digits: 49.1s vs 65.2s coefficient engine (4.2x over BSGS); crossover between engines is around 7000 digits, so engine defaults to :auto (:values from 8000 digits). Exact agreement with the coefficient engine and BSGS in all tests. Co-Authored-By: Claude Fable 5 --- gamma_mp_check.rb | 1 + lib/bigdecimal/math/gamma_multipoint.rb | 231 +++++++++++++++++++++--- test/bigdecimal/test_bigmath.rb | 2 + 3 files changed, 213 insertions(+), 21 deletions(-) diff --git a/gamma_mp_check.rb b/gamma_mp_check.rb index 720e6b78..cca6bb9d 100644 --- a/gamma_mp_check.rb +++ b/gamma_mp_check.rb @@ -9,6 +9,7 @@ MP = BigMath.const_get(:Gamma)::Multipoint MP.eval_mode = ENV['MP_EVAL'].to_sym if ENV['MP_EVAL'] +MP.engine = ENV['MP_ENGINE'].to_sym if ENV['MP_ENGINE'] abort 'multipoint is disabled (Integer::GMP_VERSION not found)' unless MP.enabled MP.min_prec = 1 # exercise the multipoint path at every precision diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index 0c20cca2..a5a3a44c 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -43,8 +43,14 @@ module Multipoint # :nodoc: @enabled = !!defined?(Integer::GMP_VERSION) @min_prec = 3000 + # :coeff = coefficient domain (barycentric pair tree + multipoint evaluation) + # :values = value domain (BGS shift of evaluation values; no tree, no eval step) + # :auto = :values above its measured crossover against :coeff (~7000 digits) + ENGINE_VALUES_MIN_PREC = 8000 + @engine = :auto + class << self - attr_accessor :eval_mode, :enabled, :min_prec + attr_accessor :eval_mode, :enabled, :min_prec, :engine end # Same full-digit criterion as the BSM/BSGS branch, applied to the shifted x. @@ -271,6 +277,119 @@ def self.fp_eval_points(poly, stride, count, keep_bits, tree = nil) [out.map {|v, e| e == emax ? v : v >> (emax - e) }, emax] end + # ---------- value-domain engine (BGS shift of evaluation values) ---------- + # Instead of polynomial coefficients, the batch transition product + # P_s(z) = prod_{t=z+1..z+s} [[den_t, 0], [num_t, num_t]] + # is represented by the values of its entries at z = u * s (u = 0..3s). + # Doubling: P_2s(z) = P_s(z) * P_s(z + s) needs P_s at u = 0..12s+3, obtained + # by shifting the value table (one convolution); then one pointwise 2x2 + # product per point. Total cost is a geometric sum over doublings instead + # of the log(m) equal-cost levels of the coefficient product tree, and no + # separate evaluation step is needed. + + # Shared kernel for shifting tables of degree d by integer a (a > d): + # reciprocals 1/(a - d + t) in fixed point, exact delta_k = prod (a + k - j) + # and d!. + def self.shift_kernel(d, a, out_len, keep_bits) + rexp = keep_bits + 64 + recips = Array.new(out_len + d) {|t| (1 << rexp) / (a - d + t) } + dfact = (1..d).reduce(1, :*) + deltas = Array.new(out_len) + delta = (a - d..a).reduce(1, :*) + out_len.times do |k| + deltas[k] = delta + delta = delta / (a + k - d) * (a + k + 1) + end + [recips, rexp, deltas, dfact] + end + + # Values Q(a), ..., Q(a + out_len - 1) of the polynomial of degree + # vals.size - 1 given by its values Q(0), ..., Q(d) + # (shift of evaluation values, Bostan-Gaudry-Schost): + # Q(a + k) = (delta_k / d!) * sum_i Q(i) * (-1)**(d-i) * C(d,i) / (a + k - i) + def self.fp_shift_values(table, kernel, keep_bits) + vals, exp = table + d = vals.size - 1 + recips, rexp, deltas, dfact = kernel + comb = 1 + svals = vals.each_with_index.map do |v, i| + sv = comb * ((d - i).odd? ? -v : v) + comb = comb * (d - i) / (i + 1) + sv + end + conv = convolve(svals, recips) + out = Array.new(deltas.size) {|k| conv[k + d] * deltas[k] / dfact } + fp_normalize(out, exp - rexp, keep_bits) + end + + def self.table_concat(t1, t2) + v1, e1 = t1 + v2, e2 = t2 + if e1 > e2 + [v1 + v2.map {|v| v >> (e1 - e2) }, e1] + elsif e2 > e1 + [v1.map {|v| v >> (e2 - e1) } + v2, e2] + else + [v1 + v2, e1] + end + end + + # Builds the value tables [D, N, M] of P_cap_s at z = u * cap_s (u = 0..3*cap_s). + # cap_s must be a power of two. + def self.batch_value_tables(xa, s2, a0, b, n1, cap_s, keep_bits) + dv = [] + nv = [] + (0..3).each do |u| + t = u + 1 + xat = xa - t * s2 + dv << xat * (t * (a0 + t)) + nv << (xat + s2) * (-b * (n1 - t)) + end + dtab = fp_normalize(dv, -keep_bits, keep_bits) + ntab = fp_normalize(nv, -keep_bits, keep_bits) + mtab = ntab + s = 1 + while s < cap_s + kernel = shift_kernel(3 * s, 3 * s + 1, 9 * s + 3, keep_bits) + dvv, de = table_concat(dtab, fp_shift_values(dtab, kernel, keep_bits)) + nvv, ne = table_concat(ntab, fp_shift_values(ntab, kernel, keep_bits)) + mvv, me = table_concat(mtab, fp_shift_values(mtab, kernel, keep_bits)) + # P_2s(u * 2s) = P_s((2u) * s) * P_s((2u + 1) * s), entrywise: + # D' = Dl * Dr, N' = Nl * Dr + Ml * Nr, M' = Ml * Mr + e1 = ne + de + e2 = me + ne + sh = e1 - e2 + l2 = 6 * s + nd = Array.new(l2 + 1) + nn = Array.new(l2 + 1) + nm = Array.new(l2 + 1) + (0..l2).each do |j| + dr = dvv[2 * j + 1] + nr = nvv[2 * j + 1] + t1v = nvv[2 * j] * dr + t2v = mvv[2 * j] * nr + nd[j] = dvv[2 * j] * dr + nn[j] = sh >= 0 ? t1v + (t2v >> sh) : (t1v >> -sh) + t2v + nm[j] = mvv[2 * j] * mvv[2 * j + 1] + end + dtab = fp_normalize(nd, 2 * de, keep_bits) + ntab = fp_normalize(nn, sh >= 0 ? e1 : e2, keep_bits) + mtab = fp_normalize(nm, 2 * me, keep_bits) + s *= 2 + end + [dtab, ntab, mtab] + end + + # Guard bits for the value-domain engine. + # Measured loss with guard = 0 is 1.35 - 1.49 * s * n1.bit_length bits over + # prec = 300 .. 10000, identical for near-node x. The feared extrapolation + # amplification of the value shifts does not appear beyond the table + # dynamic range (the polynomial itself grows at the same rate outside the + # sampled window). 2 * s * (bit_length + 4) keeps a ~1.9x margin. + def self.guard_bits_values(s, n1) + 2 * s * (n1.bit_length + 4) + 256 + end + # Guard bits on top of the target precision. # Measured loss with guard = 0 is 2.9 - 3.3 * m * n1.bit_length bits over # prec = 300 .. 10000, identical for both eval modes and for near-node x: @@ -285,6 +404,10 @@ def self.guard_bits(m, n1) # (m odd so that the barycentric reconstruction keeps positive sign); # slightly wider than the symmetric b-l .. b+l, which only adds accuracy. def self.gamma_lagrange(x, prec) # :nodoc: + if engine == :values || (engine == :auto && prec >= ENGINE_VALUES_MIN_PREC) + return gamma_lagrange_values(x, prec) + end + shift = x < 2 * prec ? 2 * prec - x.floor : 0 x += shift x = BigDecimal(x) - 1 @@ -392,29 +515,95 @@ def self.gamma_lagrange(x, prec) # :nodoc: sum = sum.mult(BigDecimal(2).power(e01 - e2, prec), prec) unless e01 == e2 prod = prod.mult(BigDecimal(2).power(m * e2, prec), prec) unless e2.zero? - # Shift product: batches of (x - i) for i = 0 ... shift, remainder handled directly - if shift > 0 - xi = (x._decimal_shift(fd).to_i << keep) / p10 - leaves = (0...m).map {|j| fp_normalize([xi - j * s2, -s2], -keep, keep) } - while leaves.size > 1 - leaves = leaves.each_slice(2).map {|p, q| q ? fp_mult(p, q, keep) : p } + prod = prod.mult(shift_prod_factor(x, fd, p10, shift, m, keep, prec), prec) if shift > 0 + + base = BigDecimal(b).power(x - a0, prec).div(prod.mult(sum, prec), prec) + [base, a0, n1 - 1, 0] + end + + # Product of (x - i) for i = 0 ... shift - 1: one product polynomial of + # mbatch leaves evaluated at multiples of mbatch, remainder factors direct. + def self.shift_prod_factor(x, fd, p10, shift, mbatch, keep, prec) + prod = BigDecimal(1) + s2 = 1 << keep + xi = (x._decimal_shift(fd).to_i << keep) / p10 + leaves = (0...mbatch).map {|j| fp_normalize([xi - j * s2, -s2], -keep, keep) } + while leaves.size > 1 + leaves = leaves.each_slice(2).map {|p, q| q ? fp_mult(p, q, keep) : p } + end + esp = leaves.first + full = shift / mbatch + if full > 0 + if eval_mode == :fast || (eval_mode == :auto && mbatch > FAST_EVAL_MIN_BATCHES) + evs, ev = fp_eval_points(esp, mbatch, full, keep) + else + evs = Array.new(full) {|k| fp_eval_int(esp, k * mbatch) } + ev = esp[1] end - esp = leaves.first - full = shift / m - if full > 0 - if fast - evs, ev = fp_eval_points(esp, m, full, keep) - else - evs = Array.new(full) {|k| fp_eval_int(esp, k * m) } - ev = esp[1] - end - full.times do |k| - prod = prod.mult(BigDecimal(evs[k]).mult(1, prec), prec) - end - prod = prod.mult(BigDecimal(2).power(full * ev, prec), prec) unless ev.zero? + full.times do |k| + prod = prod.mult(BigDecimal(evs[k]).mult(1, prec), prec) end - (full * m...shift).each {|i| prod = prod.mult(x - i, prec) } + prod = prod.mult(BigDecimal(2).power(full * ev, prec), prec) unless ev.zero? end + (full * mbatch...shift).each {|i| prod = prod.mult(x - i, prec) } + prod + end + + # Same contract as gamma_lagrange, value-domain engine (engine = :values). + # Nodes are A .. A + n1 - 1 with n1 = S * G + 1, S = 2**kappa =~ sqrt(2l). + # S is even, so the barycentric reconstruction sign (-1)**(n1 - 1) is + # positive without any parity constraint on the node count. + def self.gamma_lagrange_values(x, prec) # :nodoc: + shift = x < 2 * prec ? 2 * prec - x.floor : 0 + x += shift + x = BigDecimal(x) - 1 + b = x.round + l = Gamma.gamma_lagrange_l(b, prec) + + kappa = (0.5 * Math.log2(2 * l)).round + kappa = 1 if kappa < 1 + s_cap = 1 << kappa + g = (2 * l + s_cap - 1) / s_cap + n1 = s_cap * g + 1 + a0 = b - l + + keep = Gamma.drop_cap_bits(prec) + guard_bits_values(s_cap, n1) + s2 = 1 << keep + fd = [x.n_significant_digits - x.exponent, 0].max + p10 = 10**fd + xa = ((x - a0)._decimal_shift(fd).to_i << keep) / p10 + + dtab, ntab, mtab = batch_value_tables(xa, s2, a0, b, n1, s_cap, keep) + dvals, ed = dtab + nvals, en = ntab + mvals, em = mtab + + pw_nd = BigDecimal(2).power(en - ed, prec) + pw_md = BigDecimal(2).power(em - ed, prec) + sum_series = BigDecimal(1) + prod = x - a0 + c_k = BigDecimal(1) + g.times do |k| + dk = BigDecimal(dvals[k]).mult(1, prec) + nk = BigDecimal(nvals[k]).mult(1, prec) + sum_series = sum_series.add(c_k.mult(nk, prec).div(dk, prec).mult(pw_nd, prec), prec) + + # Same-value invariant: the batch factor of prod is derived from the + # same computed D_k used in the sum denominator, so near-node errors + # cancel exactly in prod * sum. BI is the exact small-integer part. + bik = 1 + t0 = k * s_cap + (1..s_cap).each {|j| bik *= (t0 + j) * (a0 + t0 + j) } + prod = prod.mult(dk, prec).div(bik, prec) + + if k < g - 1 + c_k = c_k.mult(BigDecimal(mvals[k]).mult(1, prec), prec).div(dk, prec).mult(pw_md, prec) + end + end + prod = prod.mult(BigDecimal(2).power(g * ed, prec), prec) unless ed.zero? + sum = sum_series.div(x - a0, prec) + + prod = prod.mult(shift_prod_factor(x, fd, p10, shift, s_cap, keep, prec), prec) if shift > 0 base = BigDecimal(b).power(x - a0, prec).div(prod.mult(sum, prec), prec) [base, a0, n1 - 1, 0] diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index e7d87549..cb028265 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -591,6 +591,8 @@ def test_gamma assert_equal(BigDecimal(24), gamma(BigDecimal(5) + BigDecimal("1e-2000"), 50)) # crosses the multipoint dispatch threshold (BSGS at 1500, multipoint at 3000 when GMP is available) assert_converge_in_precision([1500, 3000]) {|n| gamma(BigDecimal(1).div(3, n * 2), n) } + # crosses the multipoint engine threshold (coefficient domain at 4000, value domain at 8000) + assert_converge_in_precision([4000, 8000]) {|n| gamma(BigDecimal(1).div(3, n * 2), n) } end def test_lgamma From cce217507e09f9f82a195b9e5c978f5764f0b2ee Mon Sep 17 00:00:00 2001 From: tompng Date: Fri, 31 Jul 2026 03:42:48 +0900 Subject: [PATCH 08/11] Move the shift product to the value domain The shift product prod (x - i) is a single-entry batch factorial, so the same value-table doubling applies with tables of degree s (a fifth of the main tables' work). The batch size is now chosen inside shift_prod_factor as a power of two, independent of the caller's batch count. This removes the last coefficient-domain component from the value engine: gamma(sqrt2, 50000) 49.2s -> 45.0s, and the whole value pipeline is now uniformly O(PREC^1.5 * log PREC). Co-Authored-By: Claude Fable 5 --- lib/bigdecimal/math/gamma_multipoint.rb | 55 ++++++++++++++++--------- 1 file changed, 36 insertions(+), 19 deletions(-) diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index a5a3a44c..1e073c1b 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -515,35 +515,52 @@ def self.gamma_lagrange(x, prec) # :nodoc: sum = sum.mult(BigDecimal(2).power(e01 - e2, prec), prec) unless e01 == e2 prod = prod.mult(BigDecimal(2).power(m * e2, prec), prec) unless e2.zero? - prod = prod.mult(shift_prod_factor(x, fd, p10, shift, m, keep, prec), prec) if shift > 0 + prod = prod.mult(shift_prod_factor(x, fd, p10, shift, keep, prec), prec) if shift > 0 base = BigDecimal(b).power(x - a0, prec).div(prod.mult(sum, prec), prec) [base, a0, n1 - 1, 0] end - # Product of (x - i) for i = 0 ... shift - 1: one product polynomial of - # mbatch leaves evaluated at multiples of mbatch, remainder factors direct. - def self.shift_prod_factor(x, fd, p10, shift, mbatch, keep, prec) + # Value table of the shift-product batches Q_s(z) = prod_{j=0..s-1} (x - z - j) + # at z = u * s (u = 0..s), by the same doubling as batch_value_tables but + # with a single entry of degree s. + def self.shift_value_table(xi, s2, cap_s, keep_bits) + tab = fp_normalize([xi, xi - s2], -keep_bits, keep_bits) + s = 1 + while s < cap_s + kernel = shift_kernel(s, s + 1, 3 * s + 1, keep_bits) + vv, e = table_concat(tab, fp_shift_values(tab, kernel, keep_bits)) + nv = Array.new(2 * s + 1) {|j| vv[2 * j] * vv[2 * j + 1] } + tab = fp_normalize(nv, 2 * e, keep_bits) + s *= 2 + end + tab + end + + # Product of (x - i) for i = 0 ... shift - 1: full batches of power-of-two + # size (chosen here, independent of the caller's batch size) from the + # value table, remainder factors direct. + def self.shift_prod_factor(x, fd, p10, shift, keep, prec) prod = BigDecimal(1) - s2 = 1 << keep - xi = (x._decimal_shift(fd).to_i << keep) / p10 - leaves = (0...mbatch).map {|j| fp_normalize([xi - j * s2, -s2], -keep, keep) } - while leaves.size > 1 - leaves = leaves.each_slice(2).map {|p, q| q ? fp_mult(p, q, keep) : p } + full = 0 + mbatch = 1 + if shift >= 4 + mbatch = 1 << [(0.5 * Math.log2(shift)).round, 1].max + full = shift / mbatch end - esp = leaves.first - full = shift / mbatch if full > 0 - if eval_mode == :fast || (eval_mode == :auto && mbatch > FAST_EVAL_MIN_BATCHES) - evs, ev = fp_eval_points(esp, mbatch, full, keep) - else - evs = Array.new(full) {|k| fp_eval_int(esp, k * mbatch) } - ev = esp[1] + s2 = 1 << keep + xi = (x._decimal_shift(fd).to_i << keep) / p10 + tab = shift_value_table(xi, s2, mbatch, keep) + if full > mbatch + 1 + kernel = shift_kernel(mbatch, mbatch + 1, full - mbatch - 1, keep) + tab = table_concat(tab, fp_shift_values(tab, kernel, keep)) end + vals, e = tab full.times do |k| - prod = prod.mult(BigDecimal(evs[k]).mult(1, prec), prec) + prod = prod.mult(BigDecimal(vals[k]).mult(1, prec), prec) end - prod = prod.mult(BigDecimal(2).power(full * ev, prec), prec) unless ev.zero? + prod = prod.mult(BigDecimal(2).power(full * e, prec), prec) unless e.zero? end (full * mbatch...shift).each {|i| prod = prod.mult(x - i, prec) } prod @@ -603,7 +620,7 @@ def self.gamma_lagrange_values(x, prec) # :nodoc: prod = prod.mult(BigDecimal(2).power(g * ed, prec), prec) unless ed.zero? sum = sum_series.div(x - a0, prec) - prod = prod.mult(shift_prod_factor(x, fd, p10, shift, s_cap, keep, prec), prec) if shift > 0 + prod = prod.mult(shift_prod_factor(x, fd, p10, shift, keep, prec), prec) if shift > 0 base = BigDecimal(b).power(x - a0, prec).div(prod.mult(sum, prec), prec) [base, a0, n1 - 1, 0] From d2984656d48498105d9ec400e15b05eb2e65b1a9 Mon Sep 17 00:00:00 2001 From: tompng Date: Fri, 31 Jul 2026 03:50:25 +0900 Subject: [PATCH 09/11] Retire the coefficient-domain engine The value-domain engine now covers the whole multipoint range: its crossover against BSGS is ~2500 digits, inside the existing min_prec = 3000 dispatch threshold, so the coefficient engine's niche is gone. Remove the barycentric pair tree, the Horner and remainder-tree evaluation modes with their subproduct/inverse-series machinery, the coefficient-polynomial helpers and the engine/eval_mode switches: the file shrinks from ~650 to 337 lines with a single error model (loss = 1.4 * S * n1.bit_length, guard = 2 * S * (bit_length + 4) + 256). The whole pipeline is O(PREC^1.5 * log PREC). Accepting ~1.3x in the 3000-8000 digit band compared to the retired engine buys one engine, one guard law and one code path. Co-Authored-By: Claude Fable 5 --- gamma_mp_check.rb | 3 - lib/bigdecimal/math/gamma.rb | 2 +- lib/bigdecimal/math/gamma_multipoint.rb | 434 ++++-------------------- test/bigdecimal/test_bigmath.rb | 2 - 4 files changed, 71 insertions(+), 370 deletions(-) diff --git a/gamma_mp_check.rb b/gamma_mp_check.rb index cca6bb9d..93674f6d 100644 --- a/gamma_mp_check.rb +++ b/gamma_mp_check.rb @@ -1,15 +1,12 @@ # Check & benchmark for the experimental multipoint gamma (gamma_multipoint.rb) # Usage: ruby -Ilib -Itmp/arm64-darwin24/stage/lib gamma_mp_check.rb [mode] # mode: acc (default) | bench | debug -# MP_EVAL=fast|horner forces the evaluation mode require 'bigdecimal' require 'bigdecimal/math' require 'bigdecimal/math/gamma' require 'benchmark' MP = BigMath.const_get(:Gamma)::Multipoint -MP.eval_mode = ENV['MP_EVAL'].to_sym if ENV['MP_EVAL'] -MP.engine = ENV['MP_ENGINE'].to_sym if ENV['MP_ENGINE'] abort 'multipoint is disabled (Integer::GMP_VERSION not found)' unless MP.enabled MP.min_prec = 1 # exercise the multipoint path at every precision diff --git a/lib/bigdecimal/math/gamma.rb b/lib/bigdecimal/math/gamma.rb index 39eee553..3695f268 100644 --- a/lib/bigdecimal/math/gamma.rb +++ b/lib/bigdecimal/math/gamma.rb @@ -10,7 +10,7 @@ module BigMath # BSM(Binary Splitting Method) version for small digit numbers, O(PREC*log(PREC)^3) # BSGS(Baby-Step Giant-Step) version for full digit numbers, O(PREC^2*log(log(PREC))) # Multipoint evaluation version (gamma_multipoint.rb) replaces BSGS for large PREC - # when Integer multiplication is GMP-backed, O(PREC^1.5*log(PREC)^2) + # when Integer multiplication is GMP-backed, O(PREC^1.5*log(PREC)) # All orders of magnitude faster than Spouge's approximation which is O(PREC^2*log(PREC)) # (Complexities assume quasi-linear multiplication, counting large-by-small products # as (n/m) * M(m) = n * log(m) bit ops. BigDecimal multiplies the small coefficients diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index 1e073c1b..b4aec51d 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -1,25 +1,28 @@ # frozen_string_literal: true -# Experimental multipoint-evaluation version of the Lagrange interpolation -# used by BigMath.gamma, targeting full-digit x. +# Experimental sub-quadratic evaluation of the Lagrange interpolation used by +# BigMath.gamma, targeting full-digit x. # -# The BSGS version in gamma.rb costs O(PREC^2 * polylog): every node needs a -# scalar multiplication against a full-precision power of x. This file evaluates -# the same barycentric sum with sqrt-size batches instead: -# - The [sum_num, mult_num, den] triple of the BSM branch is lifted to -# polynomials in the batch offset z. One triple tree describes all batches. -# - Polynomial arithmetic runs on fixed-point coefficients via Kronecker -# substitution onto Integer multiplication, so it needs quasi-linear Integer -# multiplication (GMP-backed Ruby). -# - The polynomials are evaluated at the arithmetic progression z = 0, mb, -# 2*mb, ... (currently by per-point Horner with small multipliers; a fast -# Newton-basis transform can replace it later). -# Polynomial work is O(PREC^1.5 * polylog). +# The BSGS branch in gamma.rb costs O(PREC^2 * polylog): every interpolation +# node needs a scalar multiplication against a full-precision power of x. +# This file evaluates the same barycentric sum in O(PREC^1.5 * log PREC) with +# sqrt-size batches in the value domain (Bostan-Gaudry-Schost style): +# - The per-node transition is the 2x2 matrix [[den_t, 0], [num_t, num_t]]. +# Its batch product P_s(z) = prod_{t=z+1..z+s} A(t) is represented by the +# values of its entries at z = u * s instead of polynomial coefficients. +# - Doubling P_2s(z) = P_s(z) * P_s(z + s) extends the value tables by +# "shift of evaluation values" (one convolution with exact binomial +# weights and a small-integer reciprocal kernel), then combines pointwise. +# Costs form a geometric sum over doublings - no product tree, no +# separate evaluation step. +# - Values are fixed-point integers with a shared per-table exponent; +# convolutions run via Kronecker substitution onto Integer multiplication, +# so quasi-linear (GMP-backed) Integer multiplication is required. # -# The batch denominator values E(z) used in prod are derived from the same -# computed F2(z) used in sum (E = F2 * (x-A-z) / (B*I) with exact integer B, I), -# so the near-node cancellation between prod and sum stays exact, like the -# batch_prod reuse in the BSGS branch. +# The batch denominator values used in prod are derived from the same computed +# D_k used in sum (E_k = D_k / (B*I) with exact integer B, I), so the near-node +# cancellation between prod and sum stays exact, like the batch_prod reuse in +# the BSGS branch. require 'bigdecimal/math/gamma' @@ -27,30 +30,15 @@ module BigMath module Gamma module Multipoint # :nodoc: - # :fast = remainder-tree multipoint evaluation (quasi-linear) - # :horner = per-point Horner (simple; the only PREC^2 term of the pipeline, - # but with a machine-word-size constant) - # :auto = :fast only when the batch count is large enough to win. - # Measured crossover on GMP-backed Ruby is around m = 700 batches, - # i.e. roughly 250000 digits of precision. - FAST_EVAL_MIN_BATCHES = 700 - @eval_mode = :auto - # Dispatch control (used by Gamma.gamma_lagrange). The multipoint path # requires GMP-backed Integer multiplication: with Toom-Cook the pipeline # is asymptotically worse than BSGS. min_prec is the measured crossover - # against BSGS (~2000 digits) with margin. enabled = false is the kill switch. + # against BSGS (~2500 digits) with margin. enabled = false is the kill switch. @enabled = !!defined?(Integer::GMP_VERSION) @min_prec = 3000 - # :coeff = coefficient domain (barycentric pair tree + multipoint evaluation) - # :values = value domain (BGS shift of evaluation values; no tree, no eval step) - # :auto = :values above its measured crossover against :coeff (~7000 digits) - ENGINE_VALUES_MIN_PREC = 8000 - @engine = :auto - class << self - attr_accessor :eval_mode, :enabled, :min_prec, :engine + attr_accessor :enabled, :min_prec end # Same full-digit criterion as the BSM/BSGS branch, applied to the shifted x. @@ -124,169 +112,19 @@ def self.convolve(a, b) unpack_signed(prod, slot_hex, out_size) end - # ---------- fixed-point polynomials ---------- - # Represented as [coeffs, exp]: sum of coeffs[d] * 2**exp * z**d. - # A single exp per polynomial (fixed-point): small coefficients keep less - # relative precision, which only affects small contributions to the value. + # ---------- fixed-point value tables ---------- + # Represented as [values, exp]: entry i holds values[i] * 2**exp. + # A single exp per table (fixed-point): small entries keep less relative + # precision, which the guard budget absorbs (dynamic-range dominated). - def self.fp_normalize(coeffs, exp, keep_bits) + def self.fp_normalize(values, exp, keep_bits) max = 0 - coeffs.each {|c| bits = c.abs.bit_length; max = bits if bits > max } + values.each {|c| bits = c.abs.bit_length; max = bits if bits > max } s = max - keep_bits - return [coeffs, exp] if s <= 0 - [coeffs.map {|c| c >> s }, exp + s] - end - - def self.fp_mult(p1, p2, keep_bits) - fp_normalize(convolve(p1[0], p2[0]), p1[1] + p2[1], keep_bits) - end - - def self.fp_add(p1, p2, keep_bits) - c1, e1 = p1 - c2, e2 = p2 - if e1 > e2 - c2 = c2.map {|c| c >> (e1 - e2) } - e = e1 - elsif e2 > e1 - c1 = c1.map {|c| c >> (e2 - e1) } - e = e2 - else - e = e1 - end - out = Array.new(c1.size > c2.size ? c1.size : c2.size, 0) - c1.each_with_index {|c, i| out[i] += c } - c2.each_with_index {|c, i| out[i] += c } - fp_normalize(out, e, keep_bits) - end - - # Multiplies an fp polynomial by an exact Integer-coefficient polynomial. - def self.fp_mult_intpoly(p, ip, keep_bits) - fp_normalize(convolve(p[0], ip), p[1], keep_bits) - end - - # Exact Horner evaluation at an integer point. Returns the Integer mantissa; - # the value is mantissa * 2**poly_exp. - def self.fp_eval_int(poly, z) - acc = 0 - poly[0].reverse_each {|c| acc = acc * z + c } - acc - end - - # ---------- fast evaluation at an arithmetic progression ---------- - # Classical remainder-tree multipoint evaluation. The subproduct moduli for - # consecutive integer points are falling-factorial-type polynomials with - # small exact Integer coefficients (about count * log2(count) bits), which - # keeps the divisions well-scaled. - - def self.fp_neg(p) - [p[0].map {|c| -c }, p[1]] - end - - def self.fp_trunc(p, n) - [p[0][0, n] || [0], p[1]] - end - - def self.fp_mult_trunc(p1, p2, n, keep_bits) - fp_normalize(convolve(p1[0], p2[0])[0, n], p1[1] + p2[1], keep_bits) - end - - # Power series inverse to the given length, by Newton iteration. - # The constant term of f must be exactly 1 (monic reversed modulus). - def self.fp_inv_series(f, terms, keep_bits) - y = [[1], 0] - len = 1 - while len < terms - len = 2 * len < terms ? 2 * len : terms - fy = fp_mult_trunc(fp_trunc(f, len), y, len, keep_bits) - y = fp_mult_trunc(y, fp_add([[2], 0], fp_neg(fy), keep_bits), len, keep_bits) - end - y - end - - # Remainder of fp polynomial r modulo a monic exact-Integer polynomial - # m_int (little-endian coefficient array), via reversal and a precomputed - # power series inverse of the reversed modulus. - def self.fp_rem(r, m_int, inv, keep_bits) - dm = m_int.size - 1 - return r if r[0].size <= dm - ql = r[0].size - dm - qrev = fp_mult_trunc([r[0].reverse, r[1]], fp_trunc(inv, ql), ql, keep_bits) - qm = fp_mult([qrev[0].reverse, qrev[1]], [m_int, 0], keep_bits) - fp_trunc(fp_add(r, fp_neg(qm), keep_bits), dm) - end - - # Tree of exact moduli prod{ t - k } over k = lo ... hi. - # Leaf nodes are [modulus]; internal nodes are [modulus, left, right, nil, nil], - # where the two trailing slots memoize the reversed-modulus inverses of the - # children (shared by all evaluations against the same point set). - def self.subproduct_tree(lo, hi) - return [[-lo, 1]] if hi - lo == 1 - mid = (lo + hi) / 2 - left = subproduct_tree(lo, mid) - right = subproduct_tree(mid, hi) - [convolve(left[0], right[0]), left, right, nil, nil] - end - - def self.eval_descend(r, node, keep_bits, out) - if node.size == 1 - out << [r[0][0] || 0, r[1]] - return - end - left = node[1] - right = node[2] - # A dividend has degree < deg(node modulus), so the inverse length needed - # for division by one child is at most the degree of the other child. - node[3] ||= fp_inv_series([left[0].reverse, 0], right[0].size - 1, keep_bits) - node[4] ||= fp_inv_series([right[0].reverse, 0], left[0].size - 1, keep_bits) - eval_descend(fp_rem(r, left[0], node[3], keep_bits), left, keep_bits, out) - eval_descend(fp_rem(r, right[0], node[4], keep_bits), right, keep_bits, out) - end - - # Values of poly at z = 0, stride, 2*stride, ..., (count-1)*stride. - # Returns [mantissas, exp] with a shared exp. - def self.fp_eval_points(poly, stride, count, keep_bits, tree = nil) - sp = 1 - coeffs = poly[0].map {|c| v = c * sp; sp *= stride; v } - scaled = fp_normalize(coeffs, poly[1], keep_bits) - return [[scaled[0][0] || 0], scaled[1]] if count == 1 - - tree ||= subproduct_tree(0, count) - m_root = tree[0] - - r = scaled - if scaled[0].size > count - # Reduce blockwise: h = h_0 + h_1 * R + h_2 * R**2 + ... (mod M_root) - # with R = t**count mod M_root. R has small coefficients (values of - # t**count at the points are at most count**count), so the only wide - # division is the final reduction of a degree < 2*count polynomial. - root_inv = (tree[5] ||= fp_inv_series([m_root.reverse, 0], count, keep_bits)) - rpow = fp_rem([Array.new(count, 0) + [1], 0], m_root, root_inv, keep_bits) - blocks = scaled[0].each_slice(count).map {|blk| [blk, scaled[1]] } - acc = blocks[0] - rp = rpow - (1...blocks.size).each do |i| - acc = fp_add(acc, fp_mult(blocks[i], rp, keep_bits), keep_bits) - rp = fp_rem(fp_mult(rp, rpow, keep_bits), m_root, root_inv, keep_bits) if i + 1 < blocks.size - end - r = fp_rem(acc, m_root, root_inv, keep_bits) - end - - out = [] - eval_descend(r, tree, keep_bits, out) - emax = out.map {|_, e| e }.max - [out.map {|v, e| e == emax ? v : v >> (emax - e) }, emax] + return [values, exp] if s <= 0 + [values.map {|c| c >> s }, exp + s] end - # ---------- value-domain engine (BGS shift of evaluation values) ---------- - # Instead of polynomial coefficients, the batch transition product - # P_s(z) = prod_{t=z+1..z+s} [[den_t, 0], [num_t, num_t]] - # is represented by the values of its entries at z = u * s (u = 0..3s). - # Doubling: P_2s(z) = P_s(z) * P_s(z + s) needs P_s at u = 0..12s+3, obtained - # by shifting the value table (one convolution); then one pointwise 2x2 - # product per point. Total cost is a geometric sum over doublings instead - # of the log(m) equal-cost levels of the coefficient product tree, and no - # separate evaluation step is needed. - # Shared kernel for shifting tables of degree d by integer a (a > d): # reciprocals 1/(a - d + t) in fixed point, exact delta_k = prod (a + k - j) # and d!. @@ -380,140 +218,70 @@ def self.batch_value_tables(xa, s2, a0, b, n1, cap_s, keep_bits) [dtab, ntab, mtab] end - # Guard bits for the value-domain engine. + # Guard bits on top of the target precision. # Measured loss with guard = 0 is 1.35 - 1.49 * s * n1.bit_length bits over - # prec = 300 .. 10000, identical for near-node x. The feared extrapolation - # amplification of the value shifts does not appear beyond the table - # dynamic range (the polynomial itself grows at the same rate outside the - # sampled window). 2 * s * (bit_length + 4) keeps a ~1.9x margin. - def self.guard_bits_values(s, n1) + # prec = 300 .. 10000, identical for near-node x. The extrapolation of the + # value shifts does not amplify errors beyond the table dynamic range (the + # polynomial itself grows at the same rate outside the sampled window). + # 2 * s * (bit_length + 4) keeps a ~1.9x margin. + def self.guard_bits(s, n1) 2 * s * (n1.bit_length + 4) + 256 end - # Guard bits on top of the target precision. - # Measured loss with guard = 0 is 2.9 - 3.3 * m * n1.bit_length bits over - # prec = 300 .. 10000, identical for both eval modes and for near-node x: - # the value dynamic range across batches dominates every other rounding. - # 4 * m * n1.bit_length keeps a ~20% multiplicative margin over that. - def self.guard_bits(m, n1) - 4 * m * n1.bit_length + 256 - end - # Same contract as Gamma.gamma_lagrange. - # Interpolation nodes are A .. A + n1 - 1 with A = b - l and n1 = m**2 - # (m odd so that the barycentric reconstruction keeps positive sign); - # slightly wider than the symmetric b-l .. b+l, which only adds accuracy. + # Interpolation nodes are A .. A + n1 - 1 with n1 = S * G + 1, S = 2**kappa + # =~ sqrt(2l): slightly wider than the symmetric b-l .. b+l, which only + # adds accuracy. S is even, so the barycentric reconstruction sign + # (-1)**(n1 - 1) is positive. def self.gamma_lagrange(x, prec) # :nodoc: - if engine == :values || (engine == :auto && prec >= ENGINE_VALUES_MIN_PREC) - return gamma_lagrange_values(x, prec) - end - shift = x < 2 * prec ? 2 * prec - x.floor : 0 x += shift x = BigDecimal(x) - 1 b = x.round l = Gamma.gamma_lagrange_l(b, prec) - m = Integer.sqrt(2 * l) + 1 - m += 1 if m.even? - n1 = m * m + kappa = (0.5 * Math.log2(2 * l)).round + kappa = 1 if kappa < 1 + s_cap = 1 << kappa + g = (2 * l + s_cap - 1) / s_cap + n1 = s_cap * g + 1 a0 = b - l - keep = Gamma.drop_cap_bits(prec) + guard_bits(m, n1) + keep = Gamma.drop_cap_bits(prec) + guard_bits(s_cap, n1) s2 = 1 << keep - - # Fixed-point mantissas (keep fractional bits) of x - a0 and x fd = [x.n_significant_digits - x.exponent, 0].max p10 = 10**fd xa = ((x - a0)._decimal_shift(fd).to_i << keep) / p10 - # Barycentric pair tree over node indices j = 0 .. m-1 (j = 0 carries the - # leading term of the series). The leaf factors decompose as - # den_j = L_j * B_j * I_j, num_j = -b * L_(j-1) * G_j - # where only L_j = (x - a0 - j) - z holds full-precision coefficients; - # B_j = z + j, I_j = a0 + j + z, G_j = n1 - j - z are small. The series - # numerator then becomes - # F01 = sum_j (Omega / L_j) * w_j, Omega = prod L_i, - # with w_j collecting only small-coefficient factors. Each node keeps - # [Omega, Phi, BI, GX] (BI = prod B_i * I_i, GX = (-b)**size * prod G_(i+1), - # both exact Integer polynomials) and merges as - # Omega_P = Omega_A * Omega_C - # Phi_P = (Phi_A * Omega_C) * BI_C + (Omega_A * Phi_C) * GX_A - # so the only wide-by-wide multiplications are with Omega (degree d), - # cheaper than merging [sum, mult, den] triples of degree-3d polynomials. - nodes = (1..m - 1).map do |i| - [ - fp_normalize([xa - i * s2, -s2], -keep, keep), - [[1], 0], - [i * (a0 + i), a0 + 2 * i, 1], - [-b * (n1 - i - 1), b] - ] - end - while nodes.size > 1 - nodes = nodes.each_slice(2).map do |na, nc| - next na unless nc - [ - fp_mult(na[0], nc[0], keep), - fp_add( - fp_mult_intpoly(fp_mult(na[1], nc[0], keep), nc[2], keep), - fp_mult_intpoly(fp_mult(na[0], nc[1], keep), na[3], keep), - keep - ), - convolve(na[2], nc[2]), - convolve(na[3], nc[3]) - ] - end - end - sub = nodes.first - f2 = fp_mult_intpoly(sub[0], sub[2], keep) - # Attach the j = 0 term (Phi = 1, GX = -b * (n1 - 1 - z)): - # its Omega_C * BI_C part is exactly F2. - l0 = fp_normalize([xa, -s2], -keep, keep) - f01 = fp_add(f2, fp_mult_intpoly(fp_mult(l0, sub[1], keep), [-b * (n1 - 1), b], keep), keep) - fast = eval_mode == :fast || (eval_mode == :auto && m > FAST_EVAL_MIN_BATCHES) - if fast - tree = subproduct_tree(0, m) - v01s, e01 = fp_eval_points(f01, m, m, keep, tree) - v2s, e2 = fp_eval_points(f2, m, m, keep, tree) - else - v01s = Array.new(m) {|k| fp_eval_int(f01, k * m) } - v2s = Array.new(m) {|k| fp_eval_int(f2, k * m) } - e01 = f01[1] - e2 = f2[1] - end + dtab, ntab, mtab = batch_value_tables(xa, s2, a0, b, n1, s_cap, keep) + dvals, ed = dtab + nvals, en = ntab + mvals, em = mtab - sum = BigDecimal(0) - prod = BigDecimal(1) + pw_nd = BigDecimal(2).power(en - ed, prec) + pw_md = BigDecimal(2).power(em - ed, prec) + sum_series = BigDecimal(1) + prod = x - a0 c_k = BigDecimal(1) - m.times do |k| - z = k * m - v01 = v01s[k] - v2 = v2s[k] - xaz = x - (a0 + z) - - term = c_k.mult(BigDecimal(v01).mult(1, prec), prec).div(BigDecimal(v2).mult(1, prec), prec).div(xaz, prec) - sum = sum.add(term, prec) + g.times do |k| + dk = BigDecimal(dvals[k]).mult(1, prec) + nk = BigDecimal(nvals[k]).mult(1, prec) + sum_series = sum_series.add(c_k.mult(nk, prec).div(dk, prec).mult(pw_nd, prec), prec) - # E(z) = prod of (x - a0 - z - j) over the batch, derived from the same - # computed F2 value: E = F2 * (x - a0 - z) / (B * I) with - # B * I = prod of (z + j) * (a0 + z + j) for j = 1 .. m-1. + # Same-value invariant: the batch factor of prod is derived from the + # same computed D_k used in the sum denominator, so near-node errors + # cancel exactly in prod * sum. BI is the exact small-integer part. bik = 1 - (1..m - 1).each {|j| bik *= (z + j) * (a0 + z + j) } - ek = BigDecimal(v2).mult(1, prec).mult(xaz, prec).div(bik, prec) - prod = prod.mult(ek, prec) + t0 = k * s_cap + (1..s_cap).each {|j| bik *= (t0 + j) * (a0 + t0 + j) } + prod = prod.mult(dk, prec).div(bik, prec) - if k < m - 1 - rnum = 1 - rden = 1 - (1..m).each do |j2| - rnum *= n1 - z - j2 - rden *= (z + j2) * (a0 + z + j2) - end - c_k = c_k.mult(rnum * (-b)**m, prec).div(rden, prec) + if k < g - 1 + c_k = c_k.mult(BigDecimal(mvals[k]).mult(1, prec), prec).div(dk, prec).mult(pw_md, prec) end end - sum = sum.mult(BigDecimal(2).power(e01 - e2, prec), prec) unless e01 == e2 - prod = prod.mult(BigDecimal(2).power(m * e2, prec), prec) unless e2.zero? + prod = prod.mult(BigDecimal(2).power(g * ed, prec), prec) unless ed.zero? + sum = sum_series.div(x - a0, prec) prod = prod.mult(shift_prod_factor(x, fd, p10, shift, keep, prec), prec) if shift > 0 @@ -538,8 +306,7 @@ def self.shift_value_table(xi, s2, cap_s, keep_bits) end # Product of (x - i) for i = 0 ... shift - 1: full batches of power-of-two - # size (chosen here, independent of the caller's batch size) from the - # value table, remainder factors direct. + # size from the value table, remainder factors direct. def self.shift_prod_factor(x, fd, p10, shift, keep, prec) prod = BigDecimal(1) full = 0 @@ -565,67 +332,6 @@ def self.shift_prod_factor(x, fd, p10, shift, keep, prec) (full * mbatch...shift).each {|i| prod = prod.mult(x - i, prec) } prod end - - # Same contract as gamma_lagrange, value-domain engine (engine = :values). - # Nodes are A .. A + n1 - 1 with n1 = S * G + 1, S = 2**kappa =~ sqrt(2l). - # S is even, so the barycentric reconstruction sign (-1)**(n1 - 1) is - # positive without any parity constraint on the node count. - def self.gamma_lagrange_values(x, prec) # :nodoc: - shift = x < 2 * prec ? 2 * prec - x.floor : 0 - x += shift - x = BigDecimal(x) - 1 - b = x.round - l = Gamma.gamma_lagrange_l(b, prec) - - kappa = (0.5 * Math.log2(2 * l)).round - kappa = 1 if kappa < 1 - s_cap = 1 << kappa - g = (2 * l + s_cap - 1) / s_cap - n1 = s_cap * g + 1 - a0 = b - l - - keep = Gamma.drop_cap_bits(prec) + guard_bits_values(s_cap, n1) - s2 = 1 << keep - fd = [x.n_significant_digits - x.exponent, 0].max - p10 = 10**fd - xa = ((x - a0)._decimal_shift(fd).to_i << keep) / p10 - - dtab, ntab, mtab = batch_value_tables(xa, s2, a0, b, n1, s_cap, keep) - dvals, ed = dtab - nvals, en = ntab - mvals, em = mtab - - pw_nd = BigDecimal(2).power(en - ed, prec) - pw_md = BigDecimal(2).power(em - ed, prec) - sum_series = BigDecimal(1) - prod = x - a0 - c_k = BigDecimal(1) - g.times do |k| - dk = BigDecimal(dvals[k]).mult(1, prec) - nk = BigDecimal(nvals[k]).mult(1, prec) - sum_series = sum_series.add(c_k.mult(nk, prec).div(dk, prec).mult(pw_nd, prec), prec) - - # Same-value invariant: the batch factor of prod is derived from the - # same computed D_k used in the sum denominator, so near-node errors - # cancel exactly in prod * sum. BI is the exact small-integer part. - bik = 1 - t0 = k * s_cap - (1..s_cap).each {|j| bik *= (t0 + j) * (a0 + t0 + j) } - prod = prod.mult(dk, prec).div(bik, prec) - - if k < g - 1 - c_k = c_k.mult(BigDecimal(mvals[k]).mult(1, prec), prec).div(dk, prec).mult(pw_md, prec) - end - end - prod = prod.mult(BigDecimal(2).power(g * ed, prec), prec) unless ed.zero? - sum = sum_series.div(x - a0, prec) - - prod = prod.mult(shift_prod_factor(x, fd, p10, shift, keep, prec), prec) if shift > 0 - - base = BigDecimal(b).power(x - a0, prec).div(prod.mult(sum, prec), prec) - [base, a0, n1 - 1, 0] - end - end end end diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index cb028265..e7d87549 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -591,8 +591,6 @@ def test_gamma assert_equal(BigDecimal(24), gamma(BigDecimal(5) + BigDecimal("1e-2000"), 50)) # crosses the multipoint dispatch threshold (BSGS at 1500, multipoint at 3000 when GMP is available) assert_converge_in_precision([1500, 3000]) {|n| gamma(BigDecimal(1).div(3, n * 2), n) } - # crosses the multipoint engine threshold (coefficient domain at 4000, value domain at 8000) - assert_converge_in_precision([4000, 8000]) {|n| gamma(BigDecimal(1).div(3, n * 2), n) } end def test_lgamma From 4e6c14699f851652720c6c481e6ac8e0cb215067 Mon Sep 17 00:00:00 2001 From: tompng Date: Fri, 31 Jul 2026 04:47:27 +0900 Subject: [PATCH 10/11] Add the incomplete gamma series as a second accelerator client Generalize batch_value_tables to take the leaf values [den_t, num_t] from a block, making the doubling driver client-independent; the gamma leaf definition moves into gamma_lagrange. incgamma_mp_check.rb computes gamma(x) for full-digit x in [0.5, 3] via gamma(a) =~ r**a * e**-r * (1/a) * (1 + sum prod r/(a+i)), reusing the layer's primitives. The constant numerator degenerates the 2x2 matrix: M_s = r**s is an exact scalar and only two degree-s tables remain, so the doubling is ~3x lighter per term than the gamma client's. Measured: exact agreement with BigMath.gamma at 200..50000 digits, and 1.4x - 2.3x faster than the Lagrange multipoint gamma (0.77s vs 1.74s at 5000 digits, 32.5s vs 44.7s at 50000). Loss law 0.26 - 0.50 * S * bl, smaller than the gamma client's (positive terms, narrower tables). The loss measurement also caught a wrong term-count estimate: the Gaussian tail approximation (1+sqrt(2))*r undercounts by ~13% (the true Poisson tail exponent gamma-(1+gamma)ln(1+gamma) gives ~2.72*r), which had cost a fixed ~28% fraction of the precision. Co-Authored-By: Claude Fable 5 --- incgamma_mp_check.rb | 151 ++++++++++++++++++++++++ lib/bigdecimal/math/gamma_multipoint.rb | 24 ++-- 2 files changed, 167 insertions(+), 8 deletions(-) create mode 100644 incgamma_mp_check.rb diff --git a/incgamma_mp_check.rb b/incgamma_mp_check.rb new file mode 100644 index 00000000..9406cbf3 --- /dev/null +++ b/incgamma_mp_check.rb @@ -0,0 +1,151 @@ +# Second client of the value-domain accelerator layer (gamma_multipoint.rb): +# Gamma via the incomplete gamma series +# gamma(a) =~ gamma_lower(a, r) = r**a * e**-r * (1/a) * S, +# S = 1 + sum_{j>=1} prod_{i=1..j} r / (a + i), r =~ prec * ln(10) +# for full-digit a in [0.5, 3]. The term ratio has a CONSTANT numerator, so the +# 2x2 matrix tables degenerate: M_s = r**s is an exact scalar and only two value +# tables (D, N) of degree s are needed - the doubling is ~3x lighter per term +# than the gamma client's (degree-3 den, three tables). +# +# Usage: ruby -Ilib -Itmp/arm64-darwin24/stage/lib incgamma_mp_check.rb [acc|loss|bench] +require 'bigdecimal' +require 'bigdecimal/math' +require 'bigdecimal/math/gamma' +require 'benchmark' + +BigMath.gamma(BigDecimal('1.5'), 20) +MP = BigMath.const_get(:Gamma)::Multipoint +G = BigMath.const_get(:Gamma) +abort 'requires GMP-backed Integer' unless MP.enabled + +$incg_guard_scale = 1 # measured loss is 0.26-0.50 * S * bl; scale 1 keeps a ~2x margin. loss mode sets 0 + +# Value tables [D, N] of P_s(z) = prod_{t=z+1..z+s} [[a+t, 0], [r, r]] +# at z = u * cap_s (u = 0..cap_s); M_s = r**s is exact and returned separately. +def incg_value_tables(xa, s2, r, cap_s, keep) + dtab = MP.fp_normalize([xa + s2, xa + 2 * s2], -keep, keep) + # Full fixed-point scale even for the exact constant: a tiny-mantissa table + # would force table_concat to rebase the extension down to integer precision. + ntab = MP.fp_normalize([r * s2, r * s2], -keep, keep) + ms = r + s = 1 + while s < cap_s + kernel = MP.shift_kernel(s, s + 1, 3 * s + 1, keep) + dvv, de = MP.table_concat(dtab, MP.fp_shift_values(dtab, kernel, keep)) + nvv, ne = MP.table_concat(ntab, MP.fp_shift_values(ntab, kernel, keep)) + # D' = Dl * Dr, N' = Nl * Dr + M_s * Nr (M_s scalar) + nd = Array.new(2 * s + 1) + nn = Array.new(2 * s + 1) + (0..2 * s).each do |j| + dr = dvv[2 * j + 1] + t1v = nvv[2 * j] * dr + t2v = ms * nvv[2 * j + 1] + nd[j] = dvv[2 * j] * dr + nn[j] = de >= 0 ? t1v + (t2v >> de) : (t1v >> -de) + t2v + end + dtab = MP.fp_normalize(nd, 2 * de, keep) + ntab = MP.fp_normalize(nn, de >= 0 ? ne + de : ne, keep) + ms *= ms + s *= 2 + end + [dtab, ntab, ms] +end + +# gamma(x) for full-digit x in [0.5, 3] via the incomplete gamma series. +def incg_gamma(x, prec) + prec2 = prec + 16 + raise ArgumentError unless x >= 0.5 && x <= 3 + + lr = (prec2 + 20) * Math.log(10) + r = (lr + 2 * Math.log(lr)).ceil + 4 + # Terms until the Poisson-like tail drops below 10**-(prec2+20): + # solve gamma - (1+gamma)*log(1+gamma) = -q. (The Gaussian approximation + # sqrt(2*r*q) underestimates the count by ~13% at q =~ r, costing a fixed + # fraction of the precision.) + q = Math.log(10) * (prec2 + 20) / r + ga = 1.8 + 5.times { ga -= (ga - (1 + ga) * Math.log(1 + ga) + q) / -Math.log(1 + ga) } + nterms = ((1 + ga) * r).ceil + 32 + kappa = [(0.5 * Math.log2(nterms)).round, 1].max + s_cap = 1 << kappa + g = (nterms + s_cap - 1) / s_cap + n_total = s_cap * g + + keep = G.drop_cap_bits(prec2) + $incg_guard_scale * s_cap * (n_total.bit_length + 4) + 256 + s2 = 1 << keep + fd = [x.n_significant_digits - x.exponent, 0].max + xa = (x._decimal_shift(fd).to_i << keep) / 10**fd + + dtab, ntab, ms = incg_value_tables(xa, s2, r, s_cap, keep) + if g > s_cap + 1 + kernel = MP.shift_kernel(s_cap, s_cap + 1, g - s_cap - 1, keep) + dtab = MP.table_concat(dtab, MP.fp_shift_values(dtab, kernel, keep)) + ntab = MP.table_concat(ntab, MP.fp_shift_values(ntab, kernel, keep)) + end + dvals, ed = dtab + nvals, en = ntab + + pw_nd = BigDecimal(2).power(en - ed, prec2) + pw_c = BigDecimal(2).power(-ed, prec2) + sum_series = BigDecimal(1) + c_k = BigDecimal(1) + g.times do |k| + dk = BigDecimal(dvals[k]).mult(1, prec2) + sum_series = sum_series.add(c_k.mult(BigDecimal(nvals[k]).mult(1, prec2), prec2).div(dk, prec2).mult(pw_nd, prec2), prec2) + c_k = c_k.mult(ms, prec2).div(dk, prec2).mult(pw_c, prec2) if k < g - 1 + end + + rpow = BigDecimal(r).power(x, prec2) + emr = BigMath.exp(BigDecimal(-r), prec2) + rpow.mult(emr, prec2).mult(sum_series, prec2).div(x, prec) +end + +def rel_err_exp(a, b, prec) + e = a.sub(b, prec + 50).div(b, 10).abs + e.zero? ? :exact : e.exponent +end + +case ARGV[0] || 'acc' +when 'acc' + [200, 500, 1000, 2000].each do |prec| + { 'sqrt2' => BigDecimal(2).sqrt(2 * prec + 50), + '1+sqrt2/3' => 1 + BigDecimal(2).sqrt(2 * prec + 50).div(3, 2 * prec + 50), + '0.5001-ish' => BigDecimal('0.5') + BigDecimal(1).div(7, 2 * prec + 50)._decimal_shift(-3) }.each do |name, x| + a = incg_gamma(x, prec) + ref = BigMath.gamma(x, prec + 50) + e = rel_err_exp(a, ref, prec) + ok = e == :exact || e <= -(prec - 1) + puts format('%s prec=%-5d %-10s rel_err_exp=%s', ok ? 'OK ' : 'FAIL', prec, name, e) + end + end +when 'loss' + $incg_guard_scale = 0 + [300, 500, 1000, 2000, 5000, 10_000].each do |prec| + x = BigDecimal(2).sqrt(2 * prec + 100) + ref = BigMath.gamma(x, prec + 100) + a = incg_gamma(x, prec) + prec2 = prec + 16 + lr = (prec2 + 20) * Math.log(10) + r = (lr + 2 * Math.log(lr)).ceil + 4 + q = Math.log(10) * (prec2 + 20) / r + ga = 1.8 + 5.times { ga -= (ga - (1 + ga) * Math.log(1 + ga) + q) / -Math.log(1 + ga) } + nterms = ((1 + ga) * r).ceil + 32 + s_cap = 1 << [(0.5 * Math.log2(nterms)).round, 1].max + n_total = s_cap * ((nterms + s_cap - 1) / s_cap) + e = a.sub(ref, prec + 100).div(ref, 10).abs + achieved = e.zero? ? prec + 100 : -e.exponent + loss = ((prec2 + 19 - achieved) * Math.log2(10)).round + puts format('prec=%-6d S=%-4d bl=%-3d loss_bits=%-6d loss/(S*bl)=%.2f', + prec, s_cap, n_total.bit_length, loss, loss.to_f / (s_cap * n_total.bit_length)) + end +when 'bench' + [5000, 10_000, 20_000, 50_000].each do |prec| + x = BigDecimal(2).sqrt(2 * prec + 50) + ti = Benchmark.realtime { @a = incg_gamma(x, prec) } + tg = Benchmark.realtime { @b = BigMath.gamma(x, prec) } + puts format('prec=%-6d incgamma=%.2fs lagrange_mp=%.2fs (%.2fx) agree=%s', + prec, ti, tg, tg / ti, rel_err_exp(@a, @b, prec)) + STDOUT.flush + end +end diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index b4aec51d..d070ea2e 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -172,16 +172,21 @@ def self.table_concat(t1, t2) end end - # Builds the value tables [D, N, M] of P_cap_s at z = u * cap_s (u = 0..3*cap_s). + # Builds the value tables [D, N, M] of the batch transition product + # P_s(z) = prod_{t=z+1..z+s} [[den_t, 0], [num_t, num_t]] + # at z = u * cap_s (u = 0..3*cap_s). This driver is client-independent: + # any series of prefix products of num_t/den_t fits, as long as den_t and + # num_t are polynomials in t. The block yields their fixed-point mantissas + # [den_t, num_t] (at exponent -keep_bits) for t = 1..4; a cubic den and a + # quadratic num are the highest degrees these four samples support. # cap_s must be a power of two. - def self.batch_value_tables(xa, s2, a0, b, n1, cap_s, keep_bits) + def self.batch_value_tables(cap_s, keep_bits) dv = [] nv = [] - (0..3).each do |u| - t = u + 1 - xat = xa - t * s2 - dv << xat * (t * (a0 + t)) - nv << (xat + s2) * (-b * (n1 - t)) + (1..4).each do |t| + d, n = yield(t) + dv << d + nv << n end dtab = fp_normalize(dv, -keep_bits, keep_bits) ntab = fp_normalize(nv, -keep_bits, keep_bits) @@ -253,7 +258,10 @@ def self.gamma_lagrange(x, prec) # :nodoc: p10 = 10**fd xa = ((x - a0)._decimal_shift(fd).to_i << keep) / p10 - dtab, ntab, mtab = batch_value_tables(xa, s2, a0, b, n1, s_cap, keep) + dtab, ntab, mtab = batch_value_tables(s_cap, keep) do |t| + xat = xa - t * s2 + [xat * (t * (a0 + t)), (xat + s2) * (-b * (n1 - t))] + end dvals, ed = dtab nvals, en = ntab mvals, em = mtab From 83685866100485e0e94147ec0b6b801c77522973 Mon Sep 17 00:00:00 2001 From: tompng Date: Sat, 5 Sep 2026 02:32:15 +0900 Subject: [PATCH 11/11] Test the near-node case on the multipoint path The rounding of x before the integer test happens in Gamma.gamma/lgamma ahead of the dispatch, so the multipoint path never sees an integer x. Inside the path the batch denominator D_k is the same computed value in the sum and in prod, so the near-node cancellation is exact. Measured against BSGS at prec 3000..4100 with x within 1e-(prec-10) of a node the error stays below 0.5 ulp. Add a test that crosses the dispatch threshold with such an x. Co-Authored-By: Claude Fable 5.1 --- test/bigdecimal/test_bigmath.rb | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index e7d87549..d54d46cd 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -591,6 +591,7 @@ def test_gamma assert_equal(BigDecimal(24), gamma(BigDecimal(5) + BigDecimal("1e-2000"), 50)) # crosses the multipoint dispatch threshold (BSGS at 1500, multipoint at 3000 when GMP is available) assert_converge_in_precision([1500, 3000]) {|n| gamma(BigDecimal(1).div(3, n * 2), n) } + assert_converge_in_precision([1500, 3000]) {|n| gamma(BigDecimal(5) + BigDecimal("1e-1500"), n) } end def test_lgamma @@ -621,6 +622,7 @@ def test_lgamma assert_converge_in_precision {|n| lgamma(BigDecimal("1e+400"), n).first } # crosses the multipoint dispatch threshold (BSGS at 1500, multipoint at 3000 when GMP is available) assert_converge_in_precision([1500, 3000]) {|n| lgamma(BigDecimal(1).div(3, n * 2), n).first } + assert_converge_in_precision([1500, 3000]) {|n| lgamma(BigDecimal(5) + BigDecimal("1e-1500"), n).first } # gamma close 1 or -1 cases assert_converge_in_precision {|n| lgamma(BigDecimal('-3.143580888349980058694358781820227899566'), n).first }