The motivation for this worksheet comes through a discussion about coding a complex Gamma function
for MPFR (which covers the real case) using a library for complex numbers. The suggestion was to use
a general solution due to Spouge - it seems to be quite simple and it works in MPFR for the real case,
some according code is given by Wilder.
However I did not understand all the steps in Wilder's code (well, MPFR needs ugly looking code, not
my personal taste) and more irritating I was not able to get reasonable results with Maple (no, I have
not tried to use MPFR and am afraid of difficulties in doing that on Windows) in a systematic way.
Spouge's formula is
/ /a - 1 \\
| |----- ||
(z + 1/2) | | \ c(k) ||
z! = (z + a) exp(-z - a) sqrt(2 Pi) |1 + | ) -----|| + epsilon
| | / z + k||
| |----- ||
\ \k = 1 //
(k - 1) (k - 1/2)
(-1) (a - k) exp(a - k)
where c(k) := ---------------------------------------
(k - 1)! sqrt(2 Pi)
for z in the right half plane, where a is an integer uniform for all z
(of size somewhat larger than digits used, a ~ 1.26 * Digits is enough).
Through looking at examples first and then using asymptotics it is shown, how many summands
are needed and how too choose an appropriate working precision to achieve a desired precision.
That is not very smooth and elegant - but more or less natural and how I stumbled to implementations.
All that is not faster than Maple (so you can stop here, if you want to see this). But it was fighting
and having fun ... at least in the end. Hope it is of interest for others as well.
Download 102_Spouge_Gamma.mws or view as pdf print