Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- 2 // 2 // 3 // ------------------------------------------- 3 // ----------------------------------------------------------------------- 4 // HEP Random 4 // HEP Random 5 // --- HepStat::gammln 5 // --- HepStat::gammln --- 6 // method implementation 6 // method implementation file 7 // ------------------------------------------- 7 // ----------------------------------------------------------------------- 8 8 9 // =========================================== 9 // ======================================================================= 10 // M. Fischler - moved the gammln from Rand 10 // M. Fischler - moved the gammln from RandPoisson to here. 01/26/00 11 // =========================================== 11 // ======================================================================= 12 12 13 #include "CLHEP/Random/Stat.h" 13 #include "CLHEP/Random/Stat.h" 14 #include <cmath> 14 #include <cmath> 15 15 16 namespace CLHEP { 16 namespace CLHEP { 17 17 18 double HepStat::gammln(double xx) { 18 double HepStat::gammln(double xx) { 19 19 20 // Returns the value ln(Gamma(xx) for xx > 0. 20 // Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for 21 // xx > 1. For 0 < xx < 1. the reflection form 21 // xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first. 22 // (Adapted from Numerical Recipes in C. Rela 22 // (Adapted from Numerical Recipes in C. Relative to that routine, this 23 // subtracts one from x at the very start, and 23 // subtracts one from x at the very start, and in exchange does not have to 24 // divide ser by x at the end. The results ar 24 // divide ser by x at the end. The results are formally equal, and practically 25 // indistinguishable.) 25 // indistinguishable.) 26 26 27 static const double cof[6] = {76.18009172947 27 static const double cof[6] = {76.18009172947146,-86.50532032941677, 28 24.01409824083091 28 24.01409824083091, -1.231739572450155, 29 0.120865097386617 29 0.1208650973866179e-2, -0.5395239384953e-5}; 30 int j; 30 int j; 31 double x = xx - 1.0; 31 double x = xx - 1.0; 32 double tmp = x + 5.5; 32 double tmp = x + 5.5; 33 tmp -= (x + 0.5) * std::log(tmp); 33 tmp -= (x + 0.5) * std::log(tmp); 34 double ser = 1.000000000190015; 34 double ser = 1.000000000190015; 35 35 36 for ( j = 0; j <= 5; j++ ) { 36 for ( j = 0; j <= 5; j++ ) { 37 x += 1.0; 37 x += 1.0; 38 ser += cof[j]/x; 38 ser += cof[j]/x; 39 } 39 } 40 return -tmp + std::log(2.5066282746310005*se 40 return -tmp + std::log(2.5066282746310005*ser); 41 } 41 } 42 42 43 } // namespace CLHEP 43 } // namespace CLHEP 44 44 45 45 46 46