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.

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

Factorial Doubling for fast calculation of large factorials:

Using Legendre duplication formula, we can calculate factorial(2n) from factorial(n) and factorial(n + 0.5).
Calculating factorial(n + 0.5) is done by 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.24s	BSM
BigMath.gamma(1.25, 100000)	6000s(estimated)	4.2s	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.09s
gamma(1.25)	10000	26.9s	0.69s	0.24s
gamma(1.25)	20000	226s	3.7s	0.57s
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