What is the Gamma approximation you use? Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-02-2013, 04:46 PM 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 Walter B Unregistered Posts: 4,587 Threads: 105 Joined: Jul 2005 08-02-2013, 05:19 PM Quote: What does the WP34S use for the Gamma function? Just look at it - it's open source. d;-) Barry Mead Unregistered Posts: 65 Threads: 14 Joined: Jun 2013 08-02-2013, 07:02 PM 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. Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-02-2013, 07:39 PM 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 Gerson W. Barbosa Unregistered Posts: 2,761 Threads: 100 Joined: Jul 2005 08-02-2013, 07:57 PM 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. Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 08-02-2013, 08:37 PM 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 Ángel Martin Unregistered Posts: 1,253 Threads: 117 Joined: Nov 2005 08-03-2013, 02:22 AM 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. Kimberly Thompson Unregistered Posts: 97 Threads: 1 Joined: Jun 2013 08-03-2013, 09:41 AM 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 Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-03-2013, 10:41 AM 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. Kimberly Thompson Unregistered Posts: 97 Threads: 1 Joined: Jun 2013 08-03-2013, 11:17 AM THANKS! ps: I am indebted to you for ALL your marvelous postings (here & your web page)! Many thanks, again! SlideRule Gerson W. Barbosa Unregistered Posts: 2,761 Threads: 100 Joined: Jul 2005 08-03-2013, 12:05 PM 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. Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-03-2013, 02:39 PM 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. Gerson W. Barbosa Unregistered Posts: 2,761 Threads: 100 Joined: Jul 2005 08-03-2013, 05:20 PM 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  Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-03-2013, 06:39 PM 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 Les Koller Unregistered Posts: 253 Threads: 20 Joined: Jun 2012 08-04-2013, 12:58 AM I took a look at your page also. Bookmarked it in fact. I even took a couple things. Fantastic site, thanks so much. Marcus von Cube, Germany Unregistered Posts: 3,283 Threads: 104 Joined: Jul 2005 08-04-2013, 06:47 AM 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"},  Olivier De Smet Unregistered Posts: 275 Threads: 41 Joined: Mar 2010 08-04-2013, 10:50 AM 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. Gerson W. Barbosa Unregistered Posts: 2,761 Threads: 100 Joined: Jul 2005 08-04-2013, 01:40 PM 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  Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-04-2013, 06:41 PM 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 Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-04-2013, 06:59 PM 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! peacecalc Unregistered Posts: 97 Threads: 9 Joined: Nov 2011 08-05-2013, 03:03 AM 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 Namir Unregistered Posts: 2,247 Threads: 200 Joined: Jun 2005 08-05-2013, 07:14 AM 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 « Next Oldest | Next Newest »

 Possibly Related Threads… Thread Author Replies Views Last Post Euler's Constant Approximation Namir 0 926 10-19-2013, 11:32 AM Last Post: Namir HP15c continued fraction for Ln(Gamma) Tom Grydeland 0 1,108 09-30-2013, 05:48 AM Last Post: Tom Grydeland HP Prime emulator: Gamma function Stephan Matthys 28 7,935 08-21-2013, 04:52 PM Last Post: Namir SandMath routine of the week: Inverse Gamma Function Ángel Martin 39 9,416 03-24-2013, 08:19 AM Last Post: peacecalc Question about trig functions approximation Namir 8 2,507 01-10-2013, 06:49 AM Last Post: Valentin Albillo [wp34s] Incomplete Gamma on the wp34s Les Wright 18 5,117 12-06-2011, 11:07 AM Last Post: Namir gamma approximation in 48/49/50? Eric Smith 19 4,640 09-25-2009, 08:06 PM Last Post: hugh steers Complex Gamma Function... revisited Ángel Martin 13 3,492 09-22-2009, 08:51 AM Last Post: Ángel Martin (O.T.) Happy Pi Approximation Day! 'e*XROOT(12,e^(-3*4)+5.6789)' (N.T.) Gerson W. Barbosa 6 1,986 07-23-2009, 03:18 AM Last Post: Thomas Radtke 12C+ Gamma Gerson W. Barbosa 8 2,415 06-30-2009, 07:51 AM Last Post: Gerson W. Barbosa

Forum Jump: