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