Experiment: full-digit gamma in O(prec^1.5 log prec) by multipoint evaluation

bigdecimal.gemspec

@@ -30,6 +30,8 @@ Gem::Specification.new do |s|
     lib/bigdecimal/ludcmp.rb
     lib/bigdecimal/math.rb
     lib/bigdecimal/math/erf.rb
+    lib/bigdecimal/math/gamma.rb
+    lib/bigdecimal/math/gamma_multipoint.rb
     lib/bigdecimal/newton.rb
     lib/bigdecimal/util.rb
     sample/linear.rb

gamma_mp_check.rb

@@ -0,0 +1,63 @@
+# Check & benchmark for the experimental multipoint gamma (gamma_multipoint.rb)
+# Usage: ruby -Ilib -Itmp/arm64-darwin24/stage/lib gamma_mp_check.rb [mode]
+#   mode: acc (default) | bench | debug
+require 'bigdecimal'
+require 'bigdecimal/math'
+require 'bigdecimal/math/gamma'
+require 'benchmark'
+
+MP = BigMath.const_get(:Gamma)::Multipoint
+abort 'multipoint is disabled (Integer::GMP_VERSION not found)' unless MP.enabled
+MP.min_prec = 1 # exercise the multipoint path at every precision
+
+def bsgs_gamma(x, prec)
+  MP.enabled = false
+  BigMath.gamma(x, prec)
+ensure
+  MP.enabled = true
+end
+
+def rel_err_exp(a, b, prec)
+  e = a.sub(b, prec + 50).div(b, 10).abs
+  e.zero? ? :exact : e.exponent
+end
+
+mode = ARGV[0] || 'acc'
+
+case mode
+when 'debug'
+  prec = 50
+  x = BigDecimal(2).sqrt(150)
+  a = BigMath.gamma(x, prec)
+  b = bsgs_gamma(x, prec + 20)
+  puts "mp  = #{a.to_s("F")[0, 60]}"
+  puts "ref = #{b.to_s("F")[0, 60]}"
+  puts "rel_err_exp = #{rel_err_exp(a, b, prec)}"
+when 'acc'
+  [100, 200, 500, 1000, 2000].each do |prec|
+    cases = {
+      "sqrt2" => BigDecimal(2).sqrt(2 * prec + 50),
+      "1/3" => BigDecimal(1).div(3, 2 * prec + 50),
+      "near-node 7+eps" => BigDecimal(7) + BigDecimal(1).div(3, prec + 50)._decimal_shift(-(prec / 2)),
+      "0.6" => BigDecimal("0.6") + BigDecimal(1).div(7, 2 * prec + 50)._decimal_shift(-3),
+      "reflect sqrt2/3" => BigDecimal(2).sqrt(2 * prec + 50).div(3, 2 * prec + 50),
+      "reflect -sqrt2" => -BigDecimal(2).sqrt(2 * prec + 50),
+    }
+    cases.each do |name, x|
+      t_mp = Benchmark.realtime { @mp = BigMath.gamma(x, prec) }
+      ref = bsgs_gamma(x, prec + 50)
+      e = rel_err_exp(@mp, ref, prec)
+      ok = e == :exact || e <= -prec
+      puts format("%s prec=%-5d %-16s rel_err_exp=%-6s mp=%.2fs", ok ? "OK  " : "FAIL", prec, name, e, t_mp)
+    end
+  end
+when 'bench'
+  [2000, 5000, 10000].each do |prec|
+    x = BigDecimal(2).sqrt(2 * prec + 50)
+    t_mp = Benchmark.realtime { @mp = BigMath.gamma(x, prec) }
+    t_ref = Benchmark.realtime { @ref = bsgs_gamma(x, prec) }
+    refhi = bsgs_gamma(x, prec + 50)
+    puts format("prec=%-6d mp=%.2fs bsgs=%.2fs (%.1fx) mp_err=%s bsgs_err=%s",
+                prec, t_mp, t_ref, t_ref / t_mp, rel_err_exp(@mp, refhi, prec), rel_err_exp(@ref, refhi, prec))
+  end
+end

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

lib/bigdecimal/math.rb

@@ -626,22 +626,8 @@ def erfc(x, prec)
   #   #=> "0.17724538509055160272981674833411e1"
   #
   def gamma(x, prec)
-    prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma)
-    x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma)
-    prec2 = prec + BigDecimal::Internal::EXTRA_PREC
-    if x < 0.5
-      raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero?
-
-      # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z)
-      pi = PI(prec2)
-      sin = _sinpix(x, pi, prec2)
-      return pi.div(gamma(1 - x, prec2).mult(sin, prec2), prec)
-    elsif x.frac.zero? && x < 1000 * prec
-      return _gamma_positive_integer(x, prec2).mult(1, prec)
-    end
-
-    a, sum = _gamma_spouge_sum_part(x, prec2)
-    (x + (a - 1)).power(x - 0.5, prec2).mult(BigMath.exp(1 - x, prec2), prec2).mult(sum, prec)
+    require 'bigdecimal/math/gamma'
+    Gamma.gamma(x, prec)
   end
 
   # call-seq:
@@ -654,106 +640,8 @@ def gamma(x, prec)
   #   #=> [0.57236494292470008707171367567653e0, 1]
   #
   def lgamma(x, prec)
-    prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma)
-    x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma)
-    prec2 = prec + BigDecimal::Internal::EXTRA_PREC
-    if x < 0.5
-      return [BigDecimal::INFINITY, 1] if x.frac.zero?
-
-      loop do
-        # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z)
-        pi = PI(prec2)
-        sin = _sinpix(x, pi, prec2)
-        log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec)
-        return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC
-
-        # Retry with higher precision if loss of significance is too large
-        prec2 = prec2 * 3 / 2
-      end
-    elsif x.frac.zero? && x < 1000 * prec
-      log_gamma = BigMath.log(_gamma_positive_integer(x, prec2), prec)
-      [log_gamma, 1]
-    else
-      # if x is close to 1 or 2, increase precision to reduce loss of significance
-      diff1_exponent = (x - 1).exponent
-      diff2_exponent = (x - 2).exponent
-      extremely_near_one = diff1_exponent < -prec2
-      extremely_near_two = diff2_exponent < -prec2
-
-      if extremely_near_one || extremely_near_two
-        # If x is extreamely close to base = 1 or 2, linear interpolation is accurate enough.
-        # Taylor expansion at x = base is: (x - base) * digamma(base) + (x - base) ** 2 * trigamma(base) / 2 + ...
-        # And we can ignore (x - base) ** 2 and higher order terms.
-        base = extremely_near_one ? 1 : 2
-        d = BigDecimal(1)._decimal_shift(1 - prec2)
-        log_gamma_d, sign = lgamma(base + d, prec2)
-        return [log_gamma_d.mult(x - base, prec2).div(d, prec), sign]
-      end
-
-      prec2 += [-diff1_exponent, -diff2_exponent, 0].max
-      a, sum = _gamma_spouge_sum_part(x, prec2)
-      log_gamma = BigMath.log(sum, prec2).add((x - 0.5).mult(BigMath.log(x.add(a - 1, prec2), prec2), prec2) + 1 - x, prec)
-      [log_gamma, 1]
-    end
-  end
-
-  # Returns sum part: sqrt(2*pi) and c[k]/(x+k) terms of Spouge's approximation
-  private_class_method def _gamma_spouge_sum_part(x, prec) # :nodoc:
-    x -= 1
-    # Spouge's approximation
-    # x! = (x + a)**(x + 0.5) * exp(-x - a) * (sqrt(2 * pi)  + (1..a - 1).sum{|k| c[k] / (x + k) } + epsilon)
-    # where c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(a - k) / (k - 1)!
-    # and epsilon is bounded by a**(-0.5) * (2 * pi) ** (-a - 0.5)
-
-    # Estimate required a for given precision
-    a = (prec / Math.log10(2 * Math::PI)).ceil
-
-    # Calculate exponent of c[k] in low precision to estimate required precision
-    low_prec = 16
-    log10f = Math.log(10)
-    x_low_prec = x.mult(1, low_prec)
-    loggamma_k = 0
-    ck_exponents = (1..a-1).map do |k|
-      loggamma_k += Math.log10(k - 1) if k > 1
-      -loggamma_k - k / log10f + (k - 0.5) * Math.log10(a - k) - BigDecimal::Internal.float_log(x_low_prec.add(k, low_prec)) / log10f
-    end
-
-    # Estimate exponent of sum by Stirling's approximation
-    approx_sum_exponent = x < 1 ? -Math.log10(a) / 2 : Math.log10(2 * Math::PI) / 2 + x_low_prec.add(0.5, low_prec) * Math.log10(x_low_prec / x_low_prec.add(a, low_prec))
-
-    # Determine required precision of c[k]
-    prec2 = [ck_exponents.max.ceil - approx_sum_exponent.floor, 0].max + prec
-
-    einv = BigMath.exp(-1, prec2)
-    sum = (PI(prec) * 2).sqrt(prec).mult(BigMath.exp(-a, prec), prec)
-    y = BigDecimal(1)
-    (1..a - 1).each do |k|
-      # c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(-k) / (k-1)! / (x + k)
-      y = y.div(1 - k, prec2) if k > 1
-      y = y.mult(einv, prec2)
-      z = y.mult(BigDecimal((a - k) ** k), prec2).div(BigDecimal(a - k).sqrt(prec2).mult(x.add(k, prec2), prec2), prec2)
-      # sum += c[k] / (x + k)
-      sum = sum.add(z, prec2)
-    end
-    [a, sum]
-  end
-
-  private_class_method def _gamma_positive_integer(x, prec) # :nodoc:
-    return x if x == 1
-    numbers = (1..x - 1).map {|i| BigDecimal(i) }
-    while numbers.size > 1
-      numbers = numbers.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a }
-    end
-    numbers.first
-  end
-
-  # Returns sin(pi * x), for gamma reflection formula calculation
-  private_class_method def _sinpix(x, pi, prec) # :nodoc:
-    x = x % 2
-    sign = x > 1 ? -1 : 1
-    x %= 1
-    x = 1 - x if x > 0.5 # to avoid sin(pi*x) loss of precision for x close to 1
-    sign * sin(x.mult(pi, prec), prec)
+    require 'bigdecimal/math/gamma'
+    Gamma.lgamma(x, prec)
   end
 
   # call-seq: