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