#include "CalRecon/gamma.h"
#include <cmath>
Go to the source code of this file.
Functions | |
double | LnGamma (double z) |
double | Gamma (double a, double x) |
double | GamCf (double a, double x) |
double | GamSer (double a, double x) |
|
Computation of the incomplete gamma function P(a,x) via its continued fraction representation. The algorithm is based on the formulas and code as denoted in Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.). --- Nve 14-nov-1998 UU-SAP Utrecht Definition at line 80 of file gamma.cxx. References LnGamma(). Referenced by Gamma().
00081 { 00082 00083 int itmax = 100; // Maximum number of iterations 00084 double eps = 3.e-7; // Relative accuracy 00085 double fpmin = 1.e-30; // Smallest double value allowed here 00086 00087 if (a <= 0 || x <= 0) return 0; 00088 00089 double gln = LnGamma(a); 00090 double b = x+1-a; 00091 double c = 1/fpmin; 00092 double d = 1/b; 00093 double h = d; 00094 double an,del; 00095 for (int i=1; i<=itmax; i++) { 00096 an = double(-i)*(double(i)-a); 00097 b += 2; 00098 d = an*d+b; 00099 if (fabs(d) < fpmin) d = fpmin; 00100 c = b+an/c; 00101 if (fabs(c) < fpmin) c = fpmin; 00102 d = 1/d; 00103 del = d*c; 00104 h = h*del; 00105 if (fabs(del-1) < eps) break; 00106 //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl; 00107 } 00108 00109 double v = exp(-x+a*log(x)-gln)*h; 00110 return (1-v); 00111 } |
|
Computation of the incomplete gamma function P(a,x) The algorithm is based on the formulas and code as denoted in Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.). --- Nve 14-nov-1998 UU-SAP Utrecht Definition at line 60 of file gamma.cxx. References GamCf(), and GamSer().
|
|
Computation of the incomplete gamma function P(a,x) via its series representation. The algorithm is based on the formulas and code as denoted in Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.). --- Nve 14-nov-1998 UU-SAP Utrecht Definition at line 124 of file gamma.cxx. References LnGamma(). Referenced by Gamma().
00125 { 00126 int itmax = 100; // Maximum number of iterations 00127 double eps = 3.e-7; // Relative accuracy 00128 00129 if (a <= 0 || x <= 0) return 0; 00130 00131 double gln = LnGamma(a); 00132 double ap = a; 00133 double sum = 1/a; 00134 double del = sum; 00135 for (int n=1; n<=itmax; n++) { 00136 ap += 1; 00137 del = del*x/ap; 00138 sum += del; 00139 if (fabs(del) < fabs(sum*eps)) break; 00140 //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl; 00141 } 00142 double v = sum*exp(-x+a*log(x)-gln); 00143 return v; 00144 } |
|
Computation of ln[gamma(z)] for all z>0. The algorithm is based on the article by C.Lanczos [1] as denoted in Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.). [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86. The accuracy of the result is better than 2e-10. --- Nve 14-nov-1998 UU-SAP Utrecht Definition at line 25 of file gamma.cxx. Referenced by GamCf(), and GamSer().
00026 { 00027 00028 if (z<=0) return 0; 00029 00030 // Coefficients for the series expansion 00031 00032 double c[7]= { 2.5066282746310005, 76.18009172947146, -86.50532032941677 00033 ,24.01409824083091, -1.231739572450155, 0.1208650973866179e-2 00034 ,-0.5395239384953e-5}; 00035 00036 double x = z; 00037 double y = x; 00038 double tmp = x+5.5; 00039 tmp = (x+0.5)*log(tmp)-tmp; 00040 double ser = 1.000000000190015; 00041 for (int i=1; i<7; i++) { 00042 y += 1; 00043 ser += c[i]/y; 00044 } 00045 double v = tmp+log(c[0]*ser/x); 00046 return v; 00047 } |