▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
Quote:
What does the WP34S use for the Gamma function?
Just look at it - it's open source.
d;-)
▼
Posts: 65
Threads: 14
Joined: Jun 2013
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.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 1,253
Threads: 117
Joined: Nov 2005
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
▼
Posts: 3,283
Threads: 104
Joined: Jul 2005
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"},
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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!
Posts: 97
Threads: 1
Joined: Jun 2013
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
▼
Posts: 97
Threads: 1
Joined: Jun 2013
THANKS!
ps: I am indebted to you for ALL your marvelous postings (here & your web page)!
Many thanks, again!
SlideRule
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
▼
Posts: 253
Threads: 20
Joined: Jun 2012
I took a look at your page also. Bookmarked it in fact. I even took a couple things. Fantastic site, thanks so much.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
▼
Posts: 275
Threads: 41
Joined: Mar 2010
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
Posts: 97
Threads: 9
Joined: Nov 2011
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
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
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
|
