What is the Gamma approximation you use?

The Museum of HP Calculators

HP Forum Archive 21

What is the Gamma approximation you use?
Message #1 Posted by Namir on 2 Aug 2013, 4:46 p.m.

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

Re: What is the Gamma approximation you use?
Message #2 Posted by Walter B on 2 Aug 2013, 5:19 p.m.,
in response to message #1 by Namir

Quote:
What does the WP34S use for the Gamma function?
Just look at it - it's open source.
d;-)

Re: What is the Gamma approximation you use?
Message #3 Posted by Barry Mead on 2 Aug 2013, 7:02 p.m.,
in response to message #2 by Walter B

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.

Re: What is the Gamma approximation you use?
Message #4 Posted by Namir on 2 Aug 2013, 7:39 p.m.,
in response to message #3 by Barry Mead

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

Re: What is the Gamma approximation you use?
Message #5 Posted by Gerson W. Barbosa on 2 Aug 2013, 7:57 p.m.,
in response to message #4 by Namir

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.

Re: What is the Gamma approximation you use?
Message #6 Posted by Ángel Martin on 3 Aug 2013, 2:22 a.m.,
in response to message #2 by Walter B

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.


Re: What is the Gamma approximation you use?
Message #7 Posted by Paul Dale on 2 Aug 2013, 8:37 p.m.,
in response to message #1 by Namir

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

Re: What is the Gamma approximation you use?
Message #8 Posted by Marcus von Cube, Germany on 4 Aug 2013, 6:47 a.m.,
in response to message #7 by Paul Dale

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"},


Re: What is the Gamma approximation you use?
Message #9 Posted by Namir on 4 Aug 2013, 6:59 p.m.,
in response to message #8 by Marcus von Cube, Germany

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!

Re: What is the Gamma approximation you use?
Message #10 Posted by Kimberly Thompson on 3 Aug 2013, 9:41 a.m.,
in response to message #1 by Namir

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

Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #11 Posted by Namir on 3 Aug 2013, 10:41 a.m.,
in response to message #1 by Namir

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.

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #12 Posted by Kimberly Thompson on 3 Aug 2013, 11:17 a.m.,
in response to message #11 by Namir

THANKS!
ps: I am indebted to you for ALL your marvelous postings (here & your web page)!
Many thanks, again!
SlideRule

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #13 Posted by Namir on 3 Aug 2013, 2:39 p.m.,
in response to message #12 by Kimberly Thompson

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.

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #14 Posted by Les Koller on 4 Aug 2013, 12:58 a.m.,
in response to message #13 by Namir

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

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #15 Posted by Gerson W. Barbosa on 3 Aug 2013, 12:05 p.m.,
in response to message #11 by Namir

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.

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #16 Posted by Gerson W. Barbosa on 3 Aug 2013, 5:20 p.m.,
in response to message #11 by Namir

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


Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #17 Posted by Namir on 3 Aug 2013, 6:39 p.m.,
in response to message #16 by Gerson W. Barbosa

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

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #18 Posted by Olivier De Smet on 4 Aug 2013, 10:50 a.m.,
in response to message #17 by Namir

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.

Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #19 Posted by Gerson W. Barbosa on 4 Aug 2013, 1:40 p.m.,
in response to message #18 by Olivier De Smet

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


Re: Gamma Approxiamtion for HP-71B, HP41C, and HP67
Message #20 Posted by Namir on 4 Aug 2013, 6:41 p.m.,
in response to message #17 by Namir

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

Re: What is the Gamma approximation you use?
Message #21 Posted by peacecalc on 5 Aug 2013, 3:03 a.m.,
in response to message #1 by Namir

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

Re: What is the Gamma approximation you use?
Message #22 Posted by Namir on 5 Aug 2013, 7:14 a.m.,
in response to message #21 by 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

[ Return to Index | Top of Index ]

Go back to the main exhibit hall