Adventures of a Programmer: Parser Writing Peril XXXVI

The Logarithm of the Gamma Function over the Complex Plane with a Smooth Principal Value

As an approximation of course (where applicable) done with the Stirling series but that information would have made the header a wee bit too long.

The problem with the Stirling series is that the number of summands depend on the input; there is an ideal number of summands for a given precision. But don’t worry, D.E.G. Hare did all the math for us.

Preliminaries

The last summand will give out all the necessary information. A term T_k , or summand if you will, is T_k =  \frac{ B_{2k} }{ 2k(2k - 1)z ^ {2k} };  and with a little bit of massage (for the details see Hare’s paper[1]. You can find copies of it at Google) we get the final

\displaystyle{ F = \frac{|T_k|}{ \cos\left(\frac{1}{2} * \arg\left(z\right)^{2 * k}\right) } }

(The above might round down to zero, please take care of that)

No, sorry, not the final, we need the reciprocal of this expression \frac{1}{F} and the digits of that result gives us the expectable accuracy of the result
The result in the calc syntax

/* 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));
}

Here epsilon() sets and gets the precision as a reciprocal of the decimal digits of precision and the default of digits() are the number of decimal digits of the integer part of fractions.
This translates to JavaScript with prototyped functions to

 function lnstirling(z,n){
   var head, z2,  zz;
   var PI   = new Complex(Math.PI,0);
   var half = new Complex(1/2,0);
   var sum = new Complex();
   head = (z.sub(half) ).mul( z.log() ).sub(z).add( PI.add(PI).log().div(2) );
   sum = 0;
   zz = z.dup();
   z2 = z.mul(z);;
   for (var k = 0; k < n; k++) {
     sum = sum.add(precomp[k].div(zz));
     zz = zz.mul(z2);
   }
   return head.add(sum);
}

The computationally most costly is the calculation of the Bernoulli numbers in the terms, so it makes sense to cache the terms without the argument once evaluated. I implemented the JavaScript version for native numbers (64 bit doubles) but not only for them, it just simplifies some things and makes them more readable. So I made two static lists, one with the Bernoulli numbers and one with the precomputed Stirling terms.

Little.bernoullitable /*[30]*/ = [
  1, 
  1/6,
 -1/30,
  1/42,
 -1/30,
  5/66,
 -691/2730,
  7/6,
 -3617/510,
  43867/798,
 -174611/330,
  854513/138,
 -236364091/2730,
  8553103/6,
 -23749461029/870,
  8615841276005/14322,
 -7709321041217/510,
  2577687858367/6,
  -26315271553053477373/1919190,
  2929993913841559/6,
 -261082718496449122051/13530,
  1520097643918070802691/1806,
 -27833269579301024235023/690,
  596451111593912163277961/282,
 -5609403368997817686249127547/46410,
  495057205241079648212477525/66,
 -801165718135489957347924991853/1590,
  29149963634884862421418123812691/798,
 -2479392929313226753685415739663229/870,
  84483613348880041862046775994036021/354,
 -1215233140483755572040304994079820246041491/56786730
];
var precomp = [
  new Complex( 1/12,0),
  new Complex(-1/360,0),
  new Complex( 1/1260,0),
  new Complex(-1/1680,0),
  new Complex( 1/1188,0),
  new Complex(-691/360360,0),
  new Complex( 1/156,0),
  new Complex(-3617/122400,0),
  new Complex( 43867/244188,0),
  new Complex(-174611/125400,0),
  new Complex( 77683/5796,0),
  new Complex(-236364091/1506960,0),
  new Complex( 657931/300,0),
  new Complex(-3392780147/93960,0),
  new Complex( 1723168255201/2492028,0),
  new Complex(-7709321041217/505920,0),
  new Complex( 151628697551/396,0),
  new Complex(-26315271553053477373/2418179400,0),
  new Complex( 154210205991661/444,0),
  new Complex(-261082718496449122051/21106800,0),
  new Complex( 1520097643918070802691/3109932,0),
  new Complex(-2530297234481911294093/118680,0),
  new Complex( 25932657025822267968607/25380,0),
  new Complex(-5609403368997817686249127547/104700960,0),
  new Complex( 19802288209643185928499101/6468)
];

The Bernoulli list is public because it is nice to have them at hand elsewhere, too.
The actual code to evaluate the Stirling sum is simple. In calc syntax:

define __CZ__lngstirling(z, n)
{
    local k head sum z2 bernterm zz;
    head = (z - 1 / 2) * ln(z) - z + (ln(2 * pi()) / 2);
    sum = 0;
    bernterm = 0;
    zz = z;
    z2 = z ^ 2;

    if (size(__CZ__precomp_stirling) < n) {
        __CZ__precomp_stirling = mat[n];
        for (k = 1; k <= n; k++) {
            bernterm = bernoulli(2 * k) / (2 * k * (2 * k - 1));
            __CZ__precomp_stirling[k - 1] = bernterm;
            sum += bernterm / zz;
            zz *= z2;
        }
    } else {
        for (k = 1; k <= n; k++) {
            bernterm = __CZ__precomp_stirling[k - 1];
            sum += bernterm / zz;
            zz *= z2;
        }
    }
    return head + sum;
}

I left the code for caching in it but won’t use it for the JavaScript version which follows now.

function lnstirling(z,n){
  var head, z2,  zz;
  var PI   = new Complex(Math.PI,0);
  var half = new Complex(1/2,0);
  var sum = new Complex();
  head = (z.sub(half) ).mul( z.log() ).sub(z).add( PI.add(PI).log().div(2) );
  sum = 0;
  zz = z.dup();
  z2 = z.mul(z);;
  for (var k = 0; k < n; k++) {
    sum = sum.add(precomp[k].div(zz));
    zz = zz.mul(z2);
  }
  return head.add(sum);
}

That is nigh unreadable, hence the additional post of the implementation in the calc syntax.

lngamma

We want the full domain. The full complex plane with only a countable infinite number of exceptions. That’s huuuge! But we can take that area in smaller bytes, although mostly still uncountable infinite!

The gamma function is symmetric: \log\Gamma\left(z\right) = \overline{\log\Gamma\left(\bar{z}\right)} that saves us half of the complex plane to work with. The negative integers together with the origin are singularities and don’t need further handling. The positive integers without the origin can be done to some extend with the factorial function.
The left real line can be done with a log-gamma function which is restricted to real arguments (needs some massage nevertheless, but not much).
The equality \log\left(\Gamma\left(z+1\right)\right) = \log\left(\Gamma\left(z\right)\right) - \log z lets us push the arguments away from the dangerous places to much greener pastures.

In quasi-pseudo-you-wouldn’t-believe-it-is-code:

/* first branch: all things real */
if(z.Im == 0){
  /* First branch, first sub: integers */
  if( isInteger(z) ){
    /* First branch, first sub, first paragraph: 
       single out the singularities */
    if(z <= 0){
      return Infinity; // or an error
    }/* First branch, first sub, first paragraph*/
    else{
      /* Use factorials up to some precision and size of z */
      /* Use a simple, real-line log-gamme up to some limit */
      /* Use log(n * log(n) +1) */
    }
  }/* First branch, first sub */
  else{ 
    /* First branch, first sub, first paragraph:
       positive real line  */
   if(z > 0){
      /* Use a simple, real-line log-gamme up to some limit */
      /* Use log(n * log(n) +1) */  
   }
    /* First branch, first sub, first paragraph:
       negative real line  */
   else {
      /* Use a simple, real-line log-gamme up to some limit */
      /* Use log(n * log(n) +1) */
      /* correct principal value of logarithm for the
         imaginary part of the result */
   }
  }
} /* first branch */
else {
  /* We have an imaginary part and use the Stirling approximation. */
  /* Imaginary part negative, make it positive but don't forget it */
  if(z.Im < 0){
    z = conjugate(z);
    is_conjugated = true;
  }
  /* Evaluate the number of terms needed */
  /* This number is for 64 bit doubles */
  decimal_digits = 15;
  /* vid. D.E.G. Hare */
  d37 = decimal_digits * .37;
  d52 = decimal_digits * .52;
  term_count = Math.ceil(d52);
  if (abs(z) >= d52  && abs(z.Im) >= d37) {
     term_count = ceil(d37);
  }
  /* We need the original value later */
  Z = z;
  /* We want to change the value (greener pastures, remember?)
     but need to remember how much */
  increasedby = 0;
  /* Another fork according to the sign of the real part */
  if(z.Re > 0){
    /* change the value according to some limits given in
       D.E.G. Hare's paper. It is a bit complicated but the
       code should make it clear */
    if (abs(z.Im) >= d37) {
      while (abs(z) < d52 + 1) {
        z++;
        increasedby++;
      }
    } else {
      /* The function R() has been described above */
      tmp = R(z, termcount);
      tmp2 = tmp;
      while (tmp2 < decimal_digits) {
        z++;
        increasedby++;
        tmp2 = R(z, term_count);
        if (tmp2 < tmp) {
          return error;
        }
      }
    }
    /* Compute Stirling sum with given values */
    ret = lngstirling(z, term_count);
    /* Undo value change */
    if (increasedby > 1) {
      for (k = 0; k < increasedby; k++) {
        ret = ret - ln(Z + (k));
      }
    }
    /* Undo conjugate and return result */
  }
  /* z.Re < 0 */ 
  else {
    /* That part is a bit more complicated */
    /* The change of value has to go in the other direction
       So instead of incrementing z and decrementing it later
       we decrement it first and increment it later
       TL;DR: change z++ to z--  */
    /* We can and will use the reflection formula */
    ret = log( Pi / sin( Pi * z)) - lngstirling(1 - z, term_count);
    /* Undo the value change, this time add the terms */
    /* calculate and aply the correction term for the
       imaginary part to go in concordance with the 
       smooth principal value*/
     corr = ceil(z.Re/2 - 3/4 - kroneckerdelta(z.Im)/4);
     ret = ret + 2 * (corr * Pi) * i;
     /* Undo conjugate and return result */
  } 
  /* Hic sunt dracones! */
}

The Mathematics

As I’m known to be too lazy to search my older posts (“Was schert mich mein Geschwätz von gestern?” [What do I care about my chitchat from yesterday?] said Konrad Adenauer once[2])

/*Kronecker delta function */
define kroneckerdelta(i, j)
{
    if (isnull(j)) {
        if (i == 0) {
            return 1;
        } else {
            return 0;
        }
    }
    if (i != j) {
        return 0;
    } else {
        return 1;
    }
}

The kroneckerdelta function is offered in two flavours: the one and the two argument function, the mathematical function on the other hand is of course the two argument function:

\displaystyle{ \delta_{ij} = \begin{cases} 0 &\text{if } i \neq j \\ 1 &\text{if } i=j, \end{cases}  }

And while we are at mathematical typography—the whole thing in one equation gets too muddled but it is quite ok in two:

\displaystyle{ \log\Gamma\left(z\right) = \begin{cases} \log\bigl(\left(z-1\right)!\bigr) &\text{if } \Im z = 0 \land z \not= 0 \land z \in\mathbb{N} \\ \log\bigl(\Gamma\left(z\right)\bigr) &\text{if } \Im z = 0 \land z > 0 \\ \Bigl( \log \frac{\pi}{ \sin\left(\pi z\right)} - \log\bigl(\Gamma\left(z\right)\bigr)\Bigr) + T_c &\text{if } \Im z = 0 \land z < 0 \\ \Bigl(\log \frac{\pi}{ \sin\left(\pi z\right)} - \log\bigl(\Gamma\left(1-z\right)\bigr) \Bigr) + 2i\biggl( \Bigl\lceil \frac{\Re z}{2} - \frac{3}{4} - \frac{\delta_{\Im z,0}}{4} \Bigr\rceil \pi \biggl) &\text{if } \Im z \not= 0 \land \Re z < 0 \end{cases}  }

With T_c being

\displaystyle{ T_c = \begin{cases}\Re z_r - \pi i & \text{if } |z| \leq \frac{1}{2}\\ \Re z_r + \bigl\lceil \left|z\right| \bigr\rceil \pi i & \text{if } \bigl\lfloor\left|\Re z\right| \bigr\rfloor \text{ odd} \\ \Re z_r + \bigl\lceil\left|z\right| - \frac{1}{2} - 1\bigr\rceil \pi i & \text{if }  \bigl\lfloor\left|\Re z\right| \bigr\rfloor \text{ even } \land  \bigl\lfloor\left| z\right| \bigr\rfloor \text{ even }\\ \Re z_r + \bigl\lceil\left|z\right| - \frac{1}{2}\bigr\rceil \pi i & \text{if }  \bigl\lfloor\left|\Re z\right| \bigr\rfloor \text{ even } \land  \bigl\lfloor\left| z\right| \bigr\rfloor \text{ odd } \end{cases} }

Here z_r =  \Bigl( \log \frac{\pi}{ \sin\left(\pi z\right)} - \log\bigl(\Gamma\left(z\right)\bigr)\Bigr)

Hope that helped, no, wait, better: I hope I got it all right 😉

BTW: tried the new WordPress editor, the blue one. Some things work better, some not, some are worse but it is a bit cleaner. It is, most probably, optimized for use with tablets and phones but one can use it with a desktop, too. It is not that bad.

[1] Hare, David Edwin George. Computing the principal branch of log-Gamma. Journal of Algorithms, 1997, 25. Jg., Nr. 2, S. 221-236.
[2] Quoted in Discussion : Mastering the Skills of Moderation (2009) by Horst Hanisch, p. 91

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s