CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: gammln.cc,v 1.3 2003/08/13 20:00:12 garren Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- HepStat::gammln --- 00007 // method implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // M. Fischler - moved the gammln from RandPoisson to here. 01/26/00 00012 // ======================================================================= 00013 00014 #include "CLHEP/Random/Stat.h" 00015 #include "CLHEP/Random/defs.h" 00016 #include <cmath> 00017 00018 using namespace std; 00019 00020 namespace CLHEP { 00021 00022 double HepStat::gammln(double xx) { 00023 00024 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 00025 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 00026 // (Adapted from Numerical Recipes in C. Relative to that routine, this 00027 // subtracts one from x at the very start, and in exchange does not have to 00028 // divide ser by x at the end. The results are formally equal, and practically 00029 // indistinguishable.) 00030 00031 static double cof[6] = {76.18009172947146,-86.50532032941677, 00032 24.01409824083091, -1.231739572450155, 00033 0.1208650973866179e-2, -0.5395239384953e-5}; 00034 int j; 00035 double x = xx - 1.0; 00036 double tmp = x + 5.5; 00037 tmp -= (x + 0.5) * log(tmp); 00038 double ser = 1.000000000190015; 00039 00040 for ( j = 0; j <= 5; j++ ) { 00041 x += 1.0; 00042 ser += cof[j]/x; 00043 } 00044 return -tmp + log(2.5066282746310005*ser); 00045 } 00046 00047 } // namespace CLHEP 00048 00049