]>
Commit | Line | Data |
---|---|---|
9419eeef | 1 | // PythiaStdlib.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2010 Torbjorn Sjostrand. | |
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. | |
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. | |
5 | ||
6 | // Function definitions (not found in the header) for the Gamma function. | |
7 | ||
8 | #include "PythiaStdlib.h" | |
9 | ||
10 | namespace Pythia8 { | |
11 | ||
12 | //========================================================================== | |
13 | ||
14 | // The Gamma function for real arguments, using the Lanczos approximation. | |
15 | // Code based on http://en.wikipedia.org/wiki/Lanczos_approximation | |
16 | ||
17 | double GammaCoef[9] = { | |
18 | 0.99999999999980993, 676.5203681218851, -1259.1392167224028, | |
19 | 771.32342877765313, -176.61502916214059, 12.507343278686905, | |
20 | -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7}; | |
21 | ||
22 | double GammaReal(double x) { | |
23 | ||
24 | // Reflection formula (recursive!) for x < 0.5. | |
25 | if (x < 0.5) return M_PI / (sin(M_PI * x) * GammaReal(1 - x)); | |
26 | ||
27 | // Iterate through terms. | |
28 | double z = x - 1.; | |
29 | double gamma = GammaCoef[0]; | |
30 | for (int i = 1; i < 9; ++i) gamma += GammaCoef[i] / (z + i); | |
31 | ||
32 | // Answer. | |
33 | double t = z + 7.5; | |
34 | gamma *= sqrt(2. * M_PI) * pow(t, z + 0.5) * exp(-t); | |
35 | return gamma; | |
36 | ||
37 | } | |
38 | ||
39 | //========================================================================== | |
40 | ||
41 | } // end namespace Pythia8 |