# The Gamma Function with Spouge’s Approximation

It was promised, by me, here, personally, to implement Spouge’s approximation of
Euler’s Gamma-function in the script language offered by David I. Bell’s
arbitrary precision calculator calc.
I did.
For a short oversight over the algorithm use, see Spouge’s_approximation over at Wikipedia. WordPress’ code colouring has no template for calc-scripting, so I just use C++, it’s similar enough.

Spouge’s approximation is:
$\displaystyle{ \Gamma(z+1) = (z+a)^{z+1/2} e^{-(z+a)} \left[ c_0 + \sum_{k=1}^{a-1} \frac{c_k}{z+k} + \varepsilon_a(z) \right] }$
Let’s see what we have here:

• A leading term $(z+a)^{z+1/2} e^{-(z+a)}$, precomputable
• The first coefficient $c_0 = \sqrt{2 \pi}$ is a constant, precomputable
• All further coefficients $c_k=\frac{(-1)^{k-1}}{(k-1)!}(-k+a)^{k-1/2}e^{-k+a}$ are constants, too, precomputable

To keep those precious coefficients that we calculated at the sweat of our brows, we put them a in a global list for easy recycling.

/*
Cache for the variable independent coefficients in the sum.
*/
global Ck = list();


We calculate the natural logarithm of the function $\log\left(\Gamma\left(z\right)\right)$ —which is different from $\log\Gamma\left(z\right)$ ! —for the right half of the complex plane, that is $\mathfrak{R}\left(z\right) > 0$.

define lngammarp(z)
{
local epsilon accuracy a  factrl fact1 sum coeffk n ret;


It seems that in some versions of Calc the placement of the left curly bracket is not as free as it should be but must be where I put it and as you can see: the variables must be declared before use in the way seen above, i.e. without commas.

Next we calculate the needed accuracy in decimal digits plus the obligatory angstzuschlag[1]. The function that can be used to read and set the current precision is epsilon(). It prints the absolute value of the smallest representable number if it is called without an argument. Another useful function is digits() that returns the number of the digits of the argument and both together digits( 1 / epsilon() ) returns the number of decimal digits, so
epsilon( 10^(- digits( 1 / epsilon() )+1 ) )  does not change anything. The coefficients are computed with a certain precision which we want to keep. The precision I mean.

  epsilon  = epsilon();
accuracy = digits(int(1/epsilon)) + 3;
if(size(Ck) == 0){
append(Ck,accuracy);
}


The calculation of $a$, as described elsewhere in detail, needs only low precision, so it might be a good idea to reduce 1/epsilon to a small number, calculate $a$ and return to the old, higher precision afterwards.

  epsilon(1e-18);
a = ceil(1.252850440912568095810522965 * accuracy);


Next switch to some higher precision. But how much is some? Well, it is a known problem, the calculation of the precision that is needed in relation to the argument z to compute enough correct digits. The main culprit here is the alternating sum. There is not much that can be done for a sum of arbitrary length. There are some ways for sums of small fixed length, see for example the way the people of the math-library at Boost ( John Maddock, Paul A. Bristow, Hubert Holin and Xiaogang Zhang according to the footer at that page) implemented the Gamma functions by way of the original Lanczos’ approximation in a more tricky way.

But in this simple implementation we just take 50% more for the wanted accuracy for
$abs(z) < 1e100$. With this implementation that means that at least digits(int(1/epsilon) are correct, starting from the beginning, that is the highest digit not the first decimal.

Please keep in mind that calc was not meant to be a highly versatile and still fast arbitrary precision floating point calculator but rather a tool for generating and manipulating primes, at least that’s how I understood it.
There are other tools that fit that description better (e.g.: PARI/GP but be aware that lngamma() in PARI/GP has a different branch-cut so our log(gamma(z)) is not necessarily equal to PARI/GP’s lngamma(z) on the complete complex plane. Yes, I just repeated myself) but the gamma-function is such a generic tool, it just needs to be in every tool-box. As are the Bessel functions, but I’ll leave that for another time and post.

Where have I been? Uh, yes, setting the right amount of precision. And when we have done that we can start to compute the constants:

  epsilon(epsilon*10^(-(digits(1/epsilon)//2))); /*integer division*/
fact1 = (ln(z + a) * (z + 0.5)) - (z + a) - ln(z); /* Leading term */
factrl = 1.0;          /* The factorial */
sum = sqrt(2*pi());    /* The first coefficient c_0  */


Calc has integer division. Normal floating point division is done with a single slash 3/2 == 1.5 and integer division with two slashes 3//2 == 1.

Before we are able to fill the list with our coefficients we have to clean up from the work done before, for it could have been done for some lower precision.

  /*
TODO: delete coefficients if precision is lower than the last call.
CLOSED: memory savings not significant, 640k is enough for
everybody!
*/
if(size(Ck) != 0 && Ck[0] < accuracy){
while(size(Ck) != 1){
remove(Ck);
}
Ck[0] == accuracy;
}


The list Ck is a linked list, remove() removes the last entry in the list. The list does not grow very large and is filled with just integers, so no headaches here or cursing of the authors.
The computation of rest of the the coefficients is rather simple and nearly straight-forward.

  for (n = 1; n < a; n++){
z++;
if(size(Ck) == n){
coeffk   = ( ( (a - n)^(n - 0.5) ) * exp(a - n) )/factrl;
append(Ck,coeffk);
} else {
coeffk = Ck[n];
}
sum    += ( coeffk / (z ) );
factrl *= -n; /* the "alternating" in "alternating sum" */
}


That was all! Now put it all together, reset the precision, and return the result.

  ret = ln(sum) + fact1; /* ln() is the natural logarithm base e, log() is base 10!  */
epsilon(epsilon);
return ret;
}


For the left half of the complex plane, some needed checks and balances and the obligatory etc. we build a wrapper. Two wrappers to be exact. The first one is $\log\left(\Gamma\left(z\right)\right)$ for $\{z\vert z\in\mathbf{C}\vee-z\not\in\mathbf{N}\}$.

define lngamma(z)
{
local ret eps;
/* z \in \mathbf{Z}*/
if(isint(z)){
/*The gamma-function has poles at zero and the negative integers*/
if(z <=0)
return 1/0;
else{
/* may hold up accuracy a bit longer, but YMMV */
if(z < 20)
return ln( (z+1)! );
else
return lngammarp(z);
}
}
else{
if(re(z) > 0){
return lngammarp(z);
}
else{


The reflection formula is known to be $\Gamma(1-z) \Gamma(z) = \frac{\pi}{\sin(\pi z)}$.

        eps = epsilon();
epsilon(eps*1e-3);
ret = ln(pi()) - ln( sin(pi() *z ) )   -  lngammarp(1-z)  ;
epsilon(eps);
return ret;
}
}
}


A very simple wrapper for the normal Gamma-function:

define gamma(z)
{
if(isint(z)){
if(z <=0)
return 1/0;
else{
/* may hold up accuracy a bit longer, but YMMV */
if(z < 20)
return (z+1)!*1.0;
else
return exp(lngammarp(z));
}
}
else{
if(re(z) > 0){
return exp(lngammarp(z));
}
else{
return pi()/( sin ( pi() * z )* exp(lngammarp(1-z) ) );
}
}
}


which can be shortened to

define gamma(z)
{
return exp(lngamma(z));
}


and we are done!

Next: The logarithmic derivatives of the Gamma-function generalized as the polygamma-function. I will implement it by way of the Hurwitz-zeta function, the two-variable zeta function, with the exception of the digamma function which will get a special treatment. No, I will not use the already computed Spouge coefficients for the digamma function, I found a different formula faster.

[1] “Angstzuschlag” is the colloquial geman word for “corrosion allowance”.