diff --git a/bigdecimal.gemspec b/bigdecimal.gemspec index 38a0f6d6..5cdb58b3 100644 --- a/bigdecimal.gemspec +++ b/bigdecimal.gemspec @@ -30,6 +30,8 @@ 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/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 new file mode 100644 index 00000000..93674f6d --- /dev/null +++ b/gamma_mp_check.rb @@ -0,0 +1,63 @@ +# 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' +require 'benchmark' + +MP = BigMath.const_get(:Gamma)::Multipoint +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 + e.zero? ? :exact : e.exponent +end + +mode = ARGV[0] || 'acc' + +case mode +when 'debug' + prec = 50 + x = BigDecimal(2).sqrt(150) + 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)}" +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), + "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 = 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) + end + end +when 'bench' + [2000, 5000, 10000].each do |prec| + x = BigDecimal(2).sqrt(2 * 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 +end 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.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..3695f268 --- /dev/null +++ b/lib/bigdecimal/math/gamma.rb @@ -0,0 +1,504 @@ +# 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))) + # Multipoint evaluation version (gamma_multipoint.rb) replaces BSGS for large PREC + # 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 + # 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) + # 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 + # 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 + +require 'bigdecimal/math/gamma_multipoint' diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb new file mode 100644 index 00000000..d070ea2e --- /dev/null +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -0,0 +1,345 @@ +# frozen_string_literal: true + +# Experimental sub-quadratic evaluation of the Lagrange interpolation used by +# BigMath.gamma, targeting full-digit x. +# +# 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 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' + +module BigMath + module Gamma + module Multipoint # :nodoc: + + # 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 (~2500 digits) with margin. enabled = false is the kill switch. + @enabled = !!defined?(Integer::GMP_VERSION) + @min_prec = 3000 + + class << self + attr_accessor :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 ---------- + + # 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) + 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) + 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. + 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_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) + end + + # ---------- 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(values, exp, keep_bits) + max = 0 + values.each {|c| bits = c.abs.bit_length; max = bits if bits > max } + s = max - keep_bits + return [values, exp] if s <= 0 + [values.map {|c| c >> s }, exp + s] + end + + # 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 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(cap_s, keep_bits) + dv = [] + nv = [] + (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) + 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 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 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 + + # Same contract as Gamma.gamma_lagrange. + # 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: + 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(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(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 + + 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 + + # 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 from the value table, remainder factors direct. + def self.shift_prod_factor(x, fd, p10, shift, keep, prec) + prod = BigDecimal(1) + full = 0 + mbatch = 1 + if shift >= 4 + mbatch = 1 << [(0.5 * Math.log2(shift)).round, 1].max + full = shift / mbatch + end + if full > 0 + 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(vals[k]).mult(1, prec), prec) + end + 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 + end + end + end +end diff --git a/test/bigdecimal/test_bigmath.rb b/test/bigdecimal/test_bigmath.rb index 05bcbc82..d54d46cd 100644 --- a/test/bigdecimal/test_bigmath.rb +++ b/test/bigdecimal/test_bigmath.rb @@ -575,11 +575,23 @@ 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)) + # 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 @@ -595,23 +607,32 @@ 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 } + # 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(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