Added some basic function and needful utils to the library of Little. Most of it from my old pragmath package. Of larger interest might be a port of GNU-coreutils’ `factor`

and a port of Boost’s Gamma functions for real arguments. I wanted to base the gamma functions for complex arguments on these but they are lacking too much of accuracy for some arguments (near one and two), mainly in the log-gamma function.

I forgot if I already wrote about it and am too lazy to look it up (as if I have written thousands of posts π ): the problem is the principal value of the function which is different from on the left complex half-plane except the real-line which can be ignored because of [1]:

Now what is the function you might ask. Well, it is in the paper π

But I have done not only the long-winded arbitrary precision function implemented in `calc`

but also a small, readable version for testing purposes. Please look under the fold.

define kroneckerdelta(i, j) { if (isnull(j)) { if (i == 0) { return 1; } else { return 0; } } if (i != j) { return 0; } else { return 1; } } define lngammasmall(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; } ret = (__CZ__lngamma_lanczos(z)); return ret; } }

Where` __CZ__lngamma_lanczos()`

is just a small low-precision complex gamma function. A short test in a 10×10 rectangle on the complex plane around the origin comparing the results to GP/PARI gave no reason to doubt the correctness. Ok, not *much* reason π

Some optimizations have been left out, e.g.: the gamma functions of all half integers are rational multiples of , in the denominator for positive values and as a factor for negative values. the factor is one at .

If you wonder about the complicated branching for the negative real-line you haven’t read the math above where it says . Yepp, it has the *absolute value* of the Gamma function on the right hand side and yepp, I didn’t read it with the necessary care the first time, too which had cost me some hours of wonder(ing ) when I tested it against GP/PARI.

You might also take a look at the file `gammazeta.py`

in `mpmath`

. Have fun trying to find out how they did it. And they even wrote a comment above it describing it! π

The main problem with every such approximation is the canceling in the sums. You need some extra precision to catch them and if you work with a fixed precison…

I have implemented the multi-precision version with the Stirling approximation in `calc`

which is quite complex but D.E.G. Hare’s paper is not only based on this approximation (no it is based on the mathematics, of course but I hope you know what I mean) it also shows every necessary information you need to implement it.

My first idea was to use Spouge’s approximation but the logarithmic version has a different branch from the Stirling approximation. Or I made something wrong in the implementation. But if you look into the above mentioned mpmath file you’ll find that they gave up, too π

Next post: a Javascript implementation of the lg-gamma function with a smoothed principal value.

[1] Hare, David Edwin George. Computing the principal branch of log-Gamma. Journal of Algorithms, 1997, 25. Jg., Nr. 2, S. 221-236.