gamma approximation in 48/49/50? « Next Oldest | Next Newest »

 ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-24-2009, 03:04 AM I've been poking around in the 48GX, and I'm curious about the approximation of the gamma function it uses. For real non-integer arguments greater than 6, it uses an approximation which I've rewritten in (ugly) C code: ```double gamma (double x) { const double k0 = -49.5920119017593; const double k1 = 3.03479606073615; const double k2 = -7.65594818140208; const double k3 = 0.933584905660377; const double k4 = -0.121142857142857; const double k5 = 0.4; const double k6 = 0.918938533204673; double tx2 = 12.0 * x * x; // 12 x squared double ln_gamma = x/(k4/(k2/((k0/(tx2+60))+tx2+k1)+tx2+k3)+tx2+k5)-x+k6+(x-0.5)*log(x); return exp (ln_gamma); } ``` (Note that in C the log(x) function is the natural log, and exp(x) returns e to the x.) This doesn't appear to exactly match any of the approximation formulas I've seen elsewhere. Is this a rearrangement of one of the well-known approximation formulas? [edited 24-SEP-2009 15:45 PDT to correct transcription of k1] Edited: 24 Sept 2009, 6:45 p.m. after one or more responses were posted ▼ Ángel Martin Unregistered Posts: 1,253 Threads: 117 Joined: Nov 2005 09-24-2009, 10:57 AM if memory serves me (which is doubtful these days) that looks rather similar to the method used in the HP-29 High-Level Math application booklet. Will check when I'm back home... maybe I just don;t recall it right. Could it be HP's used the same thing ever since the 29C days? Edited: 24 Sept 2009, 10:58 a.m. ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-24-2009, 03:50 PM It's possible. I haven't yet figured out how the real gamma function in the 34C works, and I haven't even started looking at the 15C gamma function code yet. Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-24-2009, 06:54 PM I found the 19C/29C Mathematics Solutions book you're referencing on the MoHPC DVD, and the section on the Gamma function references an HP-25 program in 65 Notes v3n10, which I found on Jake's PPC CD-ROM. It does look similar, though the 48GX code uses more terms, presumably to get more digits of accuracy. It's unfortunate that the authors of that program don't give a derivation or reference. (Then again, maybe it is obvious to someone with a better mathematics education than I have.) ▼ Ángel Martin Unregistered Posts: 1,253 Threads: 117 Joined: Nov 2005 09-25-2009, 01:23 AM Interestingly HP seemed to have abandoned this approximation in other solutions books published afterwards (see the 67 and 41 HL math for instance) just to return to it (allegedly) on the saturn machines... vagaries of corporate design teams perhaps? hugh steers Unregistered Posts: 536 Threads: 56 Joined: Jul 2005 09-24-2009, 11:56 AM using this approx on 7.22! gives, 7877.71704167 but my hp50 gives 7877.71717234 which is correct. do you get the former on your 48gx? ▼ Frank Rottgardt Unregistered Posts: 115 Threads: 16 Joined: Jul 2007 09-24-2009, 12:05 PM My 48GX Singapore 3438S02321 give me the same result as you 50g. ▼ hugh steers Unregistered Posts: 536 Threads: 56 Joined: Jul 2005 09-24-2009, 12:24 PM So Eric, is this really the formula used, or could it be to do with the number system in the 48 vs the "C" code? ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-24-2009, 03:48 PM It's possible that I didn't transcribe one of the constants from the 48GX correctly, as I was doing it by hand, and I haven't got time to double-check it. (I forgot to write down the address of the code in the 48GX ROM, so I'd have to find it again.) Anyhow, unless I screwed it up, my C code is essentially the logic used in the 48GX factorial function. They use other code for dealing with integers, and with values less than 6. I assume that they're using the standard reflection formula to handle negative numbers, but I didn't study the code for it. Presumably they chose 6 as the lower bound of the domain because the approximation is less accurate for smaller values. My point wasn't the C code would be useful for anything. I just wondered if anyone was familiar with this approximation. ▼ Tim Wessman Unregistered Posts: 1,278 Threads: 44 Joined: Jul 2007 09-24-2009, 06:26 PM I just looked in the source and there wasn't any clue there at all. Here are the constants pulled out of there to check your transcription. NIBHEX 395710911029594 NIBHEX 516370606974303 NIBHEX 802041818495567 NIBHEX 773066509485339 NIBHEX 758241758241121 NIBHEX 000000000000004 NIBHEX 376402335839819 For those not familiar with saturn encoding, like all HP stuff it is sdrawkcab. . . :-) TW Edited: 24 Sept 2009, 6:29 p.m. ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-24-2009, 06:49 PM Thanks, Tim! I did drop a digit in the constant I designated k1, so I've updated the code. On my x86-64 Linux system, it still doesn't perform as well as the actual HP-48GX. IEEE double precision should have just over 15 digits of precision, so I don't think that's the problem. Perhaps the natural log and exponential functions aren't as accurate. ▼ hugh steers Unregistered Posts: 536 Threads: 56 Joined: Jul 2005 09-24-2009, 09:37 PM Something subtle is missing. these constants, they don't much change magnitude are there some exponents missing? ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-24-2009, 11:53 PM I think I have the exponents correct, but it would be helpful for someone else to study the actual HP-48 (or -49 or 50) code to see if they spot any obvious mistakes in my C rewrite. ▼ Tim Wessman Unregistered Posts: 1,278 Threads: 44 Joined: Jul 2007 09-25-2009, 12:11 AM I looked at it quickly when I grabbed those constants. I didn't see anything that looked out of place. I didn't look long though, so I easily could have missed something. TW Pal G. Unregistered Posts: 260 Threads: 21 Joined: Nov 2006 09-25-2009, 12:33 AM I do believe Eric's constant k1 and Tim's hexbin line 2 do not match the other constant's regarding their backwardness. ie.: ```constant k0 = ABCDEFG : nibhex GFEDCBA constant k1 = GFEDCBA : nibhex ABCDEFG constant k2 = ABCDEFG : nibhex GFEDCBA constant k3 = ABCDEFG : nibhex GFEDCBA constant k4 = ABCDEFG : nibhex GFEDCBA constant k5 = ABCDEFG : nibhex GFEDCBA constant k6 = ABCDEFG : nibhex GFEDCBA ``` Does that matter? ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-25-2009, 04:11 AM I appreciate your taking it look at it, but either I don't understand what you are suggesting, or you are mistaken. I double-checked and am fairly certain that all of my k[i] constants are in the reverse digit order of the NIBHEX strings Tim posted. I did have a dropped digit in the middle of k1, which I've since corrected, but that didn't improve the accuracy to match the actual 48GX. ▼ Pal G. Unregistered Posts: 260 Threads: 21 Joined: Nov 2006 09-25-2009, 11:41 AM I apologize Eric. My brain seems to have failed last night. I'm looking at it now and I can't figure out where my head went wrong. I thought k1 was mirrored wrong compared to the other entries. I will move to another thread and MYOB ;) Steven Thomas Smith Unregistered Posts: 2 Threads: 0 Joined: Sep 2009 09-25-2009, 11:46 AM That's simply the asymptotic expansion of log(Gamma(x)), rearranged as a continued fraction. For example, you can recognize the k6 constant as log(2*pi)/2 = 0.91893853320467267. The coefficients are all Bernoulli numbers, which are spelled out explicitly in the asymptotic expansion for Gamma(z). A better implementation would use Lanczos's approximation. ▼ Eric Smith Unregistered Posts: 2,309 Threads: 116 Joined: Jun 2005 09-25-2009, 05:25 PM Thanks! I knew someone here would recognize it! Quote: A better implementation would use Lanczos's approximation. Is that still a better method than the Nemes approximation? Or the Stieltjes approximation, which along with many others is described here? ▼ hugh steers Unregistered Posts: 536 Threads: 56 Joined: Jul 2005 09-25-2009, 08:06 PM maybe it really does use this formula and there is no mistake here. except there is an outer "driver" routine that conditions the input. say we want 12 accurate digits. I have the (alleged) 48 logic as: ```double _gamma48gx(double x) { const double k0 = -49.5920119017593; const double k1 = 3.03479606073615; const double k2 = -7.65594818140208; const double k3 = 0.933584905660377; const double k4 = -0.121142857142857; const double k5 = 0.4; const double k6 = 0.918938533204673; double tx2 = 12.0 * x * x; // 12 x squared double ln_gamma = x/(k4/(k2/((k0/(tx2+60))+tx2+k1)+tx2+k3)+tx2+k5)-x+k6+(x-0.5)*log(x); return exp (ln_gamma); } ``` this barely gets about 8 digits for numbers starting at 6. But, if you bump all inputs up to around 100, things are ok. eg ```double gamma48gx(double x) { double g; double p = 1.0; while (x < 100) { p *= x; x = x + 1; } return _gamma48gx(x)/p; } ``` ok, its not very good because you have to perform all those multiplies, but it does get the right answer. so i had a look at what also might be done with the same amount of work (ie 4 divides, some adds plus log/exp), i am finding the Stieltjes forms quite good. for example, i get this: ```double gammaStieltjes4(double x) { // get 12 digits x > 8 const double a = 1.0/12; const double b = 1.0/30; const double c = 53.0/210; const double d = 129.0/250; // adapted // would be a constant double e = sqrt(2*PI); double s; s = a/(x + b/(x + c/(x + d/x))); return e*exp((x-.5)*log(x) - x + s); } ``` this starts getting 12 correct digits somewhere between 7 and 8 and would be easy to bump up in the same was suggested above, but with only a few extra multiples (maybe just one). Either that or use the next best version of the approx (ie one more term). in this example, i got surprisingly accurate results from this 4 divide logic. actually this was a mistake, i had a typo in the last coefficient that turned out to make it *more* accurate! After some experimentation, it turns out to be a benefit to adjust the last coefficient term to take account of the fact that you are stopping there (can compute the adjustment from the next term ignored). I did this and arrived at my "adapted" term here. i still don't really know what the 48 actually does.

 Possibly Related Threads… Thread Author Replies Views Last Post HP Prime... NOT meant to replace HP48,49,50 ? Chris Pem10 21 5,642 11-18-2013, 03:30 PM Last Post: Chris Smith Euler's Constant Approximation Namir 0 922 10-19-2013, 11:32 AM Last Post: Namir HP15c continued fraction for Ln(Gamma) Tom Grydeland 0 1,099 09-30-2013, 05:48 AM Last Post: Tom Grydeland ANN: A MacOSX Kermit program for HP48/49/50 calculators Paul Onions 2 1,603 09-08-2013, 04:55 AM Last Post: Paul Onions HP Prime emulator: Gamma function Stephan Matthys 28 7,881 08-21-2013, 04:52 PM Last Post: Namir What is the Gamma approximation you use? Namir 21 5,250 08-05-2013, 07:14 AM Last Post: Namir HP-50 on Raspberry Pi? The HP-67 come true? Matti Övermark 10 3,334 07-30-2013, 09:40 PM Last Post: Matti Övermark Nibble reverse (HP-48,49,50g) Gerson W. Barbosa 44 12,075 07-28-2013, 10:10 PM Last Post: Gerson W. Barbosa (OT) Pandigital expression (HP-48,49,50g) Gerson W. Barbosa 19 5,475 07-19-2013, 05:59 AM Last Post: Gilles Carpentier Precision DC-50 Mic 3 1,519 05-26-2013, 01:27 PM Last Post: DavidM

Forum Jump: