Lanczos Gamma formula (Re: AW: AW: Increasing the accuracy of Gamma for double float)

Just for fun, I implemented Lanczos formula for the factorial function,
as used in our gamma-lanczos routine and as described in
http://my.fit.edu/~gabdo/gamma.txt.

I've attached a very simple and straightforward implementation described
in the link.  Since maxima is a symbolic math system, all of the
computations are symbolic, of course. :-)

The coefficients used in our implementation are given by
lanczos_coeff(607/128,15).  Actually, the coefficients are scaled, so to
get the actual coefficients, we need to compute
lanczos_coeff(607/128,15)*exp(607/128)/sqrt(2*%pi).

To compute the cofficients accurately, we need to use bfloat arithmetic.
 I set fpprec = 100, and convert the result back to float.  Thus
float(bfloat(lanczos_coeff(607/128,15)*exp(607/128)/sqrt(2*%pi))) gives

                          [    1.000000000000001    ]
                          [                         ]
                          [    57.15623566586287    ]
                          [                         ]
                          [   - 59.5979603554632    ]
                          [                         ]
                          [    14.13609797406266    ]
                          [                         ]
                          [   - .4919137997996221   ]
                          [                         ]
                          [  3.378175738533551e-5   ]
                          [                         ]
                          [  4.822676972285958e-5   ]
                          [                         ]
                          [ - 1.0731596448288737e-4 ]
(%o303)                   [                         ]
                          [  1.9020507251734867e-4  ]
                          [                         ]
                          [ - 2.910571832253896e-4  ]
                          [                         ]
                          [  3.6107115856194106e-4  ]
                          [                         ]
                          [ - 3.4385753209517444e-4 ]
                          [                         ]
                          [  2.387331021497593e-4   ]
                          [                         ]
                          [ - 1.1299295736990386e-4 ]
                          [                         ]
                          [  3.2452750586817404e-5  ]
                          [                         ]
                          [ - 4.255725190331955e-6  ]

(Using more precision doesn't change these results, so I think they're
correct).  Note that the first 5 coefficents match the code in csimp2,
but the rest are quite different.  Random spot checks show that the
computed values are correct.

Not sure what the conclusion is, but this was kind of fun.  I hope I
didn't make a mistake somewhere.

Ray