(41) Γ(x+1) [HP-41C]

[Gerson W. Barbosa]

As you may have noticed this Γ(x+1) implementation on the HP-41C is not finished yet. It’s accurate to about nine digits only for x equal or greater than 2. Also, it won’t work for negative arguments. These shortcomings can be easily addressed.
The test HP-75C program below does that by using the little trick in line 25 and by applying Euler’s reflection formula in line 50. Negative integer arguments should return a division by zero error, but because of numerical limitations this won’t occur. This can be fixed, but my primary goal is a compact algorithm using only integer constants, as they don’t take much memory space which are scarce on the HP-41C. I’m not sure that goal has been met as I haven’t included these changes in the RPN program yet. The lack of guard digits might cause some loss of precision on the HP-41C when compared to the HP-75C results.

————-

10 INPUT X
15 P=1 @ Q=0
20 IF X<0 THEN Q=1 @ W=X @ X=-X
25 IF X<3 THEN X=X+3 @ P=X*(X-1)*(X-2)
30 Y=2*X
35 C=EXP(SQR((3*X*(X-5)+55)/330))
40 Z=X+(Y*(7*(C+Y))+6*C)/(Y*(42*(C+Y)-7)+29*C+6)
45 F=SQR(2*Z*PI)*(X/EXP(1))^X/P
50 IF Q THEN F=W*PI/(F*SIN(W*ACOS(-1)))
55 D=DISP F

————-

-71.06 -> -1.08421623(305)E-99
-2.5 -> 2.36327180(084) [4/3×√π]
-2.0 -> 7.55190541728E12 [+∞]
-1.5 -> -3.54490780(301) [-2√π]
-1.0 -> -1.51038108473E13 [-∞]
-0.5-> 1.77245385(344) [√π]
0 -> .999999998537 [1]
.5-> .88622692(4190) [√(π/4)]
1 -> .999999999158 [1]
2 -> 2.00000000000
3 -> 5.99999999122 [6]
4 -> 23.9999999798 [24]
5 -> 120.000000000
6 -> 720.000000(166)
7 -> 5040.00000(104)
8 -> 40320.00000(54)
9 -> 362880.0000(28)
10-> 3628800.000(14)
11-> 39916800.000(6)
12-> 479001600.0(14)
13-> 6227020800.(15)
14-> 8717829120(2.0)
15-> 1.3076743680(3)E12
20-> 2.432902008(30)E18
30-> 2.65252859(798)E32
50-> 3.041409320(56)E64
60-> 8.32098711(335)E81
69-> 1.711224524(19)E98
69.95 -> 9.68284767(215)E99

————-


Edited to fix a few typos in the results table