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

[Gerson W. Barbosa]

Slightly more accurate and smaller (201 bytes on the HP-71B). The HP-71B return a division by zero warning for negative arguments as expected.


————-

10 INPUT X
15 P=1 @ Q=0
20 IF X<0 Q=1 @ W=X @ X=-X
25 IF X<4 THEN X=X+4 @ P=X*(X-1)*(X-2)*(X-3)
30 Y=2*X
35 C=41/30-X/38+X*X/92
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.084216232(57)E-99
-2.5 -> 2.36327180(095) [4/3×√π]
-2.0 -> 9.99999999999E499 [+∞]
-1.5 -> -3.54490770(077) [-2√π]
-1.0 -> -9.99999999999E499 [-∞]
-0.5-> 1.772453850(37) [√π]
0 -> .999999999412 [1]
.5-> .886226925(723) [√(π/4)]
1 -> 1.0000000000(4)
2 -> 2.0000000000(4)
3 -> 6.000000000(42)
4 -> 23.9999999859 [24]
5 -> 120.0000000(48)
6 -> 720.000000(145)
7 -> 5040.000000(35)
8 -> 40320.000000(6)
9 -> 362879.999999 [362880]
10-> 3628799.99998 [3628800]
11-> 39916799.9995 [39916800]
12-> 479001600.00(7)
13-> 6227020800.0(6)
14-> 87178291200.(5)
15-> 1.30767436800E12
20-> 2.4329020081(9)E18
30-> 2.65252859(780)E32
50-> 3.041409320(39)E64
60-> 8.32098711(301)E81
69-> 1.711224524(14)E98
69.95 -> 9.68284767(189)E98

————-


P.S.:

In order to make the program more faithful to the initial HP-41C version, the line 40 should be changed to

40 Z=X+1/(6-1/(Y+6/(7-6/(Y+C))))

This decreases the HP-71B byte count to 190 bytes and has no effect on the overall accuracy, except for very occasional differences of one or two units in the last significant digit:

3 -> 6.000000000(40)
7 -> 5040.000000(34)