On the Principal Branch of the Log-Gamma Function [updated]

[UPDATE at Mo 29. Jul 17:09:30 CEST 2013] fixed a small bug.
We did a fast glimps—albeit with working code—at the implementation of the gamma function with Spouge’s formula lastly. This code implements the logarithm of the gamma function $\log (\Gamma z)$ with the principal branches of the complex logarithm. Another function, named $\log\Gamma z$ uses different principal branches that are more useful. You will find that it is this function that is implemented in most programs that offer arbitrary precision.
To go from $\log (\Gamma z)$ to $\log\Gamma z$ is not easy:
It is equivalent in the strip

$\displaystyle{ \log\Gamma z = \log(\Gamma z) \quad \text{with} \quad 0 < \Re z \le 2 \vee \lvert \Im z\rvert \le \frac{7}{2} }$

The imaginary part of negative values are simply multiples of $\pi$

$\displaystyle{ \log\Gamma x = \log (\Gamma z) + 2 i \pi \biggl\lceil\frac{x}{2}-1\biggr\rceil \quad \text{with} \quad x<0 \vee -x\not\in\mathbb{N}}$

The interesting part comes with the rest of the complex plane:

$\displaystyle{\log\Gamma x = \log (\Gamma z) +2 i \pi k(z)}$

The function $k(z)$ returns an integer. To compute these size of the aforementioned integer the following simple integral has to be solved:

$\displaystyle{k(z) = \int_0^z \theta \bigl( -\Re(\Gamma t) \bigr) \bigl\lvert \Im\bigl(\Gamma(t)\psi(t) \bigr) \bigr\rvert \delta\bigl( \Im(\Gamma t) \bigr) \mathrm{d}t}$

with $\theta (z)$ the unit step function

$\displaystyle{\theta (x) = \begin{cases} 0& x< 0\\ 1&<\ge 0 \end{cases} }$

with $\psi (z)$ the first logarithmic derivative of the Gamma function

$\displaystyle{ \psi (z) = \frac{\mathrm{d}}{\mathrm{d}z}\Gamma z = \frac{\Gamma '(z)}{\Gamma(z)} }$

and finally $\delta (x)$ a kind of a derivative of the Heaviside step function and has the following limit

$\displaystyle{\delta(x)\lim_{\epsilon\to 0}\frac{\sin\frac{x}{\epsilon}}{\pi x} }$

All that information should suffice to handle the integral. Not that it wouldn’t be a good training, but others have solved the problem already, namely D. E. G. Hare in “Computing the Principal Branch of log-Gamma” published in the Journal of Algorithms 25(2):221-236 in 1997.

Attached to a simple Gamma function with the Lanczos coefficients

define kroneckerdelta(i,j){
if(isnull(j)){
if(i==0) return 1;
else return 0;
}
if(i!=j) return 0;
else return 1;
}
define __CZ__lngamma_lanczos(z){
local a k g zghalf lanczos;
mat lanczos[15] = {
9.9999999999999709182e-1,
5.7156235665862923516e1,
-5.9597960355475491248e1,
1.4136097974741747173e1,
-4.9191381609762019978e-1,
3.3994649984811888638e-5,
4.6523628927048576010e-5,
-9.8374475304879566105e-5,
1.5808870322491249322e-4,
-2.1026444172410489480e-4,
2.1743961811521265523e-4,
-1.6431810653676390482e-4,
8.4418223983852751308e-5,
-2.6190838401581411237e-5,
3.6899182659531626821e-6
};
g = 607/128;
z = z-1;
a = 0;
for(k = 12; k >= 1; k--){
a += lanczos[k]/(z+k);
}
a += lanczos[0];
zghalf = z + g + .5;
return ( ln(sqrt(2*pi())) + ln(a)  -zghalf  ) + (z+.5)*ln( zghalf );
}
define lngamma(z){
local ret eps flag corr;
flag = 0;
if(isint(z)){
if(z <=0)
return newerror("gamma(z): z is a negative integer");
else{
/* may hold up accuracy a bit longer, but YMMV */
if(z < 20)
return ln((z-1)!);
else
return __CZ__lngamma_lanczos(z);
}
}
else{
if(re(z)<0){
if(im(z) && im(z)<0){
z = conj(z);
flag = 1;
}

ret = ln( pi() / sin(pi() *z ) ) -__CZ__lngamma_lanczos(1-z);

if(!im(z)){
if(abs(z) <= 1/2)
ret = re(ret) - pi()*1i;
else if( isodd(floor(abs(re(z)))) ){
ret = re(ret) + ( ceil(abs(z))   * pi() * 1i);
}
else if( iseven(floor(abs(re(z)))) ){
/* < n+1/2 */
if(iseven(floor(abs(z)))){
ret = re(ret) + ( ceil(abs(z)-1/2 -1 )   * pi() * 1i);
}
else{
ret = re(ret) + ( ceil(abs(z) -1/2 )   * pi() * 1i);
}
}
}
else{
corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
ret = ret +2*(corr *pi() )*1i;
}
if(flag == 1)
ret = conj(ret);
return ret;
}
return __CZ__lngamma_lanczos(z);
}
}


My problem is, that this trick does not work with the Spouge coefficients. At least not for me. I didn’t find out if it is a bug in my implementation (most probably) or a more technical problem. The paper mentioned above talks about the Stirling series. It uses Bernoulli numbers which are quite expensive to compute and the series is not convergent which makes it quite tricky to find the right number of terms. D.E.G. Hare showed in section 4.1 a trick to evaluate the error in a simple way for $z\in\mathbb{C}\setminus\mathbb{R}$.

/* Compute error for Stirling series for "z" with "k" terms*/
define R(z,k){
local a b ab stirlingterm;
stirlingterm = bernoulli(2*k)/( 2*k*(2*k-1)*z^(2*k));
a = abs( stirlingterm );
b = (cos(.5*arg(z)^(2*k)));
ab = a/b;
/* a/b might round to zero */
if( abs(ab) < epsilon()){
return digits(1/epsilon());
}
return digits(1/(a/b));
}