More data concerning the precision of bigfloat calculation of gamma

I have caculated some more data to show the precision of the calculation of the
Gamma function. 

The function is tested for 150 values between 1/2 through 75 in steps of 1/2.
For the different precisions the mean error, the minimum and the maximum error
are tabulated. 

The tests are done up to a precision of fpprec : 4096 and with GCL 2.6.8.

In addition the time needed on my system is given. This is not the time for
calucalation 150 values of the Gamma function, but the time needed for the whole
test function.

I think it is nice to see, that the values of minerr and maxerr are close
together and we get any desired accuracy with increasing fpprec. Perhaps we can
even increase the accuracy a bit, so all digits are significant. As mentioned
earlier on this mailing list for double float precision we have to get a new set
of coefficients for the algorithm of gamma-lanczos to get more significant
digits.

By the way I have detected a bug in the calculation of the Gamma function in
bigfloat precision: When increasing fpprec the first call to the Gamma function
do not give the full precision given by fpprec. There seems to be a bug in bffac
or in $bern which causes this error. I will give a bug report. 

...................................................
Test of precision for double float

meanerr =  9.28E-15
minerr  =   0.00E+0 for value 45/2
maxerr  =  4.96E-14 for value 131/2

Evaluation took 0.0400 seconds (0.0400 elapsed)

...................................................
Test of precision with fpprec : 16

meanerr = 2.05b-15
minerr  = 1.59b-17 for value 30
maxerr  = 9.16b-15 for value 129/2

Evaluation took 0.9700 seconds (0.9700 elapsed)

...................................................
Test of precision with fpprec : 32

meanerr = 3.37b-31
minerr  = 4.45b-33 for value 10
maxerr  = 1.39b-30 for value 70

Evaluation took 1.5300 seconds (1.5300 elapsed)

...................................................
Test of precision with fpprec : 64

meanerr = 1.45b-62
minerr  = 6.56b-65 for value 35
maxerr  = 4.85b-62 for value 75

Evaluation took 2.4700 seconds (2.4700 elapsed)

...................................................
Test of precision with fpprec : 128

meanerr = 2.95b-126
minerr  = 2.56b-129 for value 137/2
maxerr  = 7.67b-126 for value 75/2

Evaluation took 4.4700 seconds (4.4700 elapsed)

...................................................
Test of precision with fpprec : 256

meanerr = 9.67b-254
minerr  = 7.21b-255 for value 15
maxerr  = 1.92b-253 for value 147/2

Evaluation took 9.7900 seconds (9.7900 elapsed)

...................................................
Test of precision with fpprec : 512

meanerr = 9.69b-510
minerr  = 7.59b-512 for value 11/2
maxerr  = 2.53b-509 for value 73

Evaluation took 21.6700 seconds (21.6700 elapsed)

...................................................
Test of precision with fpprec : 1024

meanerr = 2.64b-1020
minerr  = 2.01b-1020 for value 11
maxerr  = 3.21b-1020 for value 141/2

Evaluation took 82.5600 seconds (82.5600 elapsed)

...................................................
Test of precision with fpprec : 2048

meanerr = 2.52b-2044
minerr  = 2.02b-2044 for value 5
maxerr  = 3.19b-2044 for value 129/2

Evaluation took 286.5200 seconds (286.5200 elapsed)

...................................................
Test of precision with fpprec : 4096

meanerr = 3.34b-4092
minerr  = 8.64b-4093 for value 11
maxerr  = 5.16b-4092 for value 71

Evaluation took 1191.0000 seconds (1191.0000 elapsed)
...................................................

This is the code for testing the bigfloat precision:

test_gamma_bfloat(prec) :=
block([ratprint : false,
       fpprec : prec,
       fpprintprec : 3,
       count : 0,
       abserr : 0.0,
       maxerr : -1.0,
       maxerrvalue: 0.0,
       minerr : 1.0,
       minerrvalue : 0.0,
       meanerr : 0.0],
  for a:1/2 thru 75 step 1/2 do
  (
    count : count + 1,

    result : gamma(bfloat(a)),
    answer : bfloat(gamma(a)),

    abserr : abs(result - answer)/abs(answer),
    meanerr : meanerr + abserr,
    maxerr : max(maxerr, abserr),
    if maxerr=abserr then maxerrvalue:a,
    minerr : min(minerr, abserr),
    if minerr=abserr then minerrvalue:a
  ),
  meanerr : meanerr/count,
  printf(true,"...................................................~%"),
  printf(true,"~&Test of precision with fpprec : ~A~%~%",fpprec),
  printf(true,"~&meanerr = ~a~%", meanerr),
  printf(true,"minerr  = ~a for value ~A~%", minerr, minerrvalue),
  printf(true,"maxerr  = ~a for value ~A~%", maxerr, maxerrvalue),
  printf(true,"...................................................~%")
);

Dieter Kaiser