gamma.h File Reference

Function Documentation

double GamCf	(	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.

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 }

double Gamma	(	double	a,
		double	x
	)

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.

00061 {
00062 
00063    if (a <= 0 || x <= 0) return 0;
00064 
00065    if (x < (a+1)) return GamSer(a,x);
00066    else           return GamCf(a,x);
00067 }

double GamSer	(	double	a,
		double	x
	)

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.

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 }

double LnGamma

(

double

)

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.

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 }

Generated on Thu Nov 29 16:38:52 2001 by

Functions
double	LnGamma (double z)
double	Gamma (double a, double x)
double	GamCf (double a, double x)
double	GamSer (double a, double x)

gamma.h File Reference

Functions

Detailed Description

Function Documentation