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.