Increasing the accuracy of Gamma for double float

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)

meanerr = 2.01b-14
maxerr  = 1.15b-13 for value 142.
...................................................

Now we are calling bffac and convert the bigfloat result to a double float:

($float (mfuncall '$bffac ($bfloat (+ j -1)) $gamma_float_fpprec))

These are the results for different values of the precision:
...................................................
Test of precision with gamma_float_fpprec : 16

meanerr = 4.09b-15
maxerr  = 4.34b-13 for value 126.
...................................................
Test of precision with gamma_float_fpprec : 17

meanerr = 5.97b-16
maxerr  = 4.34b-13 for value 126.
...................................................
Test of precision with gamma_float_fpprec : 18

meanerr = 8.28b-17
maxerr  = 4.72b-16 for value 138.
...................................................
Test of precision with gamma_float_fpprec : 19

meanerr = 3.8b-17
maxerr  = 1.06b-16 for value 145.
...................................................
Test of precision with gamma_float_fpprec : 20

meanerr = 3.78b-17
maxerr  = 1.06b-16 for value 145.
...................................................
Test of precision with gamma_float_fpprec : 21

meanerr = 3.78b-17
maxerr  = 1.06b-16 for value 145.
...................................................
Test of precision with gamma_float_fpprec : 22

meanerr = 3.78b-17
maxerr  = 1.06b-16 for value 145.
...................................................

With 20 digits in the call to bffac we get a max error of about 1.0e-16 and a
mean error of about 3.8e-17 for values between 0.5 and 150.0. The precision do
not change anymore with increasing digits. Thus we have an improvement of about
3 digits for double float precision.

But using bffac to calculate gamma in double float precision we lost about a
factor 50 in speed.

I have done some further work and tested two different sets of coefficients for
gamma lanczos which claim to be accurate within 16 and 32 digits. Furthermore I
have tried to improve bffac. To do this, I have extended the routine $bern with
a function $bern_numerical which return the desired numerical approximation of
the bernoulli number and stores the results for a next call in an array. But we
do not gain speed in the bigfloat calculation of bffac.

The best result with a gamma-lanczos algorithm I have obtained with an 
implementation in 32 digits bigfloat precision. We get the desired precision for
double floats but lost again a factor of about 40 in speed.

This would be the code gamma-lanczos32 in bigfloat arithmetic:

(defun gamma-lanczos32 (z)
  (let (($fpprec 32)
        (c (make-array 22 :initial-contents
            '(((BIGFLOAT SIMP 109) 564220969521603723634314552515685 -32)
             ((BIGFLOAT SIMP 109) 497590237669501352603200443655143 1)
             ((BIGFLOAT SIMP 109) -472185635031721359996050261461607 4)
             ((BIGFLOAT SIMP 109) 406192754490738304492679256416322 6)
             ((BIGFLOAT SIMP 109) -419175866985132316842978815219736 7)
             ((BIGFLOAT SIMP 109) 578782181229217435525800744870256 7)
             ((BIGFLOAT SIMP 109) -564644313579207123184023795214178 7)
             ((BIGFLOAT SIMP 109) 400767491587913288793896252396625 7)
             ((BIGFLOAT SIMP 109) -420002345461945986241865118378769 6)
             ((BIGFLOAT SIMP 109) 326440153160455391086338992610886 5)
             ((BIGFLOAT SIMP 109) -375012098974408151314567252241494 3)
             ((BIGFLOAT SIMP 109) 629766444628063115246938177974192 0)
             ((BIGFLOAT SIMP 109) -379228030153090439281885993255360 -2)
             ((BIGFLOAT SIMP 109) 636796808749833654894427708162278 -6)
             ((BIGFLOAT SIMP 109) -358022988271483531754979749441415 -9)
             ((BIGFLOAT SIMP 109) 509478510919234905167759855642288 -14)
             ((BIGFLOAT SIMP 109) -423235757251797647448160467134330 -19)
             ((BIGFLOAT SIMP 109) 364804891535155920838404979294590 -25)
             ((BIGFLOAT SIMP 109) -544554409365465828007543356200845 -33)
             ((BIGFLOAT SIMP 109) 523131301064851815510611164937661 -42)
             ((BIGFLOAT SIMP 109) -377116395533824638758034339754730 -53)
             ((BIGFLOAT SIMP 109) 512721769168209434341912770400693 -69)))))
    (let* ((z ($bfloat (- z 1)))
           (zh (add z `((rat) 1 2)))
           (zgh (div (add zh ($bfloat 22.618910)) ($bfloat '$%e)))
           (zp (power zgh zh))
           (ss 
             (do ((sum ($bfloat 0.0))
                  (pp (1- (length c)) (1- pp)))
                 ((< pp 1)
                   sum)
                 (setq sum (add sum (div (aref c pp) (add z pp)))))))
      ($float 
        (mul 2
          (power (div ($bfloat '$%e) ($bfloat '$%pi)) `((rat) 1 2))
          zp
          (add ss (aref c 0)))))))

Thus, if we want to have more accurate double float results for the gamma
function it seems to be the best to use bffac.

One problem is that with GCL we do not get an overflow error when converting a
bigfloat number which does not fit into a double float. We get wrong results.
See bug report SF[2180110) on Sourceforge.net. But this is an problem of the
compiler and not of the algorithm.

So, what do you think? Should we improve the accuracy and do a call to bffac
with an preciscion of fpprec:20? Is it worth to lost a speed of 50 to gain the
missing 3 digits in double float precision? What about the bug we have with GCL?

Dieter Kaiser