Faster gamma calculation

[tompng]

This is a proof of concept implementation of faster gamma calculation using localized Lagrange interpolation of b**x / x!

Approach in this PR

Directly interpolating gamma(x) with equidistant nodes fails due to its singularities and Runge's phenomenon, so interpolate the following scaled reciprocal function instead:
$$f(x) = \frac{b^x}{x!}$$
Scaling by b**x shapes the function into a nearly symmetric bell curve centered at b.
By performing Lagrange interpolation at x = b-l, b-l+1, ..., b+l, we achieve highly accurate localized interpolation. The center b is dynamically determined based on the input x.
The formula is simple, node values f(integer) can be easily calculated, and it opens up room for various optimizations.

Design of f(x)

To calculate gamma(integer + small_rational) with Binary Splitting Method, f(n+1) needs to be calculated easily from f(n), such as making f(n+1)/f(n) a small rational.
Candidates of f(x) are x!/b**x and b**x/x!. x!/b**x requires constant times more interpolation nodes.

Other functions don't have the ability to calculate f(n+1) from f(n) easily:

• x!/x**x, x!/x**b: Simple, but hard to optimize
• x!*exp(x)/x**(x+1/2): Gives a formula similar to Spouge's approximation
• (x+a-1)!*exp(-x-a)/(x+a)**(x+1/2) at nodes -a-1..-1: Gives exactly the same formula as Spouge's approximation when modified to converge at x → ∞

Algorithm overview:

Lagrange interpolation of f(x) = b**x / x!

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 magnitude of order faster than Spouge's approximation which costs O(PREC^2*log(PREC))
Requires fast calculation of factorial(nearly_x_integer).

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 precision. Measured time grows like PREC^2.

Factorial Doubling for fast calculation of large factorials:

Using Legendre duplication formula, we can calculate factorial(2n) from factorial(n) and factorial(n + 0.5).
Calculating factorial(n + 0.5) is done by an optimized BSM version of Lagrange interpolation in quasi-linear time.
This will drastically reduce the cost of calculating large factorials.
O(PREC*log(PREC)^3*log(factorial_argument))

Stirling's approximation with Bernoulli numbers

Only used in lgamma when x is extremely large, such as:

BigMath.lgamma(10000000000000000, 10000);

Benchmark

Calculation	master branch	This PR	note
BigMath.gamma(1.25, 10000)	49s	0.17s	BSM
BigMath.gamma(1.25, 100000)	6000s(estimated)	2.3s	BSM
BigMath.gamma(BigDecimal(1).div(3, 10000), 10000)	56s	8.2s	Full-digit, BSGS
BigMath.gamma(BigDecimal(1).div(3, 100000), 100000)	7000s(estimated)	803s	Full-digit, BSGS
BigMath.gamma(10**17, 10000)	71s	1.1s	Factorial Doubling
BigMath.gamma(10**17, 100000)	9000s(estimated)	13s	Factorial Doubling

Comparison with mpmath(gmpy-backend)

Calculation	digits	mpmath first run	mpmath second run (cached)	This PR
gamma(1.25)	5000	3.2s	0.12s	0.08s
gamma(1.25)	10000	26.9s	0.69s	0.17s
gamma(1.25)	20000	226s	3.7s	0.34s
gamma(1/3)	5000	3.2s	0.12s	2.2s
gamma(1/3)	10000	24s	0.68s	8.2s
gamma(1/3)	20000	226s	3.8s	33s

Correctness

Error of lagrange interpolation on integer nodes will be maximum at the midpoint of the nodes, e.g.: x = n + 0.5
BigMath.gamma 0.5, 10000000 compared with pi.sqrt(10000000 + 100) gives 10000000 + 15 digits correct.
It internally calculates 10000000 + 16 digits, 10000000 + 15 digits correct, final result is rounded with 10000000 digits.

Other cases:
BigMath.gamma 10000000000.5, 1000000: Factorial doubling
BigMath.gamma(BigMath.sqrt(2, 100000), 100000): BSGS
Both gives prec + 15 accurate digits compared with gamma(x, prec + 100)