I have compared the algorithms implemented for gammln() in the 1st and 2nd
editions of 'Numerical Recipes in C', and found that they are different.

I strongly suspect that the 2nd edition version is incorrect. The difference
is a removal of 'xx-1.0' and having it replaced by a '/x' in the final 'log()'
call. However, in doing this, the value of 'tmp' is incorrectly computed with
'x' instead of 'x-1'.

The program below confirms this error.

It first computes log(x!), and then gammln() using various versions of
gammln(): The 1st edition version; my optimized version based on the 1st
edition code; and finally the 2nd edition version.

It then outputs the differences between log(x!) and the versions of gammln().

Specifically, the output consists of:

1st Ed. The error between gammln() and log(x!) for the 1st edition code
Mine The error between gammln() and log(x!) for the my code
Mine-1st The difference in the errors between mine and the 1st edition code
2nd Ed. The error between gammln() and log(x!) for the 2nd edition code
2nd-1st The difference in the errors between the 1st and 2nd edition code

A sample of the program's output follows the code listing.

Examining it shows that the 2nd edition version of gammaln() produces a
greater error from log(x!) than the 1st edition code. (The output of my
optimized version is consistent with that of the 1st edition version.)

#include <math.h>

/* Control function - computes x! */
double
fact(double x)
{
if (x == 1.0) {
return (x);
}
return (x * fact(x - 1.0));
}


/* gammln as implemented in the
* first edition of Numerical Recipes in C */
double gammln1(double xx)
{
double x,tmp,ser;
static double cof[6]={76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2,-0.5395239384953e-5};
int j;

x=xx-1.0;
tmp=x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<=5;j++) {
x += 1.0;
ser += cof[j]/x;
}
return -tmp+log(2.5066282746310005*ser);
}


/* My optimizations to gammln as implemented in the
* first edition of Numerical Recipes in C */
double
my_gammln(double xx)
{
double tmp, ser;

tmp = xx + 4.5;
tmp -= (xx - 0.5) * log(tmp);

ser = 1.000000000190015
+ (76.18009172947146 / xx)
- (86.50532032941677 / (xx + 1.0))
+ (24.01409824083091 / (xx + 2.0))
- (1.231739572450155 / (xx + 3.0))
+ (0.1208650973866179e-2 / (xx + 4.0))
- (0.5395239384953e-5 / (xx + 5.0));

return (log(2.5066282746310005 * ser) - tmp);
}


/* gammln as implemented in the
* second edition of Numerical Recipes in C */
double gammln2(double xx)
{
double x,y,tmp,ser;
static double cof[6]={76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2,-0.5395239384953e-5};
int j;

y=x=xx;
tmp=x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<=5;j++) ser += cof[j]/++y;
return -tmp+log(2.5066282746310005*ser/x);
}


int
main(int argc, char **argv)
{
int ii;
double lf;
double my;
double g1, g2;

printf(" xx\t 1st Ed.\t Mine\t\t Mine-1st\t 2nd Ed.\t 2nd-1st\n");

for (ii=2; ii < 10; ii++) {
lf = log(fact((double)(ii - 1.0)));
g1 = gammln1((double)ii)- lf;
my = my_gammln((double)ii) -lf;
g2 = gammln2((double)ii) -lf;

printf("%3d:\t% 13.7e\t% 13.7e\t% 13.7e\t% 13.7e\t% 13.7e\n",
ii, g1, my, my - g1, g2, g2 - g1);
}

for (ii=10; ii <= 170; ii += 5) {
lf = log(fact((double)(ii - 1.0)));
g1 = gammln1((double)ii)- lf;
my = my_gammln((double)ii) -lf;
g2 = gammln2((double)ii) -lf;

printf("%3d:\t% 13.7e\t% 13.7e\t% 13.7e\t% 13.7e\t% 13.7e\n",
ii, g1, my, my - g1, g2, g2 - g1);
}

return (0);
}


Program output:

xx 1st Ed. Mine Mine-1st 2nd Ed. 2nd-1st
2: 0.0000000e+00 -4.4408921e-16 -4.4408921e-16 -4.4408921e-16 -4.4408921e-16
3: -3.3306691e-16 -3.3306691e-16 0.0000000e+00 -5.5511151e-16 -2.2204460e-16
4: -4.4408921e-16 0.0000000e+00 4.4408921e-16 4.4408921e-16 8.8817842e-16
5: 2.2204460e-16 2.2204460e-16 0.0000000e+00 -5.5511151e-16 -7.7715612e-16
6: -4.4408921e-16 -4.4408921e-16 0.0000000e+00 -1.0547119e-15 -6.1062266e-16
7: -1.3322676e-15 -1.3322676e-15 0.0000000e+00 1.8818280e-14 2.0150548e-14
8: 1.8651747e-14 1.8651747e-14 0.0000000e+00 8.3544283e-14 6.4892536e-14
9: 8.3266727e-14 8.3044682e-14 -2.2204460e-16 2.3353541e-13 1.5026869e-13
10: 2.3470115e-13 2.3492319e-13 2.2204460e-16 4.9293902e-13 2.5823788e-13
15: 2.8750335e-12 2.8750335e-12 0.0000000e+00 3.8162806e-12 9.4124708e-13
20: 8.6646246e-12 8.6646246e-12 0.0000000e+00 1.0089707e-11 1.4250823e-12
25: 1.6363799e-11 1.6364021e-11 2.2204460e-16 1.8017143e-11 1.6533441e-12
30: 2.4798608e-11 2.4798608e-11 0.0000000e+00 2.6509017e-11 1.7104096e-12
35: 3.3292924e-11 3.3292924e-11 0.0000000e+00 3.4966696e-11 1.6737722e-12
40: 4.1483483e-11 4.1483483e-11 0.0000000e+00 4.3064663e-11 1.5811796e-12
45: 4.9198867e-11 4.9198867e-11 0.0000000e+00 5.0664362e-11 1.4654944e-12
50: 5.6336269e-11 5.6336269e-11 0.0000000e+00 5.7678307e-11 1.3420376e-12
55: 6.2983174e-11 6.2983174e-11 0.0000000e+00 6.4228622e-11 1.2454482e-12
60: 6.9041661e-11 6.9041661e-11 0.0000000e+00 7.0248696e-11 1.2070345e-12
65: 7.4625639e-11 7.4625417e-11 -2.2204460e-16 7.5714546e-11 1.0889067e-12
70: 7.9813045e-11 7.9813045e-11 0.0000000e+00 8.0784268e-11 9.7122310e-13
75: 8.4595664e-11 8.4595220e-11 -4.4408921e-16 8.5506713e-11 9.1104901e-13
80: 8.9005914e-11 8.9005914e-11 0.0000000e+00 8.9842800e-11 8.3688612e-13
85: 9.3033359e-11 9.3033359e-11 0.0000000e+00 9.3884012e-11 8.5065288e-13
90: 9.6851638e-11 9.6851860e-11 2.2204460e-16 9.7522435e-11 6.7079675e-13
95: 1.0029999e-10 1.0029999e-10 0.0000000e+00 1.0102053e-10 7.2053474e-13
100: 1.0356582e-10 1.0356582e-10 0.0000000e+00 1.0421886e-10 6.5303318e-13
105: 1.0665158e-10 1.0665158e-10 0.0000000e+00 1.0726353e-10 6.1195493e-13
110: 1.0948842e-10 1.0948842e-10 0.0000000e+00 1.0997026e-10 4.8183679e-13
115: 1.1216694e-10 1.1216672e-10 -2.2204460e-16 1.1261791e-10 4.5097259e-13
120: 1.1456280e-10 1.1456236e-10 -4.4408921e-16 1.1515944e-10 5.9663385e-13
125: 1.1699397e-10 1.1699397e-10 0.0000000e+00 1.1744294e-10 4.4897419e-13
130: 1.1913825e-10 1.1913803e-10 -2.2204460e-16 1.1958745e-10 4.4919624e-13
135: 1.2115131e-10 1.2115087e-10 -4.4408921e-16 1.2171997e-10 5.6865623e-13
140: 1.2321344e-10 1.2321366e-10 2.2204460e-16 1.2363355e-10 4.2010839e-13
145: 1.2507728e-10 1.2507728e-10 0.0000000e+00 1.2534107e-10 2.6378899e-13
150: 1.2685453e-10 1.2685453e-10 0.0000000e+00 1.2713031e-10 2.7577940e-13
155: 1.2846091e-10 1.2846091e-10 0.0000000e+00 1.2874679e-10 2.8588243e-13
160: 1.3003054e-10 1.3003043e-10 -1.1102230e-16 1.3036772e-10 3.3717473e-13
165: 1.3149182e-10 1.3149182e-10 0.0000000e+00 1.3183055e-10 3.3872904e-13
170: 1.3298052e-10 1.3298040e-10 -1.1102230e-16 1.3338042e-10 3.9990233e-13

Thanks for your report.

I have looked over the gammln() code in both the first and second editions. Both seem to me to be correct, although they will give slightly different results -- both within the claimed accuracy of the approximation.

Both codes are purporting to compute the function ln(gamma(xx)), starting with an approximation for gamma(xx+1). The first edition code does this with:

ln(gamma(xx)) = ln(gamma((xx-1)+1) = ln(gamma(x+1))
where x=xx-1

The second edition code does it with:

ln(gamma (xx)) = ln(gamma(xx+1)/xx)

Moreover, the coefficients in the second edition were adjusted to be accurate in double precision, while the first edition code is accurate only to single precision.

According to your table, the results of the two codes does, in fact, agree to within the limitations of double precision (provided that the double precision coefficients from the second edition are used in the first edition code).

The question remains, of course, whether both of these programs are right or both wrong. The table attempts to find this out by recording the differences between the codes and a computed value of log((n-1)!). The differences for both codes are as large as about 10^(-10) or so, and that just the claimed accuracy of the approximation for gamma(xx+1). So... I think we would have to say that both programs are right, in the sense that they both deliver the advertised accuracy.

The fact that your "optimized" code agrees more closely with our first edition code is not really relevant, since it implements the first edition form of gamma(xx). Had you implemented the second edition approximation of gamma(xx), your conclusions would have been reversed.

Originally posted by Sherlock
I have looked over the gammln() code in both the first and second editions. Both seem to me to be correct, although they will give slightly different results -- both within the claimed accuracy of the approximation.

The first edition's code does implement the formula presented in the text; the second edition's code does not. The adjustment in the code to remove the 'xx-1.0' at the beginning and replace it with '/x' near the end was done improperly (i.e., the 'tmp' term was not adjusted accordingly), and this (by definition) constitutes a bug.

As to accuracy, my results show that the first edition code is more accurate than the second edition code. Had it been the other way around, one could argue that the 'bug' might be considered an improvement.

While the overall results might be within the stated accuracy, the improper implementation of the formula and the decreased accuracy of the results show that the change to the second edition code for gammln() was not only not an improvement, but is a bug.

I believe I have experienced a problem with the incompete gamma function for alpha=2 and x=3 for an example. the enclosed pdf describes it in detail. I have tried the version of gammln from the first edition but that does not fix it. I hope you can read the file complaint.pdf

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "ig.h"
double f(double t)
{
double f = exp ( -t ) * t ;
return (f);
}
double simpson( double (*d)[2], double (*f)(double d), int n)
{
int m = ( ((n+1)&(~1)) == n) ?n+1:n;
/*even number of intervals, odd number of terms */
double h = ( (*d)[1] - (*d)[0]) / (double) (m-1);
double s;
assert( m > 1);
{
double b = (*d)[0];
double e = (*d)[1];
double a;
int i;
int r;
for (i=0, a=b, s=0.0;i<m;i++,a+=h)
{
double fa = (*f)( a );
if ( i == 0 || i == m-1)
r = 1;
else if ( r == 1)
r = 4;
else if (r == 4)
r = 2;
else if (r == 2)
r = 4;
//fprintf(stdout,"%ld %lf %lf %ld\n", i, a ,fa,r);
s += fa * (double) r;
}
s *= h;
s /= 3;
}
return (s);
};
main(int argc,char **argv)
{
double a=2;
double x=3;
double Fx =gamp(a,x);
fprintf(stdout,"a=%lf x=%lf P(a,x)=%lf\n", a,x,Fx);

/* but simpson, for function */
{
double d[2]={.00000001,x};
double CDF2 = simpson( &d, f,200);
fprintf(stdout,"S: %lf\n", CDF2);
}
}
#if 0
panic
a=2.000000 x=3.000000 P(a,x)=1.000000
S: 0.800852
#endif
~

Hi,

I am unable to reproduce your problem. The NR code gives the correct answer, whether gser or gcf is called by gammp.

Saul Teukolsky

Please cut and paste the tiny complete program brief.c that is exactly the gcf function, and run it. I used Microsoft Visual Studio 8. Then read the pdf complaint2.pdf. This problem is not restricted to the two arguments a=2,x=3, but also for larger x.

Dear Lawrence,

Your brief.c is not at all "exactly the gcf function." The for-loop is changed, the line d = an+ d + b; has a * (multiply) in the original, and there may be other differences. The NR code gives the correct answer, I assure you.

Saul Teukolsky

Thanks for your patience. Apparently I can't copy code accurately. The complaint3.pdf shown applies the changes you pointed out. Thanks again. Now that I have blown by credibility with you, I hope you will read my posts in the future. Thanks again.