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

[Albert Chan]

Another approach is not to use asymptotic formula at all, and use Lanczos algorithm.
see The Log Gamma Function with C#, for comparison of different methods.

I don't understand how it work, but Lanczos is amazing !


Translated code from above link, for HP-71B:

Code:

10 DIM C(6) @ K=LN(2*PI)/2-5
20 C(1)=76.1800917295  @ C(2)=-86.5053203294  @ C(3)=24.0140982408
30 C(4)=-1.23173957245 @ C(5)=.00120865097387 @ C(6)=-5.3952393849E-6
40 DEF FNS(X) @ S=1.00000000019
50 FOR I=1 TO 6 @ S=S+C(I)/(X+I) @ NEXT I
60 FNS=S/X @ END DEF
70 DEF FNL(X)=K+LN(FNS(X))+(X+.5)*(LN(X+5.5)-1)    ! = lgamma
80 DEF FNG(X)=EXP(K-X-.5)*FNS(X)*(X+5.5)^(X+.5)    ! = tgamma

>RUN
>FOR X=10 TO 90 STEP 10 @ X, FNL(X), LN(GAMMA(X))-RES @ NEXT X
10       12.8018274802      -.0000000001
20       39.3398841871       .0000000001
30       71.2570389671       .0000000001
40       106.63176026        .000000001
50       144.565743946       0
60       184.533828862      -.000000001
70       226.190548324       0
80       269.291097651       0
90       313.65282995        0

Below is doing *negative* fractional factorials !

>FOR X=.1 TO .9 STEP .1 @ X, FNG(X), GAMMA(X)-RES @ NEXT X
.1       9.51350769862       .00000000005
.2       4.59084371199       .00000000001
.3       2.99156898768       .00000000001
.4       2.21815954376       0
.5       1.7724538509        .00000000001
.6       1.48919224879       .00000000002
.7       1.29805533263       .00000000002
.8       1.16422971372       .00000000001
.9       1.06862870212       0