libbpm: bpmnr/nr_gammln.c Source File

00001 
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_nr.h>
00007 
00016 double nr_gammln( double xx ) {
00017 
00018   double x,y,tmp,ser;
00019   int j;
00020   static double cof[6]={76.18009172947146, -86.50532032941677,
00021                         24.01409824083091, -1.2317395572450155,
00022                         0.1208650973866179e-2,-0.5395239384953e-5};
00023 
00024   if( xx == 0. ) {
00025     bpm_error( "Argument is 0 in nr_gammln(...)", __FILE__, __LINE__ );
00026     return -DBL_MAX;
00027   }
00028 
00029   // Function also undefined for negative integers
00030   if( xx < 0. && nr_is_int( xx ) ) {
00031     bpm_error( "Function domain error for nr_gammln(...)", 
00032                __FILE__, __LINE__ );
00033     return -DBL_MAX;
00034   }
00035 
00036   y=x=xx;
00037   tmp=x+5.5;
00038   tmp -= (x+0.5)*log(tmp);
00039   for (j=0; j<=5;j++) ser += cof[j]/++y;
00040 
00041   return -tmp+log(2.5066282746310005*ser/x);
00042 }