Increasing the accuracy of Gamma for double float

Dieter Kaiser wrote:
> I have done some further work to improve the precision of the Gamma function for
> double float calculation.
> 
> This is the actual result for gamma-lanczos:
> ...................................................
> Test of precision for gamma-lanczos
> gamma((float(z)) against gamma(bfloat(z)), fpprec:32
> z: 0.5 thru 150.0 step 0.25 (values about 170.0 give an overflow)

I modified gammafloat in csimp2.lisp to:

(defun gammafloat (a)
  (slatec::dgamma (float a)))

A quick test says

meanerr = 2.35b-14
maxerr  = 1.27b-13 at 142.75

So gamma-lanczos is actually slightly more accurate than the SLATEC
routine.  I'm a bit surprised.  I would have expected the SLATEC routine
to do better.

Also, the slatec routine does pretty well until x = 10 at which point
the errors start to grow.  x = 10 happens to be the point where the
slatec routine switches from a series to something like Stirling's
formula.  The code claims the approximation there is accurate to
1.28d-31 (which we can't actually achieve since we only have
double-float precision).  Perhaps the computation of (x-1/2*log(x)-x
causes the loss of precision we are seeing.

Ray