Accuracy of numerical calculation of the Gamma and Beta function
Subject: Accuracy of numerical calculation of the Gamma and Beta function
From: Raymond Toy
Date: Fri, 10 Oct 2008 16:49:28 -0400
Dieter Kaiser wrote:
> 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
It's unfortunate that we only have about 40 bits of precision with
53-bit floats.
>
> 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
This makes me wonder if we shouldn't replace gamma lanczos with the same
algorithm as bffac. These look pretty good.
> To get a bit more precision for gamma-lanczos we have to get a new set of
> parameters for the approximation algorithm.
I used to have some code to generate new coefficients. But it was
written in some weird language like bc or dc. The link describing the
algorithm describes pretty well how to generate new coefficients.
I vaguely remember trying different parameters, but it's tricky.
Convergence is better when the parameter gets larger, but then the
accuracy is reduced because the coefficients become very large, causing
more roundoff or something. With a smaller parameter, the coefficients
are more uniform, but then you need more terms, causing more roundoff.
Anyway, these are interesting results.
Ray