1 /**
2 @file
3 @ingroup nr
4 */
5 #include <bpm/bpm_messages.h>
6 #include <bpm/bpm_nr.h>
7
8 /**
9 Calculates the logaritm of the gamma function
10 ln[gamma(xx)]. NR C6.1, p 214
11 supposed to be correct to double precision
12
13 @param xx the argument
14 @return the value of ln[gamma(xx)]
15 */
nr_gammln(double xx)16 double nr_gammln( double xx ) {
17
18 double x,y,tmp,ser;
19 int j;
20 static double cof[6]={76.18009172947146, -86.50532032941677,
21 24.01409824083091, -1.2317395572450155,
22 0.1208650973866179e-2,-0.5395239384953e-5};
23
24 if( xx == 0. ) {
25 bpm_error( "Argument is 0 in nr_gammln(...)", __FILE__, __LINE__ );
26 return -DBL_MAX;
27 }
28
29 // Function also undefined for negative integers
30 if( xx < 0. && nr_is_int( xx ) ) {
31 bpm_error( "Function domain error for nr_gammln(...)",
32 __FILE__, __LINE__ );
33 return -DBL_MAX;
34 }
35
36 y=x=xx;
37 tmp=x+5.5;
38 tmp -= (x+0.5)*log(tmp);
39 for (j=0; j<=5;j++) ser += cof[j]/++y;
40
41 return -tmp+log(2.5066282746310005*ser/x);
42 }
43