This is the mail archive of the
gcc-help@gcc.gnu.org
mailing list for the GCC project.
simple question
- From: cowen <cowen at attbi dot com>
- To: gcc-help at gcc dot gnu dot org
- Date: Wed, 27 Nov 2002 10:11:18 -0800
- Subject: simple question
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.