Half-precision Γ(x+1) [HP-12C]

[Albert Chan]

Hi, Gerson

With the help of Wolfram Alpha, I have improved on your surd-based CF correction.
Turning series of (Gamma(1+x)/(x/e)^x)^2/(2*pi), x=inf to continued fractions:

\(\large \left({x! \over (x/e)^x}\right)^2 / (2\pi) ≈ x
+ \Large \frac{1}{6-\frac{1}{\left(2x+\frac{77}{90}\right)
+\frac{1}{\left(\frac{16200}{5929}x + \frac{86748660}{246071287}\right) }}} \)

But, we can do much better with corrections to lgamma(x) - ln(stirling(x)):

XCas> series(lgamma(x) - ln(sqrt(2*pi/x)*(x/e)^x), x=inf, 10, polynom)

\(\large {(1/x) \over 12}
- {(1/x)^3 \over 360}
+ {(1/x)^5 \over 1260}
- {(1/x)^7 \over 1680}
+ {(1/x)^9 \over 1188} \)

Even better, coefficients can be easily generated = \(\large{B_{2k} \over (2k)(2k-1)}\)
Example: first coefficient = B₂/(2*1) = 1/12

Above series is usable as-is, but its continued fraction form is better:

\(\frac{1}{12\cdot x
+\large \frac{1}{\frac{5}{2} \cdot x
+\frac{1}{\frac{84}{53} \cdot x
+\frac{1}{\frac{2809}{2340} \cdot x
+\frac{1}{\frac{1003860}{1218947} \cdot x}}}}}\)

Lets compare all 3 continued fraction corrections (c1 = Gerson original version)
Note: "bottom" coefficient adjusted to gives better factorial approximation, for small x.
Since all CF are based from asymptotic series, the adjustment will not affect big x.

For fair comparison, we define x! = Γ(x+1), with same sized x.
(this also gives better factorial approximations, since all assumed big x)

Code:

function c1(x) return 1/(6-1/(2*x+6/(7-5/(3*x+11/20)))) end
function c2(x) return 1/(6-1/(2*x+77/90+1/(x/.366+.46))) end
function gamma(c,x) return sqrt(2*pi*(x+c(x))) * x^(x-1) * exp(-x) end

function c3(x) return 1/(12*x+1/(2.5*x+1/(84/53*x+1/(x/.833+.953/x)))) end
function gamma_c3(x) return sqrt(2*pi) * x^(x-.5) * exp(c3(x)-x) end

function fac_c1(x) return gamma(c1,1+x) end
function fac_c2(x) return gamma(c2,1+x) end
function fac_c3(x) return gamma_c3(1+x) end

lua> fmt = '%g\t%-18.10g%-18.10g%-18.10g' -- show 10 digits precision
lua> for x = 1,10 do print(fmt:format(x, fac_c1(x), fac_c2(x), fac_c3(x))) end

1       1.000015628      0.9999985323      1
2       2.000008196      1.999999694       1.999999998
3       6.000008831      5.999999848       5.999999999
4       24.00001549      23.99999988       24
5       120.0000389      119.9999999       120
6       720.000129       719.9999999       720
7       5040.000538      5040              5040
8       40320.00271      40320             40320
9       362880.0161      362880            362880
10      3628800.11       3628800           3628800


Update: Instead of adjusting c3 coefficients, I kept everything exact, and add a correction.
To keep code short, I turn all constants to integer.

lua> function c3(x) return 1/(12*x+2/(5*x+53/(42*x+1170/(53*x+22999/(429*x+470/x))))) end
lua> for x=1,10 do io.write(('%.10g '):format(fac_c3(x))) end; print()

1 2 6 24 120 720 5040 40320 362880 3628800