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