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

[tompng]

Status: experimental, not intended to be merged.
This draft PR exists only as a record of an experiment, so that the code, the
measurements and the design notes stay findable from the repository. It is
stacked on the Lagrange interpolation gamma branch (#524); only the
commits after that branch belong to this PR.

Summary

BigMath.gamma / BigMath.lgamma for a full-digit argument (an x with as
many significant digits as the requested precision) currently cost
O(prec² · log log prec) with the BSGS branch of the Lagrange interpolation.
This branch evaluates the same barycentric interpolation sum in
O(prec^1.5 · log prec) and dispatches to it automatically when the
precision is at least 3000 digits and Integer multiplication is GMP-backed.
Results agree with the BSGS branch to the last digit in every test.

Measured on Apple Silicon with GMP 6.3.0, gamma(sqrt(2), prec):

digits	BSGS	this branch	speedup
2,000	0.28s	0.28s	crossover
5,000	2.0s	1.4s	1.5x
10,000	7.9s	4.0s	2.0x
50,000	207s	45s	4.6x

The measured local exponent between 10,000 and 50,000 digits drops from 2.03
(BSGS) to 1.53.

How it works

The Lagrange interpolation of f(x) = b**x / x! at consecutive integer nodes
is a sum whose consecutive term ratio is a rational function of the node
index: rho(t) = -b * (n1 - t) / (t * (a0 + t)) times the ratio of two node
distances. Such a sum is a "matrix factorial": the partial sum state advances
by the 2x2 matrix [[den_t, 0], [num_t, num_t]] per node, and the product of
s consecutive matrices is a polynomial in the batch start z. This is the
setting of Chudnovsky-Chudnovsky / Bostan-Gaudry-Schost (BGS).

The implementation is the value-domain variant of BGS ("shift of evaluation
values"):

• Split the ~2l nodes into S ≈ sqrt(2l) batches of S nodes.
• Represent the batch product P_s(z) by the values of its three entries
  (D, N, M) at z = u * s, not by coefficients.
• Double P_2s(z) = P_s(z) * P_s(z + s): extend the value tables by shifting
  evaluation values (one convolution with exact binomial weights and a
  small-integer reciprocal kernel), then combine pointwise. There is no product
  tree and no separate multipoint evaluation step; the total cost is a
  geometric sum over the doublings.
• Values are fixed-point Integers with one shared exponent per table.
  Convolutions run by Kronecker substitution onto Integer multiplication,
  which is why GMP is required: with Toom-Cook the pipeline is asymptotically
  worse than BSGS, and it stays off automatically without
  Integer::GMP_VERSION.
• The remaining full-precision BigDecimal work is O(S) operations: one
  c_k * N_k / D_k per batch for the sum and one D_k factor per batch for
  the product.
• The shift product prod (x - i) that moves x above 2 * prec uses the
  same doubling with a single-entry table.

Numerical stability

The fixed-point tables lose precision through the dynamic range of the values
across batches, not through the extrapolation of the value shift: the
measured loss with zero guard bits is 1.35 .. 1.49 * S * n1.bit_length
bits, stable over 300 .. 10,000 digits and identical for x next to a node.
The guard is set to 2 * S * (n1.bit_length + 4) + 256 bits.

An x within 10**-q of an interpolation node cancels q digits in one
batch denominator. The batch factor of prod is derived from the same
computed D_k that the sum divides by, so this cancellation is exact in
prod * sum, the same contract as the batch_prod reuse in the BSGS branch.
Measured against BSGS at 3000 .. 4100 digits with x within
10**-(prec-10) of a node, the error stays below 0.5 ulp.

What is in the branch

• lib/bigdecimal/math/gamma_multipoint.rb (345 lines): the engine, plus
  Multipoint.use?, Multipoint.enabled (kill switch) and
  Multipoint.min_prec (3000; the measured crossover against BSGS is
  ~2500 digits).
• lib/bigdecimal/math/gamma.rb: a one-line dispatch at the top of
  Gamma.gamma_lagrange. Reflection (x < 0.5), lgamma and the factorial
  doubling all reach the engine through this single point.
• test/bigdecimal/test_bigmath.rb: convergence tests that cross the dispatch
  threshold, including an x next to a node.
• gamma_mp_check.rb: accuracy sweep and benchmark against BSGS
  (Multipoint.enabled toggle).
• incgamma_mp_check.rb: a second client of the same doubling driver, the
  incomplete gamma series gamma(a) ≈ r**a * e**-r / a * (1 + sum prod r/(a+i))
  for full-digit a in [0.5, 3]. Its constant numerator degenerates the 2x2
  matrix to two tables, and it runs 1.4x .. 2.3x faster than the Lagrange
  engine at 5,000 .. 50,000 digits with exact agreement. Left as a script; it
  suggests that a future lane split (incomplete gamma series for small x,
  Lagrange only for the duplication formula at large x) would be simpler and
  faster than this PR's single engine.

History of the experiment

The commit sequence records the path, including the parts that were retired:

1. Coefficient-domain engine: batch polynomials built by a product tree,
   evaluated per point by Horner (O(prec^1.5 · log² prec)).
2. Remainder-tree multipoint evaluation; measured crossover against Horner was
   ~250,000 digits, so it never won in practice.
3. Kronecker pack/unpack constants (borrow-folded signed packing) and a guard
   formula fitted to measured loss.
4. Barycentric pair tree: only two wide-by-wide products per merge instead of
   four (50,000 digits: 92s to 70s).
5. Dispatch wiring into Gamma.gamma_lagrange.
6. Value-domain engine (this PR's final form), then the shift product moved to
   the value domain, then the coefficient-domain engine removed: one engine,
   one error model, one code path (650 to 345 lines).
7. Incomplete gamma series as a second client of the driver.

Why this is not being merged

• It only helps full-digit arguments above ~2500 digits on GMP builds. The
  BSM and BSGS branches cover the common cases and stay unchanged.
• The engine is a fixed-point Integer pipeline with an empirically fitted
  guard formula. It is verified over 300 .. 50,000 digits, not proved, and the
  guard has not been re-checked for S > 1000 (above ~10^6 digits).
• Memory: the value tables hold ~prec^1.5 digits, hundreds of MB at 10^5
  digits.
• I do not have the capacity to maintain it. The design notes and
  measurements are in the commit messages, so the branch can serve as a
  starting point if someone wants to pick it up.

Try it

ruby -Ilib -e 'require "bigdecimal/math"; require "benchmark";
  x = BigMath.sqrt(BigDecimal(2), 20000)
  puts Benchmark.realtime { BigMath.gamma(x, 20000) }
  BigMath.const_get(:Gamma)::Multipoint.enabled = false
  puts Benchmark.realtime { BigMath.gamma(x, 20000) }'
ruby -Ilib gamma_mp_check.rb acc     # accuracy sweep against BSGS
ruby -Ilib gamma_mp_check.rb bench   # timings