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

[Gerson W. Barbosa]

and about half the memory space.

Code:


01- ENTER 
02- ENTER 
03- ENTER 
04- 8
05- 4
06- 0
07- ×
08- 3  
09- 1
10- 4
11- +
12- ×
13- 6
14- 6
15- + 
16- STO 1
17- R↓
18- 5
19- 0
20- 4
21- 0
22- ×
23- 1
24- 4
25- 6
26- 4
27- +
28- ×
29- 4 
30- 1
31- 9
32- +
33- STO/ 1
34- R↓
35- RCL 1
36- +
37- 7
38- 1
39- 0
40- ×
41- 1
42- 1
43- 3
44- /
45- √
46- R↓
47- R↓
48- LN
49- 1
50- -
51- ×
52- eᵡ
53- ×
54- GTO 00

0.1 R/S -> 0.9(392743098)
0.5 R/S -> 0.8862(420073)
1 R/S -> 1.0000(93310)
2 R/S -> 2.0000(31342)
3 R/S -> 6.0000(24847)
4 R/S -> 24.0000(3634)
5 R/S -> 120.0000(823)
10 R/S -> 3628800.(316)
12 R/S -> 4790016(319)
13 R/S -> 622702(1117)
14 R/S -> 8.717829(639) 10
15 R/S -> 1.307674(440) 12
20 R/S -> 2.432902(140) 18
50 R/S -> 3.041409(537) 64
60 R/S -> 8.3209871(88) 81
69 R/S -> 1.7112245(67) 98
69.95 R/S -> 9.68284767(3) 99

[Albert Chan]

Very nice and compact formula, where did you get it ?
I am curious, your thread named half-precision, does a full-precision version exist ?

I rewrite below term, reduced 9 steps, all done in the stack.

\(\large c={840 x^2 + 314x + 66 \over 5040x^2 + 1464x + 419} =
{1 \over {6 - {1 \over \Large 2x + {360x+66 \over 420x-23} }}}\)

\(\Gamma(x+1) ≈ (x/e)^x \sqrt{2\pi(x+c)}\)

Code:

01- ENTER 
02- ENTER 
03- ENTER 
04- 3
05- 6
06- 0
07- ×
08- 6  
09- 6
10- +
11- X<>Y
12- 4
13- 2
14- 0
15- ×
16- 2
17- 3
18- − 
19- /
20- +
21- +
22- 1/X
23- CHS
24- 6
25- +
26- 1/X
27- +
28- 7
29- 1
30- 0
31- ×
32- 1
33- 1
34- 3
35- /
36- √
37- R↓
38- R↓
39- LN
40- 1
41- -
42- ×
43- eᵡ
44- ×
45- GTO 00

[Gerson W. Barbosa]

(02-16-2020 05:47 PM)Albert Chan Wrote:  Very nice and compact formula, where did you get it ?
I am curious, your thread named half-precision, does a full-precision version exist ?

I rewrite below term, reduced 9 steps, all done in the stack.

\(\large c={840 x^2 + 314x + 66 \over 5040x^2 + 1464x + 419} =
{1 \over {6 - {1 \over \Large 2x + {360x+66 \over 420x-23} }}}\)

\(\Gamma(x+1) ≈ (x/e)^x \sqrt{2\pi(x+c)}\)

That was just a test. I have seen correction factors to Stirling's approximation outside the surd, so I decided to place a continued fraction inside it, hoping it to eventually assume a generalized form. If such an infinite continued fraction existed, then it would be easy to get full accuracy with minimum memory for whatever programmable calculator or computer we would like to implement it on. Anyway, here is what I have so far, in case someone is interested in finding the elusive pattern, if any:

\[x! \approx \sqrt{2\pi \left ( x + \frac{1}{6-\frac{1}{2x+\frac{6}{7-\frac{5}{3x+\frac{11}{20}}}}} \right )}\left ( \frac{x}{e} \right )^{x}\]

Here, the last denominator (11/20) is somewhat arbitrary.

I had written it using the plain continued fraction above, in 47 steps or so, but since it's not definitive, I decided to use the Horner form version. In the linked Wikipedia article there is an even more compact and more accurate formula. Scroll down to "Versions suitable for calculators" and look for Gergő Nemes's formula.

Back to the continued fraction, that was the latest one I was trying yesterday:

\[x + \frac{1}{6-\frac{1}{2x+\frac{6}{7-\frac{5}{3x+\frac{11}{8-\frac{9}{4x+\frac{16}{23-\frac{13}{5x+\frac{21}{10}}}}}}}}}\]

I think the first four or five terms are correct. The next terms don't appear to improve the approximation, so they are partially or completely wrong.

[Gerson W. Barbosa]

More than two-third precision and one-third the program memory space using Gergő Nemes's formula mentioned above:

Code:


01- 1
02- +
03- ENTER 
04- ENTER 
05- 12×
06- x⇆y
07- 1
08- 0
09- ×
10- 1/x
11- -
12- 1/x
13- +
14- 1
15- eᵡ
16- /
17- x⇆y
18- yᵡ
19- x⇆y
20- 3
21- 9
22- 6
23- LN
24- 8
25- ×
26- 5
27- 8
28- √
29- /
30- x⇆y
31- /
32- √
33- ×
34- GTO 00

0.1 R/S -> 0.951(1141433)
0.5 R/S -> 0.886(1706366)
1 R/S -> 0.9999(830883)
2 R/S -> 1.99999(5097)
3 R/S -> 5.99999(6379)
4 R/S -> 23.99999(521)
5 R/S -> 119.99999(02)
10 R/S -> 3628799.99(3)
12 R/S -> 479001599.(2)
13 R/S -> 6227020(775)
14 R/S -> 8.7178291(15) 10
15 R/S -> 1.307674368 12
20 R/S -> 2.4329020(19) 18
50 R/S -> 3.0414093(04) 64
60 R/S -> 8.32098711(0) 81
69 R/S -> 1.7112245(50) 98
69.67 R/S -> 2.9433192(47) 99

[Gamo]

This version used two Registers: R0 and R1

Store Pi to R0

Γ(x) for 1 ≤ x ≤ 64

---------------------

Example: [FIX] 2

Γ(5.25)
5.25 [R/S] 35.21

7! = Γ(8)
8 [R/S] 5040.00

2.34! = Γ(3.34)
3.34 [R/S] 2.80
--------------------
Program:

Code:


01  5
02  +
03 ENTER
04 1/x
05 ENTER
06  x
07 ENTER
08 ENTER
09  3
10  0
11  ÷
12  1
13  -
14  x
15  1
16  2
17  ÷
18 X<>Y
19 LN
20  -
21  x
22  +
23 CHS
24 e^x
25 X<>Y
26 ENTER
27  +
28 RCL 0
29  x
30 √x
31  x
32 STO 1
33 CLx
34  5
35  -
36 STO÷1
37  1
38  +
39 X≤Y
40 GTO 36
41 RCL 1
42 GTO 00

Remark: This routine is adapted from HP-55 Statistic Programs P.14

Gamo

[Gerson W. Barbosa]

(02-16-2020 07:39 PM)Gerson W. Barbosa Wrote:  If such an infinite continued fraction existed, then it would be easy to get full accuracy with minimum memory for whatever programmable calculator or computer we would like to implement it on.

Most likely there is no such a continued fraction, otherwise it would be known by now. If it existed short programs for gamma function would be possible, like the one below for the HP-75C. The results won’t become any better by increasing the number of iterations (I).

10 INPUT X
15 Y=2*X
20 I=4
25 N=5*I-4
30 D=I+5
35 C=Y+(N+5)/(D+1)
40 FOR K=1 TO I
45 C=Y+N/(D-N/C)
50 N=N-5
55 D=D-1
60 NEXT K
70 F=SQR(2*PI*(C-X))*(X/EXP(1))^X
75 DISP F

.5-> .88662(7414717)
1 -> 1.0000(5537894)
2 -> 2.00000(622165)
3 -> 6.00000(235186)
4 -> 24.00000(15634)
5 -> 120.000000(852)
6 -> 719.999996708 (720)
7 -> 5039.99996836 (5040)
8 -> 40319.9997698 (40320)
9 -> 362879.998278 (362880)
10-> 3628799.98602 (3628800)
11-> 39916799.8757 (39916800)
12-> 479001598.800 (479001600)
13-> 6227020787.22 (6227020800)
14-> 87178291052.2 (87178291200)
15-> 1.30767436(614)E12
20-> 2.43290200(663)E18
50-> 3.041409320(28)E64
60-> 8.320987112(84)E81
69-> 1.711224524(11)E98


P.S.: The results will be consistently less then the exact ones (except when x is close to 0) if line 20 is changed to

20 I=2

[Gerson W. Barbosa]

41 steps on my HP-41C (but you can always make it shorter).

01 LBL "GXP1"
02 ENTER
03 ENTER
04 336
05 *
06 224
07 +
08 X<>Y
09 *
10 99
11 ST+ Y
12 X<> L
13 ENTER
14 ENTER
15 56
16 *
17 42
18 +
19 X<>Y
20 *
21 18
22 ST+ Y
23 X<> L
24 X<> Z
25 /
26 X<>Y
27 +
28 ST+ X
29 LASTX
30 RCL X
31 -1
32 E^X
33 *
34 X<>Y
35 Y^X
36 X<>Y
37 PI
38 *
39 SQRT
40 *
41 END



.5-> 0.8862(859731)
1 -> 1.00000(8281)
2 -> 2.00000(1930)
3 -> 6.000000(983)
4 -> 24.000000(68)
5 -> 120.0000000
6 -> 719.9999952 (720)
7 -> 5039.999965 (5040)
8 -> 40319.99985 (40320)
9 -> 362879.9989 (362880)
10-> 3628799.990 (3628800)
11-> 39916799.89 (39916800)
12-> 479001598.8 (479001600)
13-> 6227020801 (6227020800)
14-> 871782912(1)E10
15-> 1.307674368E12
20-> 2.4329020(11)E18
50-> 3.0414093(32)E64
60-> 8.3209871(07)E81
69-> 1.71122452(1)E98
69.95 -> 9.68284767(8)E99

P.S.:

37 steps

01 LBL "GXP1"
02 ENTER
03 ENTER
04 +
05 RCL X
06 1.5
07 +
08 6
09 X<>Y
10 /
11 CHS
12 7
13 +
14 6
15 X<>Y
16 /
17 +
18 1/X
19 CHS
20 6
21 +
22 1/X
23 +
24 ST+ X
25 X<>Y
26 R^
27 -1
28 E^X
29 *
30 X<>Y
31 Y^X
32 X<>Y
33 PI
34 *
35 SQRT
36 *
37 END

[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

[Gerson W. Barbosa]

Very good job! You can call it full-precision Γ(x+1) now.

[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

[Gerson W. Barbosa]

(09-13-2020 10:29 PM)Albert Chan Wrote:  Another approach is not to use asymptotic formula at all, and use Lanczos algorithm.

Lanczos approximation is great. The problem is the constants take up too much memory space on the 12C and even on the HP-41C, I fear.

[bshoring]

(02-16-2020 12:30 AM)Gerson W. Barbosa Wrote:  and about half the memory space.

Code:


01- ENTER 
02- ENTER 
03- ENTER 
04- 8
05- 4
06- 0
07- ×
08- 3  
09- 1
10- 4
11- +
12- ×
13- 6
14- 6
15- + 
16- STO 1
17- R↓
18- 5
19- 0
20- 4
21- 0
22- ×
23- 1
24- 4
25- 6
26- 4
27- +
28- ×
29- 4 
30- 1
31- 9
32- +
33- STO/ 1
34- R↓
35- RCL 1
36- +
37- 7
38- 1
39- 0
40- ×
41- 1
42- 1
43- 3
44- /
45- √
46- R↓
47- R↓
48- LN
49- 1
50- -
51- ×
52- eᵡ
53- ×
54- GTO 00

Thanks for this. I’m having fun with it on my HP-38C emulator (RPN-38CX). Here’s my listing. If anyone has RPN-38CX, you can paste this directly into the display in program mode. By the way, I changed storage register 1 to I. That way no numbered or Financial registers are used.

[code]
01 - 31 ENTER
02 - 31 ENTER
03 - 31 ENTER
04 - 8 8
05 - 4 4
06 - 0 0
07 - 61 ×
08 - 3 3
09 - 1 1
10 - 4 4
11 - 51 +
12 - 61 ×
13 - 6 6
14 - 6 6
15 - 51 +
16 - 21 26 STO I
17 - 25 33 R↓
18 - 5 5
19 - 0 0
20 - 4 4
21 - 0 0
22 - 61 ×
23 - 1 1
24 - 4 4
25 - 6 6
26 - 4 4
27 - 51 +
28 - 61 ×
29 - 4 4
30 - 1 1
31 - 9 9
32 - 51 +
33 - 21 71 26 STO ÷ I
34 - 25 33 R↓
35 - 22 51 26 RCL + I
36 - 7 7
37 - 1 1
38 - 0 0
39 - 61 ×
40 - 1 1
41 - 1 1
42 - 3 3
43 - 71 ÷
44 - 24 21 √x
45 - 25 33 R↓
46 - 25 33 R↓
47 - 25 23 LN
48 - 1 1
49 - 41 −
50 - 61 ×
51 - 25 22 eˣ
52 - 61 ×
53 - 25 7 00 GTO 00
[code]