Adventures of a Programmer: Parser Writing Peril XXXV

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 \log\Gamma() function which is different from \log\left(\Gamma\left(\right)\right) on the left complex half-plane except the real-line which can be ignored because of \Re\left(\log\Gamma\left(\right)\right) =  \log\left(|\Gamma\left(\right)|\right) [1]:

\Im\left(\log\Gamma\left(\right)\right) = \Im\left( \log\left(\Gamma\left(\right)\right) \right) + 2k\left(x\right)\pi = \begin{cases} 2k\left(x\right)\pi & \text{if } \lfloor x\rfloor \text{ is even}\\ \left(2k\left(x\right)\right)\pi & \text{if } \lfloor x\rfloor \text{ is odd} \end{cases}

Now what is the function k 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 \sqrt{\pi} , in the denominator for positive values and as a factor for negative values. the factor is one at x = 0.5.

If you wonder about the complicated branching for the negative real-line you haven’t read the math above where it says \Re\left(\log\Gamma\left(\right)\right) =  \log\left(|\Gamma\left(\right)|\right) . 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.

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