216 Chapter 6.Special Functions which is related to the gamma function by B(z,w)= I(z)T(w) I(z+w) (6.1.9) hence #include <math.h> float beta(float z,float w) Returns the value of the beta function B(z,w). 83 float gammln(float xx); granted for 19881992 return exp(gammln(z)+gammln(w)-gammln(z+w)); Ca1-800-72 (including this one) /Cambridge users to make one paper NUMERICAL RECIPES IN CITED REFERENCES AND FURTHER READING: (Nort server Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by America computer e University Press. THE Dover Publications,New York),Chapter 6. ART Lanczos,C.1964,SIAM Journal on Numerical Analysis,ser.B.vol.1,pp.86-96.[1] 9 6.2 Incomplete Gamma Function,Error Function,Chi-Square Probability Function, Cumulative Poisson Function OF SCIENTIFIC COMPUTING (ISBN 198918920 The incomplete gamma function is defined by 10-521 Pa,x)=a,=1 e-ta-ldt 43106 r(a) -I(a)Jo (a>0) (6.2.1) Numerical Recipes It has the limiting values (outside North Software. P(a,0)=0andP(a,o)）=1 (6.2.2) The incomplete gamma function P(a,z)is monotonic and (for a greater than one or visit website so)rises from“near-zero”to“near-unity”in a range of x centered on about a-l, and of width about va (see Figure 6.2.1). The complement of P(a.is also confusingly called an incomplete gamma function. oa1-Pa=筒=高厂e- (a>0)(6.2.3)
216 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). which is related to the gamma function by B(z,w) = Γ(z)Γ(w) Γ(z + w) (6.1.9) hence #include <math.h> float beta(float z, float w) Returns the value of the beta function B(z, w). { float gammln(float xx); return exp(gammln(z)+gammln(w)-gammln(z+w)); } CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), Chapter 6. Lanczos, C. 1964, SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp. 86–96. [1] 6.2 Incomplete Gamma Function, Error Function, Chi-Square Probability Function, Cumulative Poisson Function The incomplete gamma function is defined by P(a, x) ≡ γ(a, x) Γ(a) ≡ 1 Γ(a) x 0 e−t t a−1dt (a > 0) (6.2.1) It has the limiting values P(a, 0) = 0 and P(a, ∞)=1 (6.2.2) The incomplete gamma function P(a, x) is monotonic and (for a greater than one or so) rises from “near-zero” to “near-unity” in a range of x centered on about a − 1, and of width about √a (see Figure 6.2.1). The complement of P(a, x) is also confusingly called an incomplete gamma function, Q(a, x) ≡ 1 − P(a, x) ≡ Γ(a, x) Γ(a) ≡ 1 Γ(a) ∞ x e−t t a−1dt (a > 0) (6.2.3)
6.2 Incomplete Gamma Function 217 TT 1.0 0.5 1.0 a=3.0 6 a=10 Permission is granted for 2 .com or call 1-800-872- (including this one) internet 0 from NUMERICAL RECIPES IN C: 0 2 68 10 12 14 to any server computer, -7423 (North America t users to make one paper 1988-1992 by Cambridge University Press. THE Figure 6.2.1.The incomplete gamma function P(a,x)for four values of a. only). It has the limiting values Programs Q(a,0)=1 and Q(a,o)=0 (6.2.4) The notations P(a,r),(a,x),and I(a,x)are standard;the notation Q(a,x)is specific to this book. There is a series development for (a,)as follows: 1881892 ART OF SCIENTIFIC COMPUTING (ISBN Y(a,r)=e-rza I(a) (6.2.5) n=0 T(a+1+n) 10-521 One does not actually need to compute a new I(a+1+n)for each n;one rather uses equation (6.1.3)and the previous coefficient. Numerical Recipes 431085 A continued fraction development for r(a,x)is (outside T(a,r)=e-=ra 11-a12-a2 x+1+x+1+x+ (x>0) (6.2.6) North Software. It is computationally better to use the even part of(6.2.6),which converges twice as fast (see 85.2): visit website I(a,z)=e-=za 11·(1-a)2.(2-a) x+1-a-x+3-a-x+5-a- (x>0) (6.2.7) It turns out that (6.2.5)converges rapidly for x less than about a +1,while (6.2.6)or(6.2.7)converges rapidly for x greater than about a+1.In these respective
6.2 Incomplete Gamma Function 217 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 0 2 4 6 8 10 12 14 0 .2 .4 .6 .8 1.0 a = 3.0 1.0 0.5 incomplete gamma function P(a,x) x a = 10 Figure 6.2.1. The incomplete gamma function P(a, x) for four values of a. It has the limiting values Q(a, 0) = 1 and Q(a, ∞)=0 (6.2.4) The notations P(a, x), γ(a, x), and Γ(a, x) are standard; the notation Q(a, x) is specific to this book. There is a series development for γ(a, x) as follows: γ(a, x) = e−xxa ∞ n=0 Γ(a) Γ(a +1+ n) xn (6.2.5) One does not actually need to compute a new Γ(a +1+ n) for each n; one rather uses equation (6.1.3) and the previous coefficient. A continued fraction development for Γ(a, x) is Γ(a, x) = e−xxa 1 x + 1 − a 1 + 1 x + 2 − a 1 + 2 x + ··· (x > 0) (6.2.6) It is computationally better to use the even part of (6.2.6), which converges twice as fast (see §5.2): Γ(a, x) = e−xxa 1 x + 1 − a − 1 · (1 − a) x + 3 − a − 2 · (2 − a) x + 5 − a − ··· (x > 0) (6.2.7) It turns out that (6.2.5) converges rapidly for x less than about a + 1, while (6.2.6) or (6.2.7) converges rapidly for x greater than about a+ 1. In these respective
218 Chapter 6.Special Functions regimes each requires at most a few times va terms to converge,and this many only near =a,where the incomplete gamma functions are varying most rapidly Thus (6.2.5)and (6.2.7)together allow evaluation of the function for all positive a and x.An extra dividend is that we never need compute a function value near zero by subtracting two nearly equal numbers.The higher-level functions that return P(a,r)and Q(a,x)are float gammp(float a,float x) Returns the incomplete gamma function P(a,z). f void gcf(float *gammcf,float a,float x,float *gln); void gser(float *gamser,float a,float x,float *gln); void nrerror(char error_text []) (including float gamser,gammcf,gIn; 19881992 if (x <0.0 1I a <=0.0)nrerror("Invalid arguments in routine gammp"); 11800 if(x<(a+1.0)){ Use the series representation. gser(&gamser,a,x,&gln); return gamser; n NUMERICAL RECIPES else Use the continued fraction representation gcf (&gammcf,a,x,&gln); return 1.0-gammcf; and take its complement. server computer, (North America to make one paper University Press. THE ART float gammq(float a,float x) Returns the incomplete gamma function Q(a,x)=1-P(a,z). Programs void gcf(float *gammcf,float a,float x,float *gln); strictly proh void gser(float *gamser,float a,float x,float *gln); void nrerror(char error_text[]); float gamser,gammcf,gln; to dir if (x 0.0 II a <0.0)nrerror("Invalid arguments in routine gammq"); if(x<(a+1.0)){ Use the series representation gser(&gamser,a,x,&gln); return 1.0-gamser; and take its complement. 1881992 OF SCIENTIFIC COMPUTING(ISBN else Use the continued fraction representation. gcf (&gammcf,a,x,&gln); v@cam return gammcf; Numerical Recipes 021 43108 The argument gin is set by both the series and continued fraction procedures (outside to the value In I(a);the reason for this is so that it is available to you if you want to modify the above two procedures to give (a,z)and I(a,x),in addition to P(a,x) North Software. and (a,x)(cf.equations 6.2.1 and 6.2.3). The functions gser and gcf which implement(6.2.5)and(6.2.7)are #include <math.h> #define ITMAX 100 #define EPS 3.0e-7 void gser(float *gamser,float a,float x,float *gln) Returns the incomplete gamma function P(a,evaluated by its series representation as gamser Also returns InI(a)as gln. float gammln(float xx);
218 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). regimes each requires at most a few times √a terms to converge, and this many only near x = a, where the incomplete gamma functions are varying most rapidly. Thus (6.2.5) and (6.2.7) together allow evaluation of the function for all positive a and x. An extra dividend is that we never need compute a function value near zero by subtracting two nearly equal numbers. The higher-level functions that return P(a, x) and Q(a, x) are float gammp(float a, float x) Returns the incomplete gamma function P(a, x). { void gcf(float *gammcf, float a, float x, float *gln); void gser(float *gamser, float a, float x, float *gln); void nrerror(char error_text[]); float gamser,gammcf,gln; if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp"); if (x < (a+1.0)) { Use the series representation. gser(&gamser,a,x,&gln); return gamser; } else { Use the continued fraction representation gcf(&gammcf,a,x,&gln); return 1.0-gammcf; and take its complement. } } float gammq(float a, float x) Returns the incomplete gamma function Q(a, x) ≡ 1 − P(a, x). { void gcf(float *gammcf, float a, float x, float *gln); void gser(float *gamser, float a, float x, float *gln); void nrerror(char error_text[]); float gamser,gammcf,gln; if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq"); if (x < (a+1.0)) { Use the series representation gser(&gamser,a,x,&gln); return 1.0-gamser; and take its complement. } else { Use the continued fraction representation. gcf(&gammcf,a,x,&gln); return gammcf; } } The argument gln is set by both the series and continued fraction procedures to the value ln Γ(a); the reason for this is so that it is available to you if you want to modify the above two procedures to give γ(a, x) and Γ(a, x), in addition to P(a, x) and Q(a, x) (cf. equations 6.2.1 and 6.2.3). The functions gser and gcf which implement (6.2.5) and (6.2.7) are #include <math.h> #define ITMAX 100 #define EPS 3.0e-7 void gser(float *gamser, float a, float x, float *gln) Returns the incomplete gamma function P(a, x) evaluated by its series representation as gamser. Also returns ln Γ(a) as gln. { float gammln(float xx);
6.2 Incomplete Gamma Function 219 void nrerror(char error_text []) int n; float sum,del,api *gln=gammln(a); if(x<=0.0){ 1f(x<0.0) nrerror("x less than 0 in routine gser"); *gamser=0.0; return; else ap=a; del=sum=1.0/a; http://www.nr. read able files Permission is for (n=1;n<=ITMAX;n++){ ++ap; del *x/api sum del: if (fabs(del)<fabs(sum)EPS){ *gamser=sum*exp(-x+a*log(x)-(*gln)) return; .com or call 1-800-872- (including this one) granted fori internet nrerror("a too large,ITMAX too small in routine gser"); return: 7423 (North America to any server computer,is t users to make one paper 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN C: THE #include <math.h> 是 #define ITMAX 100 Maximum allowed number of iterations Programs #define EPS 3.0e-7 Relative accuracy. #define FPMIN 1.0e-30 Number near the smallest representable floating-point number void gcf(float *gammcf,float a,float x, float *gln) Returns the incomplete gamma function Q(a,x)evaluated by its continued fraction represen- to dir Copyright (C) tation as gammcf.Also returns InI(a)as gln. float gammln(float xx); void nrerror(char error_text []) int i; ART OF SCIENTIFIC COMPUTING(ISBN 0-521 float an,b,c,d,del,h; *gln=gammln(a); b=x+1.0-a; Set up for evaluating continued fraction c=1.0/FPMIN; by modified Lentz's method ($5.2) 餐8质 1988-1992 by Numerical Recipes d=1.0/b; with bo =0. Numerical Recipes books or -431085 h=d; for (i=1;i<=ITMAX;i++) Iterate to convergence an =-i(i-a); b+=2.0; d=an*d+b; if (fabs(d)<FPMIN)d=FPMIN: c=b+an/c; (outside North America) Software. if (fabs(c)<FPMIN)c=FPMIN; d=1.0/d; del=d*c; h=del; if (fabs(del-1.0)<EPS)break; if (i>ITMAX)nrerror("a too large,ITMAX too small in gcf"); *gammcf=exp(-x+a*log(x)-(*gln))*h; Put factors in front
6.2 Incomplete Gamma Function 219 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). void nrerror(char error_text[]); int n; float sum,del,ap; *gln=gammln(a); if (x <= 0.0) { if (x < 0.0) nrerror("x less than 0 in routine gser"); *gamser=0.0; return; } else { ap=a; del=sum=1.0/a; for (n=1;n<=ITMAX;n++) { ++ap; del *= x/ap; sum += del; if (fabs(del) < fabs(sum)*EPS) { *gamser=sum*exp(-x+a*log(x)-(*gln)); return; } } nrerror("a too large, ITMAX too small in routine gser"); return; } } #include <math.h> #define ITMAX 100 Maximum allowed number of iterations. #define EPS 3.0e-7 Relative accuracy. #define FPMIN 1.0e-30 Number near the smallest representable floating-point number. void gcf(float *gammcf, float a, float x, float *gln) Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation as gammcf. Also returns ln Γ(a) as gln. { float gammln(float xx); void nrerror(char error_text[]); int i; float an,b,c,d,del,h; *gln=gammln(a); b=x+1.0-a; Set up for evaluating continued fraction by modified Lentz’s method (§5.2) with b0 = 0. c=1.0/FPMIN; d=1.0/b; h=d; for (i=1;i<=ITMAX;i++) { Iterate to convergence. an = -i*(i-a); b += 2.0; d=an*d+b; if (fabs(d) < FPMIN) d=FPMIN; c=b+an/c; if (fabs(c) < FPMIN) c=FPMIN; d=1.0/d; del=d*c; h *= del; if (fabs(del-1.0) < EPS) break; } if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf"); *gammcf=exp(-x+a*log(x)-(*gln))*h; Put factors in front. }
220 Chapter 6. Special Functions Error Function The error function and complementary error function are special cases of the incomplete gamma function,and are obtained moderately efficiently by the above procedures.Their definitions are 2 erf(x)= e-dt (6.2.8) and 2 erfc(x)≡l-erf(x)= e-dt (6.2.9) The functions have the following limiting values and symmetries: erf(0)=0 erf(oo)=1 erf(-x)=-erf(x) (6.2.10) erfc(0)=1 erfc(oo)=0 erfc(-x)=2-erfc(x) (6.2.11) 令 They are related to the incomplete gamma functions by Press. e=P(（2） (x≥0) (6.2.12) ART and 9 Program erfc(z）= (x≥0) (6.2.13) We'll put an extra"f into our routine names to avoid conflicts with names already to dir in some C libraries: float erff(float x) Returns the error function erf(x). float gammp(float a,float x); return x 0.0 -gammp(0.5,x*x):gammp(0.5,x*x); @cambridge.org(outside North America). Numerical Recipes books or 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING (ISBN 10-621 431085 float erffc(float x) Returns the complementary error function erfc(x). Software. float gammp(float a,float x); float gammg(float a,float x); ying of machine return x 0.0 1.0+gammp(0.5,x*x):gammq(0.5,x*x); 2 If you care to do so,you can easily remedy the minor inefficiency in erff and erffc,namely that I(0.5)=Vr is computed unnecessarily when gammp or gammq is called.Before you do that,however,you might wish to consider the following routine,based on Chebyshev fitting to an inspired guess as to the functional form:
220 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Error Function The error function and complementary error function are special cases of the incomplete gamma function, and are obtained moderately efficiently by the above procedures. Their definitions are erf(x) = 2 √π x 0 e−t2 dt (6.2.8) and erfc(x) ≡ 1 − erf(x) = 2 √π ∞ x e−t2 dt (6.2.9) The functions have the following limiting values and symmetries: erf(0) = 0 erf(∞)=1 erf(−x) = −erf(x) (6.2.10) erfc(0) = 1 erfc(∞)=0 erfc(−x)=2 − erfc(x) (6.2.11) They are related to the incomplete gamma functions by erf(x) = P 1 2 , x2 (x ≥ 0) (6.2.12) and erfc(x) = Q 1 2 , x2 (x ≥ 0) (6.2.13) We’ll put an extra “f” into our routine names to avoid conflicts with names already in some C libraries: float erff(float x) Returns the error function erf(x). { float gammp(float a, float x); return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x); } float erffc(float x) Returns the complementary error function erfc(x). { float gammp(float a, float x); float gammq(float a, float x); return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x); } If you care to do so, you can easily remedy the minor inefficiency in erff and erffc, namely that Γ(0.5) = √π is computed unnecessarily when gammp or gammq is called. Before you do that, however, you might wish to consider the following routine, based on Chebyshev fitting to an inspired guess as to the functional form: