Implementation of the Logarithm of the Gamma function

Some days ago I have reported some work on the implementation of the Incomplete
Beta function. By the way I have recognized that for the numerical evaluation of
the Beta function the transformation to the Gamma function is used: 

beta(a,b)=gamma(a)*gamma(b)/gamma(a+b).

Especially for float arguments I think it is better to use the Logarithm of the
Gamma function for evaluation:

beta(a,b)=exp(gamma_log(a)+gamma_log(b)-gamma(a+b)).

So in a first step I have rewritten the algorithm of gamma-lanczos, bffac and
cbffac to calculate the Logarithm of the Gamma function. With these functions we
can calculate the quotient of Gamma functions for a much wider range of values
and more accurate. The most complicated part was to find a consistent
calculation for a negative realpart of the argument.

I have committed a first implementation of the function gamma_log(z).

For a full implementation of the Incomplete Beta function the Beta function has
to be extended. 

So what do you think? Should we give more support to the Beta function using the
Logarithm of the Gamma function?

Dieter Kaiser