simple question

[cowen <cowen at attbi dot com>]

Hello all.

This is perhaps a stupid newbie question. I'm implementing gammln, one of 
the numerical recipes
(available from 
http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf.html, chapter 6).

gammln computes the natural log of the gamma function, which should allow 
one to compute
interesting things, like factorials.

Here's the code:

***************** utils.h **************************************

#ifndef UTILS_H
#define UTILS_H

double gammln ( double xx );

#endif // UTILS_H

**************** utils.c **************************************

#include <math.h>

double gammln ( double xx )
{
        int j;
        double x,y,tmp,ser;
        static double cof[6] = {
                76.18009172947146,
                -86.50532032941677,
                24.01409824083091,
                -1.231739572450155,
                0.1208650973866179e-2,
                -0.5395239384953e-5
        };

        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 );
}

**************************** test.c *******************************************
#include <stdio.h>
#include <math.h>
#include "utils.h"

int main ()
{
    int i;
    for ( i = 1; i <= 20; i++ ) {
        printf( "%3d!: %.0f\n", i, exp( gammln(  i + 1.0 ) ) );
    }
    return 0;
}

Here's the output:

  1!: 1
  2!: 2
  3!: 6
  4!: 24
  5!: 120
  6!: 720
  7!: 5040
  8!: 40320
  9!: 362880
 10!: 3628800
 11!: 39916800
 12!: 479001600
 13!: 6227020800
 14!: 87178291200
 15!: 1307674368006
 16!: 20922789888126
 17!: 355687428098601
 18!: 6402373705783494
 19!: 121645100410059440
 20!: 2432902008204834800


**********************************************************

This is obviously wrong. (15! is evidence enough.)

The algorithm, however, is correct -- I lifted it
directly from the numerical recipe, and I'm fairly confident that
the authors have not made an error. I've built and run the example
on different platforms. The result is always the same.

Here is a primitive Makefile:

all: utils.o test.o test.exe

%.o: %.c
	gcc -c -g -Wall -o $@ $^

test.exe:
	gcc -o $@ -g -Wall utils.o test.o


Any suggestions would be most helpful.

Thanks in advance,

Charles.