Accuracy of numerical calculation of the Gamma and Beta function

I have done some tests concerning the numerical precision of calculating the
Gamma and Beta function. The main idea is to test the floating point evaluation
against the result for an integral or half integral value. All tests with GCL
2.6.8.

First the Gamma function for double float precison (this is a test of the
function gamma-lanczos):

Range of values for gamma(a) is 1/2 ... 75 step 1/2
Compare r1 : float(gamma(a)) with r2 : gamma(float(a)) 

The maximum of the relative error abs(r1-r2)/abs(r1):

maxerr = 4.96e-14

Now the Gamma function for bigfloat precsion (this is a test of the function
bffac):

fpprec  maxerr

  16    9.16b-15
  32    1.39b-30
  64    4.85b-62
 128    7.67b-126
 256    1.92b-253
 512    2.53b-509
1024    3.21b-1020

Second the Beta function (we are using the alogrithm of the functions
gamma-lanczos, log-gamma-lanczos, bffac and bfloat-log-gamma):

Range of values for beta(a,b) is 1/2 ... 75 step 2 for a and b
Compare r1 : float(gamma(a)*gamma(b)/gamma(a+b)) with
r2 : gamma(float(a))*gamma(float(b))/gamma(float(a+b)) and
r3 : exp(log_gamma(float(a))+log_gamma(float(b))-log_gamma(float(a+b)).

Results for double float precision:

maxerr = 1.88e-13 (formula r2)
maxerr = 2.29e-13 (formula r3)

Result for bigfloat precision:

fpprec  maxerr r2  maxerr r3
16      3.24b-14   3.58b-14
32      3.84b-30   3.99b-30
64      7.87b-62   9.48b-62

Conclusion:

There is no significant difference beetween the different algorithm. The
algoritm using log_gamma seems to be sligthly less accurate. Perhaps something
in the routines for log_gamma can be done better.

In bigfloat precision we can achieve any desired accuracy. We have a systematic
difference of 2 to 4 digits in accuracy. Perhaps we can give the function bffac
a bit extra precision to get the missing digits.

To get a bit more precision for gamma-lanczos we have to get a new set of
parameters for the approximation algorithm.

Dieter Kaiser