# HP Forums

Full Version: gamma approximation in 48/49/50?
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

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

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.

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?

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

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?

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.

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.

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.

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.

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.)

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

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.

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

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?

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?

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.

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 ;)

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.

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?

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.