The Museum of HP Calculators

HP Forum Archive 19

[ Return to Index | Top of Index ]

gamma approximation in 48/49/50?
Message #1 Posted by Eric Smith on 24 Sept 2009, 3:04 a.m.

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?
Message #2 Posted by Ángel Martin on 24 Sept 2009, 10:57 a.m.,
in response to message #1 by Eric Smith

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?
Message #3 Posted by Eric Smith on 24 Sept 2009, 3:50 p.m.,
in response to message #2 by Ángel Martin

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?
Message #4 Posted by Eric Smith on 24 Sept 2009, 6:54 p.m.,
in response to message #2 by Ángel Martin

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?
Message #5 Posted by Ángel Martin on 25 Sept 2009, 1:23 a.m.,
in response to message #4 by Eric Smith

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?
Message #6 Posted by hugh steers on 24 Sept 2009, 11:56 a.m.,
in response to message #1 by Eric Smith

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?
Message #7 Posted by Frank Rottgardt on 24 Sept 2009, 12:05 p.m.,
in response to message #6 by hugh steers

My 48GX Singapore 3438S02321 give me the same result as you 50g.

                  
Re: gamma approximation in 48/49/50?
Message #8 Posted by hugh steers on 24 Sept 2009, 12:24 p.m.,
in response to message #7 by Frank Rottgardt

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?
Message #9 Posted by Eric Smith on 24 Sept 2009, 3:48 p.m.,
in response to message #8 by hugh steers

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?
Message #10 Posted by Tim Wessman on 24 Sept 2009, 6:26 p.m.,
in response to message #9 by Eric Smith

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?
Message #11 Posted by Eric Smith on 24 Sept 2009, 6:49 p.m.,
in response to message #10 by Tim Wessman

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?
Message #12 Posted by hugh steers on 24 Sept 2009, 9:37 p.m.,
in response to message #11 by Eric Smith

Something subtle is missing. these constants, they don't much change magnitude are there some exponents missing?

                                                
Re: gamma approximation in 48/49/50?
Message #13 Posted by Eric Smith on 24 Sept 2009, 11:53 p.m.,
in response to message #12 by hugh steers

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?
Message #14 Posted by Tim Wessman on 25 Sept 2009, 12:11 a.m.,
in response to message #13 by Eric Smith

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?
Message #15 Posted by Pal G. on 25 Sept 2009, 12:33 a.m.,
in response to message #13 by Eric Smith

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?
Message #16 Posted by Eric Smith on 25 Sept 2009, 4:11 a.m.,
in response to message #15 by Pal G.

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?
Message #17 Posted by Pal G. on 25 Sept 2009, 11:41 a.m.,
in response to message #16 by Eric Smith

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?
Message #18 Posted by Steven Thomas Smith on 25 Sept 2009, 11:46 a.m.,
in response to message #1 by Eric Smith

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?
Message #19 Posted by Eric Smith on 25 Sept 2009, 5:25 p.m.,
in response to message #18 by Steven Thomas Smith

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?
Message #20 Posted by hugh steers on 25 Sept 2009, 8:06 p.m.,
in response to message #19 by Eric Smith

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.

[ Return to Index | Top of Index ]

Go back to the main exhibit hall