# HP Forums

Full Version: What is the Gamma approximation you use?
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

If you calculate the gamma function and have to use an approximation (instead of built in functions like with Excel and Matlab), what approximation would you use?

What does the WP34S use for the Gamma function?

Namir

Quote:
What does the WP34S use for the Gamma function?

Just look at it - it's open source.

d;-)

If you haven't looked at the source before you can grab your
own copy with this command:

svn checkout svn://svn.code.sf.net/p/wp34s/code wp34s

I took a look at the code. It is in wp34s/trunk/decn.c file.

It seems to be using a function that computes the Natural Log of Gamma and then raises e to that power.

The function seems similar to the Spouge Gamma Approximation using
21 precomputed constants. It is a very accurate approximation
(seems better than most desktop math library functions that I have compared it to.)

Edited: 2 Aug 2013, 7:05 p.m.

I made a small survey for the most popular approximations for the gamma function. The Spouge method came first. You can calculate the constants somewhat easily and dynamically, unlike many other approximations. The Lanczos approximation appearing in Numerical Recipes. This method required storing an array of constants.

Namir

Edited: 2 Aug 2013, 10:01 p.m. after one or more responses were posted

Quote:
The Spouge method came first. You can calculate the constants somewhat easily and dynamically, unlike many other approximations.

Like Lanczos, does it work in the complex domain also? I
recently used the Lancoz Approximation on the HP-48G/GX and HP-28S and was pleased with the result.

Gerson.

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

Gamma:
%%HP: T(3)A(D)F(.);
\<< 0 { 76.1800917295
-86.5053203294
24.0140982408
-1.23173957245
1.20865097387E-3
-5.39523938495E-6 }
1
\<< 3 PICK NSUB + /
+
\>> DOSUBS
1.00000000019 +
OVER / 2 \pi * \v/ *
SWAP 5.5 + DUP LN
OVER 5 - * SWAP -
EXP *
\>>
Gamma:
« 0 { 76.1800917295
-86.5053203294
24.0140982408
-1.23173957245
1.20865097387E-3
-5.39523938495E-6 }
1 6
FOR i DUP i GET 4
PICK i + / ROT +
SWAP
NEXT DROP
1.00000000019 + OVER
/ 2 ‡ * ƒ * SWAP 5.5
+ DUP LN OVER 5 - *
SWAP - EXP *
»
‡ --> pi
ƒ --> sqrt


P.S.: Never mind! The answer to my question can be found in Pugh's thesis (linked in Paul Dale's post below), starting at page 36.

Edited: 2 Aug 2013, 9:15 p.m.

I'm not certain which approximation we use now, I think it is Lanczos. Marcus did some work to improve the accuracy of gamma a year or two back and I don't remember if the algorithm was changed or just the constants. The history is in subversion if anyone really wants to check.

The reference I used for the first implementation was Pugh's thesis on the gamma function. I originally used the table of constants on page 126 which we later found out weren't entirely correct.

the algorithm used works for real and complex arguments using the same series.

- Pauli

Quote:
Just look at it - it's open source.

Well that sure helps the sharing and debating nature of this forum...

I used a Lanczos approximation with 6 terms, for both the SandMath and 41Z implementations. It works fine (accurate to the 9th decimal digit for real arguments and the 8th for complex at worst) - with the reduced precision range in the platform but you guys are in the stratospheric accuracy range so I'm sure have used more terms or yet a better approach.

Namir

The algorithim employed varies w/ the artifact employed (different routines for different machines). I find the assorted routines on Viktor T. Toth's web site

http://www.rskey.org/CMS/index.php/exhibit-hall/95

very interesting: thoughts?

ps: view the page for a model to see the attendant gamma routine for that model.

SlideRule

Here are the source codes for the gamma function using Spouge's approximation for the HP-71B, HP41C, and HP-67:

HP-71B Implementation
======================
10 DESTROY ALL
20 REM GAMMA USING SPOUGE APPROX
30 INPUT "ENTER X? ";X @ X=X-1
40 A=12.5 @ P2=SQR(2*PI)
50 S=1 @ S1=1
60 FOR I= 1 TO 12
70 C=S1/P2/FACT(I-1)*(A-I)^(I-0.5)*EXP(A-I)
80 S=S+C/(X+I)
90 S1=-S1
100 NEXT I
110 X1=X+A
120 G=X1^(X+.5)/EXP(X1)*P2*S
130 DISP "GAMMA=";G
HP-41-C Implementation
======================
R00 = x and = x-1
R01 = a
R02 = CHS
R03 = Sum
R04 = I
R05 = sqrt(2*pi)
R06 = integer part of I
LBL "GAMMA"
LBL A
"X?"
PROMPT
1
-
STO 00
12.5
STO 01	# a = 12.5
1
STO 02	# CHS = 1
STO 03	# Sum = 1
2
PI
*
SQRT
STO 05	# sqrt(2*pi)
1.012
STO 04	# set up loop control variable
LBL 00	# start the loop
RCL 02
RCL 05
/
RCL 04
INT
STO 06
1
-
FACT
/
RCL 01
RCL 06
-
RCL 06
0.5
-
Y^X
*
RCL 01
RCL 07
-
EXP
*
RCL 06
RCL 00
+
/
STO+ 03	# Sum = Sum + C/(X+I)
RCL 02
CHS
STO 02	# CHS = -CHS
ISG 04	# end of loop
GTO 00
RCL 00
RCL 01
+
STO 06
RCL 00
0.5
+
Y^X
RCL 06
EXP
/
RCL 05
*
RCL 03
*
"GAMA="
ARCL X
PROMPT
GTO A
HP-67 Implementation
====================
R0 = x and = x-1
R1 = a
R2 = CHS
R3 = Sum
R4 = Integer part of I, x+a
R5 = sqrt(2*pi)
RI = I
LBL A
1
-
STO 0
12.5
STO 1	# a = 12.5
1
STO 3	# Sum = 1
CHS
STO 2	# CHS = -1
2
PI
*
SQRT
STO 5	# sqrt(2*pi)
12
CHS
STI	# set up loop control variable
LBL 0	# start the loop
RCL 2
RCL 5
/
RCI
ABS
STO 4
1
-
N!
/
RCL 1
RCL 4
-
RCL 4
0.5
-
Y^X
*
RCL 1
RCL 4
-
EXP
*	# calculate C
RCL 4
RCL 0
+
/
STO+ 3	# Sum = Sum + C/(X+I)
RCL 2
CHS
STO 2	# CHS = -CHS
ISZ	# end of loop
GTO 0
RCL 0
RCL 1
+
STO 4
RCL 0
0.5
+
Y^X
RCL 4
EXP
*
RCL 5
*
RCL 3
*
R/S
GTO A


In the case of the 67 and 41C, enter the value for x and press the [A] key to get the gamma function value.

Namir

Edited: 3 Aug 2013, 10:42 a.m.

THANKS!

ps: I am indebted to you for ALL your marvelous postings (here & your web page)!

Many thanks, again!

SlideRule

Here is a non-optimal HP-48G/GX version:

%%HP: T(3)A(D)F(.);
\<< 1 - 12.5 \pi 2 * \v/
1 \-> x a p s
\<< 1 1 12
FOR i '(-1)^(i+
1)/p/(i-1)!*(a-i)^(
i-.5)*EXP(a-i)'
EVAL x i + / +
NEXT a x + DUP
x .5 + ^ SWAP EXP /
* p *
\>>
\>>
(1, 2)  -->  (.151904002653, 198048801563E-2)

As a comparison, the HP-50g GAMMA function returns
             (.15190400267,  1980488015619E-2)


Very nice!

P.S.: The local variable s is not necessary. The following should be slightly faster:

%%HP: T(3)A(D)F(.);
\<< 1 - 12.5 \pi 2 * \v/
\-> x a p
\<< 1 1 12
FOR i '1/p/(i-1
)!*(a-i)^(i-.5)*EXP
(a-i)' EVAL x i + /
+ NEG
NEXT a x + DUP
x .5 + ^ SWAP EXP /
* p *
\>>
\>>


Edited: 3 Aug 2013, 1:43 p.m.

You are most welcome. Sharing code here is fun. You can find more code for calculators and some programming languages on my web site. Please click here.

Only 8 terms appear to give more accurate results on the HP-71B and HP-48:

10 DESTROY ALL
20 DISP "    GAMMA(N)         N"
30 FOR N=1 TO 15 @ X=N-1
40 A=8.5 @ P2=SQR(2*PI)
50 S=1
60 FOR I=1 TO 8
70 C=1/P2/FACT(I-1)*(A-I)^(I-.5)*EXP(A-I)
80 S=-S-C/(X+I)
90 NEXT I
100 X1=X+A
110 G=X1^(X+.5)/EXP(X1)*P2*S
120 DISP G,
130 DISP USING 150;N;
140 NEXT N
150 IMAGE DD
RUN
GAMMA(N)         N
1.00000000002        1
1.00000000006        2
2.00000000014        3
6.00000000051        4
24.0000000017        5
120.000000010        6
720.000000110        7
5040.00000097        8
40320.0000039        9
362880.000055       10
3628800.00049       11
39916800.0045       12
479001600.110       13
6227020799.76       14
87178291207.6       15

When the lines 40 and 60 are changed to
40 A=12.5 @ P2=SQR(2*PI)
60 FOR I=1 TO 12

the output is
    GAMMA(N)         N
0.99999999999        1
1.00000000015        2
1.99999999980        3
6.00000000551        4
24.0000000236        5
120.000000091        6
720.000000533        7
5040.00001102        8
40320.0000823        9
362880.000676       10
3628800.00388       11
39916800.1351       12
479001600.603       13
6227020811.48       14
87178291351.0       15


Very interesting. Spouge's methods calculates the upper limit of the summation (that needs the FOR loop) as the integer(ceiling(A))-1 which gives 12 for A=12.5.

Using an upper limit of 8 must be causing the accuracy of the 71B to give better results. I will keep that in mind! Thanks!

Namir

I took a look at your page also. Bookmarked it in fact. I even took a couple things. Fantastic site, thanks so much.

Pauli is right, it's some time ago when I had a closer look at Gamma to make it fit for double precision. I used Pugh's thesis and information from Victor T. Toth's site. The list of constants can be found in the file compile_consts.c, but here you are:

// Gamma estimate constants
{ DFLT, "gammaR",	"23.118910" },
{ DFLT, "gammaC00",	"2.5066282746310005024157652848102462181924349228522"},
{ DFLT, "gammaC01",	"18989014209.359348921215164214894448711686095466265"},
{ DFLT, "gammaC02",	"-144156200090.5355882360184024174589398958958098464"},
{ DFLT, "gammaC03",	"496035454257.38281370045894537511022614317130604617"},
{ DFLT, "gammaC04",	"-1023780406198.473219243634817725018768614756637869"},
{ DFLT, "gammaC05",	"1413597258976.513273633654064270590550203826819201"},
{ DFLT, "gammaC06",	"-1379067427882.9183979359216084734041061844225060064"},
{ DFLT, "gammaC07",	"978820437063.87767271855507604210992850805734680106"},
{ DFLT, "gammaC08",	"-512899484092.42962331637341597762729862866182241859"},
{ DFLT, "gammaC09",	"199321489453.70740208055366897907579104334149619727"},
{ DFLT, "gammaC10",	"-57244773205.028519346365854633088208532750313858846"},
{ DFLT, "gammaC11",	"12016558063.547581575347021769705235401261600637635"},
{ DFLT, "gammaC12",	"-1809010182.4775432310136016527059786748432390309824"},
{ DFLT, "gammaC13",	"189854754.19838668942471060061968602268245845778493"},
{ DFLT, "gammaC14",	"-13342632.512774849543094834160342947898371410759393"},
{ DFLT, "gammaC15",	"593343.93033412917147656845656655196428754313318006"},
{ DFLT, "gammaC16",	"-15403.272800249452392387706711012361262554747388558"},
{ DFLT, "gammaC17",	"207.44899440283941314233039147731732032900399915969"},
{ DFLT, "gammaC18",	"-1.2096284552733173049067753842722246474652246301493"},
{ DFLT, "gammaC19",	".0022696111746121940912427376548970713227810419455318"},
{ DFLT, "gammaC20",	"-.00000079888858662627061894258490790700823308816322084001"},
{ DFLT, "gammaC21",	".000000000016573444251958462210600022758402017645596303687465"},


To get better results:

- do loop the other way (8 downto 1) starting with small factors to keep accuracy

- do sum two at a time to avoid loss of digits too early in the loop

10 DESTROY ALL
20 DISP "    GAMMA(N)         N"
30 FOR N=1 TO 19 @ X=N-1
40 A=8.5 @ P2=SQR(2*PI)
50 S=1 @ S0=0
60 FOR I=1 TO 8
70 C=1/P2/FACT(I-1)*(A-I)^(I-.5)*EXP(A-I)
72 S=-S-C/(X+I)
74 IF MOD(I,2)=0 THEN 90
76 I0=9-I @ I1=I0-1
78 C0=(A-I0)^(I0-.5)*EXP(A-I0)/FACT(I0-1)
80 C1=(A-I1)^(I1-.5)*EXP(A-I1)/FACT(I1-1)
82 C1=C1/(X+I1) @ C0=C0/(X+I0) @ C2=C1-C0
84 S0=S0+C2
86 IF I=7 THEN S0=S0+P2
90 NEXT I
100 X1=X+A
110 G=X1^(X+.5)/EXP(X1)*P2*S
112 G0=X1^(X+.5)/EXP(X1)*S0
120 DISP G, G0
140 NEXT N
150 IMAGE DD


It gives better results most of the time (rounding can be tricky)

 1.00000000002        1.00000000001        1
1.00000000006        1.00000000002        1
2.00000000014        2.00000000015        2
6.00000000051        6.00000000056        6
24.0000000017        24.0000000019        24
120.00000001         120.000000002        120
720.00000011         720.000000111        720
5040.00000097        5040.00000083        5040
40320.0000039        40320.0000025        40320
362880.000055        362879.999993        362880
3628800.00049        3628800.00034        3628800
39916800.0045        39916800.0035        39916800
479001600.11         479001600.007        479001600
6227020799.76        6227020800.96        6227020800
87178291207.6        87178291211.9        87178291200
1.30767436817E12     1.30767436809E12     1.30767436800E12
2.09227898902E13     2.09227898901E13     2.09227898880E13
3.55687428249E14     3.55687428098E14     3.55687428096E14
6.40237370629E15     6.40237370589E15     6.40237370573E15


Olivier

Edited: 4 Aug 2013, 11:28 a.m.

It appears the left column results are more accurate, but perhaps the sample is too small. You might want to expand it. K=8 in line 35 gives more exact answer when compared to other even K.

Regards,

Gerson.

10 DESTROY ALL
15 S1=0 @ S2=0
20 DISP TAB(15);"GAMMA(N)                    N"
30 FOR N=1 TO 14 @ X=N-1
35 K=8
40 A=K+.5 @ P2=SQR(2*PI)
50 S=1 @ S0=0
60 FOR I=1 TO K
70 C=1/P2/FACT(I-1)*(A-I)^(I-.5)*EXP(A-I)
72 S=-S-C/(X+I)
74 IF MOD(I,2)=0 THEN 90
76 I0=K+1-I @ I1=I0-1
78 C0=(A-I0)^(I0-.5)*EXP(A-I0)/FACT(I0-1)
80 C1=(A-I1)^(I1-.5)*EXP(A-I1)/FACT(I1-1)
82 C1=C1/(X+I1) @ C0=C0/(X+I0) @ C2=C1-C0
84 S0=S0+C2
86 IF I=K-1 THEN S0=S0+P2
90 NEXT I
100 X1=X+A
110 G=X1^(X+.5)/EXP(X1)*P2*S
112 G0=X1^(X+.5)/EXP(X1)*S0
115 S1=S1+G @ S2=S2+G0
120 DISP G,G0,
130 DISP USING 150;N
140 NEXT N
145 DISP ABS(S1-6749977114),ABS(S2-6749977114)
150 IMAGE DD
RUN
GAMMA(N)                    N
1.00000000002        1.00000000001        1
1.00000000006        1.00000000002        2
2.00000000014        2.00000000015        3
6.00000000051        6.00000000056        4
24.0000000017        24.0000000019        5
120.000000010        120.000000002        6
720.000000110        720.000000111        7
5040.00000097        5040.00000083        8
40320.0000039        40320.0000025        9
362880.000055        362879.999993       10
3628800.00049        3628800.00034       11
39916800.0045        39916800.0035       12
479001600.110        479001600.007       13
6227020799.76        6227020800.96       14
0.12000000000        0.97000000000


Thanks Olivier and Gerson for your input. As Voltaire once said, "The better is the enemy of the good!"

I tried to do a curve fit between 1/gamma(x) and a tenth order polynomial. I also tired a Pade approximation using fifth order polynomials in the numerator and denominator. Neither attempts yielded good results.

Namir

The Nemes approximations that Viktor Toth mentions are good ones, but not as good as the Spouge and Lanczos approximations. The Nemes approximation DO come third in my little study!

Hello Namir,

20 years ago I programmed the gamma-function in "turbo pascal 6" with assembler routines coded for the 387 coprocessor. The Stirling-formula was used (for arguments > 10),

for smaller arguments the recursion (gamma(x) = gamma(x+1)/x)).

For negative Arguments the equation:

gamma(x) = pi/(sin(pi*x)*gamma(1-x)).

The function was only usefull for real numbers, and it was a luck, that I didn't had to earn my money with programming....

Greetings peacecalc

Interesting that you mentioned Turbo Pascal. I remember implementing the gamma function in Turbo Pascal in the late eighties. I used the series expansion that employs 26 constants too implement a polynomial approximation for 1/Gamma(x) for 1<=x<=2. I used recursion for arguments that were greater than 2.

I made a living then by writing books about programming in Turbo Pascal, and then switched to Visual Basic and Visual C++.

Namir