 What is the Gamma approximation you use? - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: What is the Gamma approximation you use? (/thread-247662.html) What is the Gamma approximation you use? - Namir - 08-02-2013 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? - Walter B - 08-02-2013 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? - Barry Mead - 08-02-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. Re: What is the Gamma approximation you use? - Namir - 08-02-2013 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? - Gerson W. Barbosa - 08-02-2013 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? - Paul Dale - 08-02-2013 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? - Ángel Martin - 08-03-2013 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? - Kimberly Thompson - 08-03-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 Gamma Approxiamtion for HP-71B, HP41C, and HP67 - Namir - 08-03-2013 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 - Kimberly Thompson - 08-03-2013 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 - Gerson W. Barbosa - 08-03-2013 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 - Namir - 08-03-2013 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 - Gerson W. Barbosa - 08-03-2013 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 - Namir - 08-03-2013 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 - Les Koller - 08-04-2013 I took a look at your page also. Bookmarked it in fact. I even took a couple things. Fantastic site, thanks so much. Re: What is the Gamma approximation you use? - Marcus von Cube, Germany - 08-04-2013 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: Gamma Approxiamtion for HP-71B, HP41C, and HP67 - Olivier De Smet - 08-04-2013 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 - Gerson W. Barbosa - 08-04-2013 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 - Namir - 08-04-2013 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? - Namir - 08-04-2013 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? - peacecalc - 08-05-2013 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? - Namir - 08-05-2013 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