From 8f1d59df614dff86ec3eab3d2e8afd47c864cdd9 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 00:17:26 +0900 Subject: [PATCH 01/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 c2f8569cd7d701fe68bb20272751d839fd399837 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 00:51:53 +0900 Subject: [PATCH 02/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 fe7df905150bd0654e62311a80bd2a1e0d3f8cc5 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 01:16:39 +0900 Subject: [PATCH 03/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 5b69342756e51dd1b287b3d40330dbb0f187c720 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 01:46:34 +0900 Subject: [PATCH 04/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 d18673b34594af32b53eda7c9c939655d7102441 Mon Sep 17 00:00:00 2001 From: tompng Date: Tue, 28 Jul 2026 03:04:35 +0900 Subject: [PATCH 05/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 18325474..81b46aab 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.) @@ -293,6 +295,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 @@ -511,3 +518,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 2261b5c3471c11714309a5be5a07e21992fc5e98 Mon Sep 17 00:00:00 2001 From: tompng Date: Wed, 29 Jul 2026 02:58:52 +0900 Subject: [PATCH 06/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 380a7a2e8c2716f25f08a65ef83eada9c155a2d0 Mon Sep 17 00:00:00 2001 From: tompng Date: Fri, 31 Jul 2026 03:42:48 +0900 Subject: [PATCH 07/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 7f7af2849c31bdcd06f2e53ea53b0b57c28b38af Mon Sep 17 00:00:00 2001 From: tompng Date: Fri, 31 Jul 2026 03:50:25 +0900 Subject: [PATCH 08/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 81b46aab..ae25d742 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 294446e526a9d25e628b4c12a94fe2950b86f88a Mon Sep 17 00:00:00 2001 From: tompng Date: Fri, 31 Jul 2026 04:47:27 +0900 Subject: [PATCH 09/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 5813465df0a0d22005502ecaa03cfd3055e13058 Mon Sep 17 00:00:00 2001 From: tompng Date: Sat, 5 Sep 2026 02:32:15 +0900 Subject: [PATCH 10/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 } From acfb5a8623be5ef93d67c773084caf8d933f21e3 Mon Sep 17 00:00:00 2001 From: tompng Date: Sun, 13 Sep 2026 13:28:39 +0900 Subject: [PATCH 11/11] Accumulate per-batch products with extra digits in the multipoint path Same fix as the BSGS branch (60829da0): sum_series, c_k and prod are rounded once per batch, so their errors grow with the batch count g. Accumulate them, and the table values they take in, with log10(g) extra digits; shift_prod_factor likewise with its own factor count. Measured at the gamma_lagrange level against a 40-digit-higher BSGS reference, x = sqrt(2): 13.4 / 19.0 / 40.7 ulp -> 2.9 / 1.0 / 2.7 ulp at 3000 / 5000 / 10000 digits, the same level as the BSGS branch (4.2 / 1.6 / 3.0 ulp). Co-Authored-By: Claude Fable 5.1 --- lib/bigdecimal/math/gamma_multipoint.rb | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/lib/bigdecimal/math/gamma_multipoint.rb b/lib/bigdecimal/math/gamma_multipoint.rb index d070ea2e..fd85aeac 100644 --- a/lib/bigdecimal/math/gamma_multipoint.rb +++ b/lib/bigdecimal/math/gamma_multipoint.rb @@ -271,10 +271,13 @@ def self.gamma_lagrange(x, prec) # :nodoc: sum_series = BigDecimal(1) prod = x - a0 c_k = BigDecimal(1) + # sum_series, c_k and prod are updated once per batch, so their rounding + # errors accumulate in proportion to the number of batches. + accumulate_prec = prec + Math.log10(g).ceil + 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) + dk = BigDecimal(dvals[k]).mult(1, accumulate_prec) + nk = BigDecimal(nvals[k]).mult(1, accumulate_prec) + sum_series = sum_series.add(c_k.mult(nk, accumulate_prec).div(dk, accumulate_prec).mult(pw_nd, accumulate_prec), accumulate_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 @@ -282,10 +285,10 @@ def self.gamma_lagrange(x, prec) # :nodoc: 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) + prod = prod.mult(dk, accumulate_prec).div(bik, accumulate_prec) if k < g - 1 - c_k = c_k.mult(BigDecimal(mvals[k]).mult(1, prec), prec).div(dk, prec).mult(pw_md, prec) + c_k = c_k.mult(BigDecimal(mvals[k]).mult(1, accumulate_prec), accumulate_prec).div(dk, accumulate_prec).mult(pw_md, accumulate_prec) end end prod = prod.mult(BigDecimal(2).power(g * ed, prec), prec) unless ed.zero? @@ -323,6 +326,9 @@ def self.shift_prod_factor(x, fd, p10, shift, keep, prec) mbatch = 1 << [(0.5 * Math.log2(shift)).round, 1].max full = shift / mbatch end + # prod takes full batch factors and up to mbatch - 1 remainder factors, + # each rounded once. + accumulate_prec = prec + Math.log10(full + mbatch).ceil + 1 if full > 0 s2 = 1 << keep xi = (x._decimal_shift(fd).to_i << keep) / p10 @@ -333,12 +339,12 @@ def self.shift_prod_factor(x, fd, p10, shift, keep, prec) end vals, e = tab full.times do |k| - prod = prod.mult(BigDecimal(vals[k]).mult(1, prec), prec) + prod = prod.mult(BigDecimal(vals[k]).mult(1, accumulate_prec), accumulate_prec) end - prod = prod.mult(BigDecimal(2).power(full * e, prec), prec) unless e.zero? + prod = prod.mult(BigDecimal(2).power(full * e, prec), accumulate_prec) unless e.zero? end - (full * mbatch...shift).each {|i| prod = prod.mult(x - i, prec) } - prod + (full * mbatch...shift).each {|i| prod = prod.mult(x - i, accumulate_prec) } + prod.mult(1, prec) end end end