simple question
cowen
cowen@attbi.com
Wed Nov 27 16:34:00 GMT 2002
Sorry -- this is an inappropriate list for this question.
I'll take it to gcc-help.
Again, my apologies.
C.
At 09:51 AM 11/27/2002 -0800, cowen wrote:
>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.
>
>
More information about the Gcc
mailing list