Implementation of the Logarithm of the Gamma function

Hello Ray,

I am not sure, but I think that the accurarcy of the numerical routine
gamma-lanczos is the main source of small errors.

One point is, that actually the algorithm of gamma-lanczos and bffac do not
calculate gamma but something like log_gamma. In the last line of the code of
the routines the exponentiation is done.

So by modifing this routins to get log_gamma we do not do an extra call the
Logarithmic function, but leave out the exponentiation.

In your example below the algorithm do something like

exp(log(exp(log_gamma(x)))+ ....

So my remark to use the log_gamma function directly mean that we do the
following calculation:

exp(log_gamma(x) + ...

Dieter Kaiser

-----Urspr?ngliche Nachricht-----
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com] 
Gesendet: Mittwoch, 8. Oktober 2008 23:49
An: Dieter Kaiser
Cc: maxima at math.utexas.edu
Betreff: Re: AW: [Maxima] Implementation of the Logarithm of the Gamma function

Dieter Kaiser wrote:
> Ray wrote:
>> Note that using logs reduces the accuracy of the result, especially for
>> large arguments that otherwise would have fit in a float.
> 
>> Don't know if this really matters or not.
> 
> Perhaps using the Logarithm of the Gamma function is only interesting for
double
> float values. Here we get an overflow for the Gamma function for values
greater
> than about 170.0. With the help of gamma_log we no longer have the problem of
an
> overflow.

It would be nice, perhaps, if we always tried the float version first.
If an overflow happened, we would use the log version.  Or maybe just
check to see if overflow will happen and then use the log version.

But that's something for another day.


> 
> I have done some tests to compare exp(gamma_log(z)) with gamma(z). We get an
> increasing absolute error for the difference of the two result with increasing
> z, but the relative error is almost constant and is within the expected error
of
> double float calculation.

Consider this:

float(beta(170,-170+1/2));
.1360409980339479

exp(log(gamma(170.0))+log(gamma(-170.0+0.5))-log(gamma(170-170+0.5)));
.1360409980339581

So we lose about 2-3 digits.  May not be important.  I didn't check to
see if the source of the errors are from the gamma function, the log
function or from using logs instead of the direct computation.

Ray