Threaded Mode | Linear Mode

Gerson W. Barbosa · (This post was last modified: 03-20-2021 08:40 PM by Gene.)

01 LBL "GXP1"
02 ENTER
03 ENTER
04 6
05 1/X
06 RCL Y
07 22
08 /
09 -
10 RCL Y
11 X^2
12 110
13 /
14 +
15 SQRT
16 E^X
17 X<>Y
18 ST+ Z
19 ST+ T
20 RDN
21 +
22 6
23 X<>Y
24 /
25 CHS
26 7
27 +
28 6
29 X<>Y
30 /
31 +
32 1/X
33 CHS
34 6
35 +
36 1/X
37 +
38 ST+ X
39 X<>Y
40 R^
41 -1
42 E^X
43 *
44 X<>Y
45 Y^X
46 X<>Y
47 PI
48 *
49 SQRT
50 *
51 END

Γ(x+1), x = 1..69

RUN
0,999998868
1,999999993
6,000000000
24,00000001
120,0000001
719,9999998
5.039,999998
40.320,00008
362.880,0005
3.628.800,004
39.916.800,02
479.001.600,0
6.227.020.813,
8,717829135+10
1,307674370+12
2,092278990+13
3,556874281+14
6,402373724+15
1,216451007+17
2,432902011+18
5,109094224+19
1,124000729+21
2,585201681+22
6,204484035+23
1,551121007+25
4,032914619+26
1,088886946+28
3,048883422+29
8,841761897+30
2,652528633+32
8,222838734+33
2,631308388+35
8,683317650+36
2,952327992+38
1,033314794+40
3,719933246+41
1,376375297+43
5,230226253+44
2,039788233+46
8,159152904+47
3,345252679+49
1,405006121+51
6,041526302+52
2,658271563+54
1,196222200+56
5,502622253+57
2,586232450+59
1,241391572+61
6,082818683+62
3,041409334+64
1,551118754+66
8,065817499+67
4,274883258+69
2,308436953+71
1,269640355+73
7,109985964+74
4,052691987+76
2,350561344+78
1,386831189+80
8,320987107+81
5,075802119+83
3,146997305+85
1,982608348+87
1,268869339+89
8,247650677+90
5,443449430+92
3,647111106+94
2,480035546+96
1,711224521+98

Γ(x+1)/x!, x = 1..69

RUN
9,999988683-01
9,999999965-01
1,000000000+00
1,000000000+00
1,000000001+00
9,999999997-01
9,999999996-01
1,000000002+00
1,000000001+00
1,000000001+00
1,000000001+00
1,000000000+00
1,000000002+00
1,000000002+00
1,000000002+00
1,000000000+00
1,000000000+00
1,000000003+00
1,000000002+00
1,000000001+00
1,000000001+00
1,000000001+00
1,000000003+00
1,000000003+00
1,000000002+00
1,000000002+00
1,000000001+00
9,999999921-01
9,999999890-01
1,000000013+00
1,000000010+00
1,000000007+00
1,000000004+00
1,000000001+00
9,999999971-01
9,999999941-01
9,999999913-01
1,000000015+00
1,000000012+00
1,000000009+00
1,000000005+00
1,000000002+00
9,999999993-01
9,999999955-01
9,999999925-01
1,000000017+00
1,000000014+00
1,000000010+00
1,000000007+00
1,000000005+00
1,000000001+00
9,999999978-01
9,999999939-01
9,999999913-01
1,000000016+00
1,000000012+00
1,000000009+00
1,000000006+00
1,000000003+00
9,999999993-01
9,999999961-01
9,999999933-01
1,000000017+00
1,000000013+00
1,000000010+00
1,000000007+00
1,000000004+00
1,000000002+00
9,999999982-01

Gerson W. Barbosa · (This post was last modified: 05-01-2020 06:01 PM by 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

pinkman · 04-30-2020, 09:58 PM

Happy to see you’re still prospecting around Pi and Gamma.
It works on Free42, the accuracy is about 1.E-4 for GX1(.5)

Having a good accuracy (and support for R-/Z-) for Gamma would be of a great help, I should analyze the code and help find the optimisations.

Gerson W. Barbosa · (This post was last modified: 05-01-2020 08:31 PM by 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)

Gerson W. Barbosa · (This post was last modified: 05-01-2020 08:47 PM by Gerson W. Barbosa.)

(04-30-2020 09:58 PM)pinkman Wrote: It works on Free42, the accuracy is about 1.E-4 for GX1(.5)

Having a good accuracy (and support for R-/Z-) for Gamma would be of a great help, I should analyze the code and help find the optimisations.

Please take a look at the latest BASIC version. Conversion to HP-41C is left as an exercise to the reader :-) You can start with the extant RPN code, although it’s not properly optimized as it was just a test. Also, if you use Free42 remember the HP-41C has no recall arithmetic. The lack of that feature might make it difficult to do it all on the stack. Anyway it won’t hurt using a numbered register or two, most important it to keep the program as compact as possible.

Albert Chan · 05-01-2020, 11:59 PM

(04-30-2020 08:35 PM)Gerson W. Barbosa Wrote: Negative integer arguments should return a division by zero error, but because of numerical limitations this won’t occur ...

My guess HP-75C were running with default RADIANS, and HP-71B were on DEGREES
If HP-71B were on RADIANS, (-71.06)! = -1.08421623308E-99, error = 308 - 244 = 64 ULP

To make it work on both angle units, do angle reduction with MOD

10 INPUT X
15 P=1 @ Q=0
20 IF X<0 THEN 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+1/(6-1/(Y+6/(7-6/(Y+C))))
45 F=SQR(2*Z*PI)*(X/EXP(1))^X/P
50 IF Q THEN F=W*PI/(F*SIN(MOD(W,2)*ACOS(-1)))
55 DISP F

>DEFAULT OFF ! turn div-by-0 as error
>RADIANS
>RUN
? -71.06
-1.08421623254E-99
>RUN
? -2
ERR L50:/Zero

Using Sinc function , Euler’s reflection formula is easy to remember: (x!)(-x)! sinc(pi*x) = 1

Gerson W. Barbosa · 05-02-2020, 11:04 AM

(05-01-2020 11:59 PM)Albert Chan Wrote:
(04-30-2020 08:35 PM)Gerson W. Barbosa Wrote: Negative integer arguments should return a division by zero error, but because of numerical limitations this won’t occur ...

My guess HP-75C were running with default RADIANS, and HP-71B were on DEGREES
If HP-71B were on RADIANS, (-71.06)! = -1.08421623308E-99, error = 308 - 244 = 64 ULP

This explains why I was getting different results for that argument. Thanks!
I only noticed it after I posted. I thought of changing the angle mode on the HP-75C but I didn’t remember the syntax is OPTION ANGLES DEGREES, so I would do it later.

(05-01-2020 11:59 PM)Albert Chan Wrote: To make it work on both angle units, do angle reduction with MOD

...

50 IF Q THEN F=W*PI/(F*SIN(MOD(W,2)*ACOS(-1)))

...

>DEFAULT OFF ! turn div-by-0 as error
>RADIANS
>RUN
? -71.06
-1.08421623254E-99
>RUN
? -2
ERR L50:/Zero

That does the trick, but it doesn’t return infinite results for negative odd integer arguments. On the HP-41C I haven’t thought of anything better than multiplying F by FRAC(W)/FRAC(W) to force the division by zero error.

The HP-75C program is just a test for a possible HP-41C version. I have the math module for the HP-71B which includes GAMMA. I wish I had the HP-75C math module, but they appear to be even harder to find. I like the HP-75C because I can place it on a desk and quickly type my programs into it.

(05-01-2020 11:59 PM)Albert Chan Wrote: Using Sinc function , Euler’s reflection formula is easy to remember: (x!)(-x)! sinc(pi*x) = 1

That’s a great mnemonic. I’ll keep it.

Gerson W. Barbosa · (This post was last modified: 05-03-2020 10:05 PM by Gerson W. Barbosa.)

01 LBL "GXP1"
02 1
03 STO 01
04 SF 05
05 X<>Y
06 X<=0?
07 X=0?
08 GTO 01
09 CF 05
10 STO 02
11 CHS
12 LBL 01
13 4
14 STO Z
15 X<=Y?
16 GTO 03
17 +
18 RCL X
19 1
20 +
21 LBL 02
22 1
23 -
24 ST* 01
25 DSE Z
26 GTO 02
27 LBL 03
28 RDN
29 ENTER
30 ENTER
31 41
32 30
33 /
34 RCL Y
35 38
36 /
37 -
38 RCL Y
39 X^2
40 92
41 /
42 +
43 X<>Y
44 ST+ Z
45 ST+ T
46 RDN
47 +
48 6
49 X<>Y
50 /
51 CHS
52 7
53 +
54 6
55 X<>Y
56 /
57 +
58 1/X
59 CHS
60 6
61 +
62 1/X
63 +
64 ST+ X
65 X<>Y
66 R^
67 1
68 E^X
69 /
70 X<>Y
71 Y^X
72 X<>Y
73 PI
74 *
75 SQRT
76 *
77 RCL 01
78 /
79 FS? 05
80 GTO 04
81 PI
82 RCL 02
83 *
84 X<>Y
85 LASTX
86 -1
87 ACOS
88 *
89 SIN
90 *
91 /
92 LBL 04
93 END

Γ(x+1), x = 0..69:

RUN
1,000000000
1,000000001
1,999999999
6,000000014
24,00000001
120,0000001
719,9999998
5.040,000012
40.320,00007
362.880,0005
3.628.800,004
39.916.800,12
479.001.601,3
6.227.020.813,
8,717829135+10
1,307674373+12
2,092278995+13
3,556874291+14
6,402373724+15
1,216451007+17
2,432902018+18
5,109094238+19
1,124000732+21
2,585201681+22
6,204484052+23
1,551121011+25
4,032914629+26
1,088886949+28
3,048883422+29
8,841762138+30
2,652528633+32
8,222838734+33
2,631308388+35
8,683317650+36
2,952327992+38
1,033314794+40
3,719933246+41
1,376375334+43
5,230226253+44
2,039788232+46
8,159152904+47
3,345252679+49
1,405006120+51
6,041526302+52
2,658271563+54
1,196222233+56
5,502622253+57
2,586232450+59
1,241391572+61
6,082818683+62
3,041409332+64
1,551118754+66
8,065817499+67
4,274883375+69
2,308437015+71
1,269640355+73
7,109985964+74
4,052691987+76
2,350561344+78
1,386831189+80
8,320987107+81
5,075802258+83
3,146997390+85
1,982608348+87
1,268869339+89
8,247650677+90
5,443449430+92
3,647111106+94
2,480035546+96
1,711224567+98

Γ(x+1)/x!, x = 0..69:

RUN
1,000000000+00
1,000000001+00
9,999999995-01
1,000000002+00
1,000000000+00
1,000000001+00
9,999999997-01
1,000000002+00
1,000000002+00
1,000000001+00
1,000000001+00
1,000000003+00
1,000000003+00
1,000000002+00
1,000000002+00
1,000000004+00
1,000000003+00
1,000000003+00
1,000000003+00
1,000000002+00
1,000000004+00
1,000000004+00
1,000000004+00
1,000000003+00
1,000000006+00
1,000000005+00
1,000000004+00
1,000000004+00
9,999999921-01
1,000000016+00
1,000000013+00
1,000000010+00
1,000000007+00
1,000000004+00
1,000000001+00
9,999999971-01
9,999999941-01
1,000000018+00
1,000000015+00
1,000000012+00
1,000000009+00
1,000000005+00
1,000000001+00
9,999999993-01
9,999999955-01
1,000000020+00
1,000000017+00
1,000000014+00
1,000000010+00
1,000000007+00
1,000000004+00
1,000000001+00
9,999999978-01
1,000000021+00
1,000000018+00
1,000000016+00
1,000000012+00
1,000000009+00
1,000000006+00
1,000000003+00
9,999999993-01
1,000000023+00
1,000000020+00
1,000000017+00
1,000000013+00
1,000000010+00
1,000000007+00
1,000000004+00
1,000000002+00
1,000000025+00

Γ(x+1), x = -19/2..17/2

RUN
-2,633521509-05
2,238493289-04
-1,678869965-03
1,091265477-02
-6,001960120-02
2,700882051-01
-9,453087196-01
2,363271799+00
-3,544907695+00
1,772453847+00
8,862269277-01
1,329340391+00
3,323350974+00
1,163172841+01
5,234277792+01
2,878852784+02
1,871254308+03
1,403440731+04
1,192924620+05

Edited to change lines 06 and 07.

Previously,

06 X!=0?
07 X>0?

Gerson W. Barbosa · (This post was last modified: 05-09-2020 05:29 PM by Gerson W. Barbosa.)

.
HP-41C:

01 LBL "GXP1"
02 1
03 STO 01
04 SF 05
05 X<>Y
06 X<=0?
07 X=0?
08 GTO 01
09 CF 05
10 STO 02
11 CHS
12 LBL 01
13 4
14 STO Z
15 X<=Y?
16 GTO 03
17 +
18 RCL X
19 1
20 +
21 LBL 02
22 1
23 -
24 ST* 01
25 DSE Z
26 GTO 02
27 LBL 03
28 X<>Y
29 4.13333
30 2.215
31 RCL Z
32 *
33 .4875
34 -
35 R^
36 X^2
37 /
38 +
39 12
40 /
41 +
42 72
43 *
44 1/X
45 6
46 1/X
47 +
48 +
49 360
50 D-R
51 *
52 SQRT
53 X<>Y
54 1
55 E^X
56 /
57 R^
58 Y^X
59 *
60 RCL 01
61 /
62 FS? 05
63 GTO 04
64 PI
65 RCL 02
66 *
67 X<>Y
68 LASTX
69 -1
70 ACOS
71 *
72 SIN
73 *
74 /
75 LBL 04
76 END

117 BYTES

———————-

-69.95 -> -1.450781639E-97 {0 ULP}
-2.5 -> 2.36327180(1) [4/3×√π] {1}
-2.0 -> DATA ERROR
-1.5 -> -3.544907(697) [-2√π] {5}
-1.0 -> DATA ERROR
-0.5-> 1.7724538(47) [√π] {4}
0 -> 1.00000000(1) {1}
.5-> .88622692(77) [√(π/4)] {22}
1 -> 1.00000000(1) {1}
2 -> 1.999999999 {1}
3 -> 6.0000000(12) {12}
4 -> 24.0000000(3) {3}
5 -> 119.999999(5) {5}
6 -> 719.999999(5) {5}
7 -> 5040.0000(10) {10}
8 -> 40320.0000(6) {6}
9 -> 362880.000(4) {4}
10-> 3628800.00(2) {2}
11-> 39916800.(11) {11}
12-> 47900160(1.1) {11}
13-> 62270208(12) {12}
14-> 87178291(32) {12}
15-> 1.3076743(73)E12 {5}
20-> 2.4329020(18)E18 {10}
30-> 2.652528(633)E32 {35}
50-> 3.0414093(32)E64 {12}
60-> 8.3209871(07)E81 {6}
69-> 1.7112245(67)E98 {43}
69.95 -> 9.682847673E99

———————-

HP-75C:

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 C=62/15-443/(200*X)-39/(80*X*X)
35 Z=(6*X+1+1/(12*X+C))/6
40 F=SQR(2*Z*PI)*(X/EXP(1))^X/P
45 IF Q THEN F=W*PI/(F*SIN(W*ACOS(-1)))
50 D=DISP F

————-

-253.11 -> -2468119655(10)E-497 {2 ULP}
-71.06 -> -1.084216232(55)E-99 {11}
-2.5 -> 2.363271801(91) [4/3×√π] {70}
-2.0 -> -9.99999999999E499 [-∞]
-1.5 -> -3.54490770(200) [-2√π] {19}
-1.0 -> -9.99999999999E499 [-∞]
-0.5-> 1.772453850(91) [√π] {76}
0 -> 1.0000000000(3) {3}
.5-> .886226925(833) [√(π/4)] {380}
1 -> 1.0000000000(2) {2}
2 -> 1.999999999(61) {39}
3 -> 5.99999999(801) {199}
4 -> 24.00000000(73) {73}
5 -> 120.0000000(24) {24}
6 -> 719.999999(860) {140}
7 -> 5039.99999(833) {167}
8 -> 40319.9999(865) {135}
9 -> 362879.999(892) {108}
10-> 3628799.999(06) {94}
11-> 39916799.99(15) {85}
12-> 479001599.9(26) {74}
13-> 6227020799.(21) {79}
14-> 87178291(190.7) {93}
15-> 1.30767436(787)E12 {13}
20-> 2.432902008(10)E18 {8}
30-> 2.65252859(778)E32 {34}
50-> 3.041409320(40)E64 {23}
60-> 8.32098711(301)E81 {27}
69-> 1.711224524(14)E98 {14}
69.95 -> 9.68284767(192)E99 {108}
253.1 -> 8.998861511(22)E499 {23}

ANGLE OPTION DEGREES

For better accuracy in radians mode change line 45 to

45 IF Q THEN F=W*PI/(F*SIN(MOD(W,2)*ACOS(-1)))

per Albert Chan’s suggestion above.

For best 12-digit results change line 30 to

30C=4.1332883071+2.1999058932/X-.3505828048/(X*X)-.2950235942/X^3
————-

Edited to fix a few typos.

Albert Chan · (This post was last modified: 09-11-2020 01:28 PM by Albert Chan.)

I tried turning Stirling's formula's correction to Gamma, in continued fraction form.
Amazingly, every divide gives me back 2 terms.

XCas> c3(x) := 1 + 1/(12x-1/2+1/(720/293*x+1/(7211316/4406147*x)))

\(x \rightarrow 1+\frac{1}{\Large 12x- \frac{1}{2}+\frac{1}{\frac{720}{293}x+\frac{1}{\frac{7211316}{4406147}x}}}\)

XCas> series(c3(x), x=inf, polynom)

\(1+\frac{1}{12}\left({1\over x}\right)
+\frac{1}{288} \left({1\over x}\right)^2
-\frac{139}{51840} \left({1\over x}\right)^3
-\frac{571}{2488320} \left({1\over x}\right)^4
+\frac{163879}{209018880} \left({1\over x}\right)^5
+\frac{5246819}{75246796800} \left({1\over x}\right)^6
\)

All terms matches correctly to Series of Gamma(x) / (sqrt(2*pi/x) * (x/e)^x), x=inf

Below, we define 3 functions, FNS(x)=sinc(pi*x) , FNF(n)=n! , FNG(x)=Γ(x)

Code:

10 DEF FNS(X)=SIN(ACOS(-1)*MOD(X,2))/(PI*X) ! = sinc(pi*x)

20 DEF FNF(N)                               ! = factorial

30 IF N<0 THEN FNF=1/(FNS(N)*FNG(1-N)) ELSE FNF=FNG(1+N)

40 END DEF

100 DEF FNG(X)                              ! = gamma

110 C=1 @ WHILE X<12 @ C=C*X @ X=X+1 @ END WHILE

120 C=SQRT(2*PI)*X^(X-.5)*EXP(-X)/C         ! stirling

130 FNG=C+C/(12*X-.5+1/(720/293*X+.609/X))  ! correction

140 END DEF

>RUN
>FNF(5), FNF(10), FNF(15)
120 3628800 1.307674368E12

>X=1.1 ! check reflection formula (FNG does not do reflection)
>FNF(X)*FNF(-X)*FNS(X), FNG(1+X)*FNG(1-X)*FNS(X)
1 1

>FOR X=1 TO 2 STEP .1 @ G=FNG(X) @ X,G,GAMMA(X)-G @ NEXT X
1 1 0
1.1 .951350769865 .000000000002
1.2 .918168742399 .000000000001
1.3 .897470696308 -.000000000002
1.4 .887263817504 -.000000000001
1.5 .886226925454 -.000000000001
1.6 .893515349285 .000000000003
1.7 .908638732849 .000000000004
1.8 .931383770985 -.000000000005
1.9 .961765831906 .000000000001
2 1 0

Albert Chan · 09-13-2020, 12:49 PM

We can easily get reciprocal of continued fraction (here, b is rest of the CF terms).

\(\large 1 ÷ \left( 1 + \Large {1 \over (a-{1\over2})\;+\;b\;} \right)
= \frac{(a-{1\over2})\;+\;b}{(a+{1\over2})\;+\;b}
= \large 1 - \Large {1 \over (a+{1\over2})\;+\;b\;}
\)

Example, this is the code for 1/Γ(x), using correction 1/c3(x)

Code:

50 DEF FNR(X)                              ! = 1/gamma

60 C=1 @ WHILE X<12 @ C=C*X @ X=X+1 @ END WHILE

70 C=C/(SQRT(2*PI)*X^(X-.5)*EXP(-X))       ! 1/stirling

80 FNR=C-C/(12*X+.5+1/(720/293*X+.609/X))  ! correction

90 END DEF