simple question

[cowen <cowen@attbi.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 (or lnGamma as I've named it) 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 lnGamma ( double xx );

#endif // UTILS_H


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

#include <math.h>

double lnGamma ( 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( lnGamma(  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.