Javia blog

The gamma function is the extension of the factorial to real numbers (n! == gamma(n+1)), and it’s a very important function in mathematics. Gamma(x) grows very quickly (similar to x**x), that’s why lgamma(), which is the logarithm of the gamma function, is often used instead. In the discussion below, I’m only considering the lgamma(x) on the positive axis (x>0).

I wanted to have an implementation of lgamma() in Java. I studied these four implementation alternatives:

1) The Lanczos approximation using g=7 and 9 terms, which is used in the GNU Scientific Library. Also the sample code in the wikipedia article uses the same parameters (g=7, 9 terms).

private static final double L9[] = {
    0.99999999999980993227684700473478,
    676.520368121885098567009190444019,
    -1259.13921672240287047156078755283,
    771.3234287776530788486528258894,
    -176.61502916214059906584551354,
    12.507343278686904814458936853,
    -0.13857109526572011689554707,
    9.984369578019570859563e-6,
    1.50563273514931155834e-7
};
private static final double LN_SQRT2PI = 0.9189385332046727418; //log(2*PI)/2

static final double lanczosLGamma9(double x) {
    if (x <= -1) return Double.NaN;
    double a = L9[0];
    for (int i = 1; i < 9; ++i) {
        a+= L9[i]/(x+i);
    }
    return (LN_SQRT2PI + Math.log(a) - 7.) + (x+.5)*Math.log((x+7.5)/Math.E);
}

2) The Lanczos approximation using g=607/128 and 15 terms, used in the Jakarta Mathematics Library and described here.

private static final double[] L15 = {
    0.99999999999999709182,
    57.156235665862923517,
    -59.597960355475491248,
    14.136097974741747174,
    -0.49191381609762019978,
    .33994649984811888699e-4,
    .46523628927048575665e-4,
    -.98374475304879564677e-4,
    .15808870322491248884e-3,
    -.21026444172410488319e-3,
    .21743961811521264320e-3,
    -.16431810653676389022e-3,
    .84418223983852743293e-4,
    -.26190838401581408670e-4,
    .36899182659531622704e-5,
};

static final double lanczosLGamma15(double x) {
    if (x <= -1) return Double.NaN;
    double a = L15[0];
    for (int i = 1; i < 15; ++i) {
        a += L15[i]/(x+i);
    }

    double tmp = x + (607/128. + .5);
    return (LN_SQRT2PI + Math.log(a)) + (x+.5)*Math.log(tmp) - tmp;
}

3) The Stirling approximation with 4 terms.

private static final double
    LC1 = 0.08333333333333333,
    LC2 = -0.002777777777777778,
    LC3 = 7.936507936507937E-4,
    LC4 = -5.952380952380953E-4;

static final double stirlingLGamma(double x) {
    final double
        r1 = 1./x,
        r2 = r1*r1,
        r3 = r1*r2,
        r5 = r2*r3,
        r7 = r3*r3*r1;
    return (x+.5)*Math.log(x) -x + LN_SQRT2PI + LC1*r1 + LC2*r3 + LC3*r5 + LC4*r7;
}

4) The lgamma implementation from SUN’s FDLIBM 5.3.
As this implementation is more complex, and notably uses 61 precomputed double constants, I don’t inline here, but you may find it in the full source code.

I was curious which one of these implementations is the most precise, and which one is the fastest.
I computed an “error table” which shows the error of these functions in ULPs over the integers from 0 to 170. This is the beginning of the error table:

Error in ULPs

x	FDLIBM	Lanczos 15	Lanczos 9	Stirling
1	0	0	32	-1000+
2	0	1	9	-1000+
3	0	-2	10	-1000+
4	0	-1	-1	-1000+
5	0	-1	3	-1000+
6	0	-2	1	-1000+
7	-1	0	2	-1000+
8	0	-1	1	-1000+
9	1	1	4	-1000+
10	1	-1	-1	-463
11	0	0	0	-98
12	1	1	0	-45
13	-1	0	0	-23
14	0	2	0	-12
15	0	2	0	-5
16	1	0	2	-4
17	0	0	0	-1
18	0	-1	2	-1
19	0	0	0	-1
20	0	0	0	1

You may also consult the full table of lgamma error.

The conclusion is that FDLIBM’s implementation is the most precise, as it always has an error of at most 1 ULP.
Lanczos with 15 terms is good, as it has an error of at most 2 ULPs. Lanczos with 9 terms is somewhat worse, especially for small values. And Stirling is bad for arguments under 16, but above 16 has an error similar to Lanczos 15.

So FDLIBM’s lgamma() is the clear winner on precision. What about speed?
Below is the time (ms) needed for computing lgamma for all the integers from 0 to 170, 20000 times.

method	time (ms)
FDLIBM	1082
Lanczos 15	1905
Lanczos 9	1524
Stirling	618

(smaller time is better)

The fastest method is Stirling (but unfortunately it’s not accurate for small arguments), followed by FDLIM’s lgamma(), which is almost twice as fast as Lanczos 15.
So the winner considering both performance and speed is FDLIBM. The downside of FDLIBM’s lgamma() is that it is a bit more complex that the other methods, and also it uses a table of 61 precomputed double values, which is 4 times bigger than the table used by Lanczos 15.

The full source code, containing the four implementations of lgamma(), the error evaluation and speed benchmark is available here:
Gamma.java

This entry was posted on Sunday, November 5th, 2006 at 03:58 and is filed under Uncategorized. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.