gamma approximation in 48/49/50? - 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: gamma approximation in 48/49/50? (/thread-156775.html) gamma approximation in 48/49/50? - Eric Smith - 09-24-2009 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 Re: gamma approximation in 48/49/50? - Ángel Martin - 09-24-2009 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. Re: gamma approximation in 48/49/50? - hugh steers - 09-24-2009 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? Re: gamma approximation in 48/49/50? - Frank Rottgardt - 09-24-2009 My 48GX Singapore 3438S02321 give me the same result as you 50g. Re: gamma approximation in 48/49/50? - hugh steers - 09-24-2009 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? Re: gamma approximation in 48/49/50? - Eric Smith - 09-24-2009 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. Re: gamma approximation in 48/49/50? - Eric Smith - 09-24-2009 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. Re: gamma approximation in 48/49/50? - Tim Wessman - 09-24-2009 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. Re: gamma approximation in 48/49/50? - Eric Smith - 09-24-2009 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. Re: gamma approximation in 48/49/50? - Eric Smith - 09-24-2009 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.) Re: gamma approximation in 48/49/50? - hugh steers - 09-24-2009 Something subtle is missing. these constants, they don't much change magnitude are there some exponents missing? Re: gamma approximation in 48/49/50? - Eric Smith - 09-24-2009 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. Re: gamma approximation in 48/49/50? - Tim Wessman - 09-25-2009 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 Re: gamma approximation in 48/49/50? - Pal G. - 09-25-2009 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? Re: gamma approximation in 48/49/50? - Ángel Martin - 09-25-2009 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? Re: gamma approximation in 48/49/50? - Eric Smith - 09-25-2009 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. Re: gamma approximation in 48/49/50? - Pal G. - 09-25-2009 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 ;) Re: gamma approximation in 48/49/50? - Steven Thomas Smith - 09-25-2009 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. Re: gamma approximation in 48/49/50? - Eric Smith - 09-25-2009 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? Re: gamma approximation in 48/49/50? - hugh steers - 09-25-2009 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.