# On Doing It Yourself

### When the output of Mathematica® is correct but of not much use

When I do some calculations on the command line, it is in most cases some quick computations involving numbers with different bases in a precision that is not very high but higher than what the normal calculator offers. Starting up PARI/GP every time is distracting and PARI/GP does not offer an easy way to enter numbers in different bases; dc and bc are slightly more complicated to use—although they are quite excellent in shell-scripts!—and all others are either insufficient in power, usage and/or price so it is calc what I use regularly.

It has some downsides, too. One of these disadvantages is the lack of functions that are included in almost every other arbitrary precision calculator because calc had been written for a special need originally. See the webpage of calc I linked to above for further information.

Another advantage of calc is that the syntax is quite C-like and allows for easy implementation of the functions and algorithms you need. At least that’s what I hope for 😉

The most common implementation of an approximation of the gamma function when arbitrary precision is needed, is the variation of James Stirling’s approximation found by John L. Spouge[1]. It is some time ago, that I implemented it and I had to look it up. The first page that came up after asking Google for some help, was the page of the entry of Literate Prgrams for an implementation in Mathematica®. No, I really cannot remember what I’ve asked Google to get that result. I found out later, that Rosetta Code has a rudimentary but working example implementation in C.

Now, where do I am going to get to with such a long-winded introduction, you might ask, and I will tell you: on the Literate-Programming page I found a short explanation of the algorithm which includes the following gem:

Conservatively omitting the factor $a^{-1/2}$ and rearranging gives us the formula…

They talk about the error and the number of coefficients needed:

(1.0) $\displaystyle{ \epsilon \le a^{-1/2} (2 \pi)^{-(a+1/2)}}$

for which they give the approximation of the parameter $a$ as:

(1.1) $\displaystyle{ a = -\frac{\log \epsilon}{\log 2 \pi}}$

If you enter the inequality into the text field over at wolframalpha.com and tell it to solve for a you’ll get an exact but unusable result. Changing the inequality to an equality you’ll get another, simpler but still unusable result involving the Lambert-W function. If you want to play with it, here is a JavaScript implementation of the approximation on the positive real line (click on it to unfold it):

Math.lambertW = function(x){
/*
On the Lambert W Function (1996)
by R. M. Corless , G. H. Gonnet , D. E. G. Hare , D. J. Jeffrey , D. E. Knuth
*/
var ret = 0.0;
var a   = Math.log(x);   //L1
var b   = Math.log(Math.log(x));//L2
var c   = a-b + b/a;
var l1  = (b*(-2 + b))/(2*a*a);
var l2  = (b*(6-9*b+2*b*b))/(6*a*a*a);
var l3  = (b*(-12+36*b-22*b*b+3*b*b*b))/(12*a*a*a*a);
var l4  = (b*(60-300*b+350*b*b-125*b*b*b+12*b*b*b*b))/(60*a*a*a*a*a*a);

return c+l1+l2+l3+l4;
};


So I decided to do it by hand. Why? Well, if you are seriously questioning it read no further, for what I do is not for you.
But I do not use anything not taught in highschool and even give the needed rules for the logarithms I use here, so you are heartily invited to give it a glimpse. You might even find an error 😉

At first make the inequality an equality because we are interested in the largest possible error only. Expand all of the exponents to make it more readable. For me at least but YMMV, as the saying goes.

(2.1) $\displaystyle{ \epsilon = \frac{1}{\sqrt{a}\,\pi^{a+1/2}}}$

The term $\sqrt{a}$ can be omitted here, that is indeed correct. I’ll give a short sketch of a proof at the end.

(2.2) $\displaystyle{ \epsilon = \frac{1}{2\pi^{a+\frac{1}{2}}}}$

To get rid of the exponent we have to take the logarithm of both sides. Any logarithm will do.

(2.3) $\displaystyle{ \log\epsilon = \log\left(\frac{1}{2\pi^{a+\frac{1}{2}}}\right)}$

Expand it, such that we get the exponents back we dropped so prematurely above.
Rules: $\log (a/b) = \log a - \log b$ and $\log (a^b) = log(a)\, b$

(2.4) $\displaystyle{ \log\epsilon = \log{1}-\log(2\pi) \left(a+\frac{1}{2}\right)}$

The logarithm of 1 (one) is 0 (zero) and can be omitted for brevity. Now move the constant $-\log{2\pi}$ to the left of the equation.

(2.5) $\displaystyle{ -\frac{\log\epsilon}{\log(2\pi)} = a + \frac{1}{2}}$

Subtract $\frac{1}{2}$ and we are done so far.

(2.6) $\displaystyle{-\frac{log\epsilon}{\log(2\pi)} - \frac{1}{2} = a}$

Because of the conditions $a \in \mathbb{N}$ together with $a > 2$ we can omit the $\frac{1}{2}$, too.

What? Uh, the sketch for the proof for dropping the poor $\sqrt{a}$, yes.
We start again from equation (2.1). To make it more readable we move the square root to the left to get

(3.1) $\displaystyle{ \epsilon\sqrt{a} = \frac{1}{2\pi^{a+\frac{1}{2}}}}$

Take the logarithm of both sides and expand
Rule: $\log(\sqrt[n]{a}) = \frac{\log a}{n}$ which is just a special case of $\log (a^b) = log(a)\, b$

(3.2) $\displaystyle{ \log\epsilon + \frac{\log a}{2} = -\log (2\pi) (a+\frac{1}{2}) }$

Uh, no, move $\sqrt{a}$ back to the right, sorry.

(3.3) $\displaystyle{ \log\epsilon = -\log (2\pi) \left(a+\frac{1}{2}\right) - \frac{\log a}{2} }$

Divide by $-\log (2\pi)$ as in equation (2.5) and expand

(3.4) $\displaystyle{ -\frac{\log\epsilon}{\log (2\pi)} = -\frac{2}{\log a\log(2\pi)} + a+\frac{1}{2} }$

Subtract $\frac{1}{2}$ and rearrange the right side to get

(3.5) $\displaystyle{ -\frac{\log\epsilon}{\log (2\pi)} -\frac{1}{2} = a - \frac{2}{\log a\log(2\pi)}}$

Now we can easily see by way of the inequality

(3.6) $\displaystyle{ a - \frac{2}{\log a\log(2\pi)} \le a}\quad\text{with}\quad a\in\mathbb{N}\, \wedge\, a>2$

that the term $\sqrt{a}$ had been dropped conservatively but rightfully $\blacksquare$

If you waited the whole time for the calc-script: you will find it in the repository if the authors of calc accept it; here, in another posting, otherwise.

[1] Spouge, John L. (1994). Computation of the Gamma, Digamma, and Trigamma Functions. SIAM Journal on Numerical Analysis 31 (3): 931–000.