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.