6.1 Gamma,Beta,and Related Functions 213 Hart,J.F.,et al.1968,Computer Approximations (New York:Wiley). Hastings,C.1955,Approximations for Digita/Computers(Princeton:Princeton University Press). Luke,Y.L.1975,Mathematical Functions and Their Approximations(New York:Academic Press). 6.1 Gamma Function,Beta Function,Factorials, Binomial Coefficients ://www.nr The gamma function is defined by the integral T(z)= t=-le-tdt (6.1.1) When the argument z is an integer,the gamma function is just the familiar factorial function,but offset by one, n!=T(n+1) (6.1.2) The gamma function satisfies the recurrence relation 2为 9 T(z+1)=2T(z) (6.1.3) If the function is known for arguments z>1 or,more generally,in the half complex 0a是g%9 plane Re(z)>1 it can be obtained for z0) You can see that this is a sort of take-off on Stirling's approximation,but with a series of corrections that take into account the first few poles in the left complex plane.The constant co is very nearly equal to 1.The error term is parametrized by e. For5,N-6,and a certain set ofc's,the error is smaller than e0
6.1 Gamma, Beta, and Related Functions 213 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). Hart, J.F., et al. 1968, Computer Approximations (New York: Wiley). Hastings, C. 1955, Approximations for Digital Computers (Princeton: Princeton University Press). Luke, Y.L. 1975, Mathematical Functions and Their Approximations (New York: Academic Press). 6.1 Gamma Function, Beta Function, Factorials, Binomial Coefficients The gamma function is defined by the integral Γ(z) = ∞ 0 t z−1e−t dt (6.1.1) When the argument z is an integer, the gamma function is just the familiar factorial function, but offset by one, n! = Γ(n + 1) (6.1.2) The gamma function satisfies the recurrence relation Γ(z + 1) = zΓ(z) (6.1.3) If the function is known for arguments z > 1 or, more generally, in the half complex plane Re(z) > 1 it can be obtained for z 0) (6.1.5) You can see that this is a sort of take-off on Stirling’s approximation, but with a series of corrections that take into account the first few poles in the left complex plane. The constant c0 is very nearly equal to 1. The error term is parametrized by . For γ = 5, N = 6, and a certain set of c’s, the error is smaller than || 0.
214 Chapter 6.Special Functions It is better to implement In T(z)than I(x).since the latter will overflow many computers'floating-point representation at quite modest values of x.Often the gamma function is used in calculations where the large values ofI()are divided by other large numbers,with the result being a perfectly ordinary value.Such operations would normally be coded as subtraction of logarithms.With(6.1.5)in hand,we can compute the logarithm of the gamma function with two calls to a logarithm and 25 or so arithmetic operations.This makes it not much more difficult than other built-in functions that we take for granted,such as sinz or e: #include float gammln(float xx) Returns the value Inr(xx)for xx 0. 鱼 18881892 Internal arithmetic will be done in double precision,a nicety that you can omit if five-figure -00 g accuracy is good enough. double x,y,tmp,ser; stat1 c doub1ecof[6]=[76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5]; int ji to make y=x=xx; tmp=x+5.5; (North America 州bMe se UnN电.t tmp -(x+0.5)*log(tmp); one paper ART ser=1.000000000190015; for (j=0;j 药 Numerical Recipes 10-621 float factrl(int n) Returns the value n!as a floating-point number. 43108 float gammln(float xx); void nrerror(char error._text[☐)； (outside static int ntop=4; static f1oata[33]=[1.0,1.0,2.0,6.0,24.0}; Fill in table only as required. Software. int j; if (n O)nrerror("Negative factorial in routine factrl"); Ame if (n 32)return exp(gammln(n+1.0)); Larger value than size of table is required.Actually,this big a value is going to overflow on many computers,but no harm in trying. while (ntop<n){ Fill in table up to desired value. j=ntop++; a[ntop]=a[j]*ntop; return a[n];
214 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). It is better to implement ln Γ(x) than Γ(x), since the latter will overflow many computers’ floating-point representation at quite modest values of x. Often the gamma function is used in calculations where the large values of Γ(x) are divided by other large numbers, with the result being a perfectly ordinary value. Such operations would normally be coded as subtraction of logarithms. With (6.1.5) in hand, we can compute the logarithm of the gamma function with two calls to a logarithm and 25 or so arithmetic operations. This makes it not much more difficult than other built-in functions that we take for granted, such as sin x or ex: #include float gammln(float xx) Returns the value ln[Γ(xx)] for xx > 0. { Internal arithmetic will be done in double precision, a nicety that you can omit if five-figure accuracy is good enough. double x,y,tmp,ser; static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5}; int j; y=x=xx; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.000000000190015; for (j=0;j float factrl(int n) Returns the value n! as a floating-point number. { float gammln(float xx); void nrerror(char error_text[]); static int ntop=4; static float a[33]={1.0,1.0,2.0,6.0,24.0}; Fill in table only as required. int j; if (n 32) return exp(gammln(n+1.0)); Larger value than size of table is required. Actually, this big a value is going to overflow on many computers, but no harm in trying. while (ntop<n) { Fill in table up to desired value. j=ntop++; a[ntop]=a[j]*ntop; } return a[n]; }
6.1 Gamma,Beta,and Related Functions 215 A useful point is that factrl will be exact for the smaller values of n,since floating-point multiplies on small integers are exact on all computers.This exactness will not hold if we turn to the logarithm of the factorials.For binomial coefficients. however,we must do exactly this,since the individual factorials in a binomial coefficient will overflow long before the coefficient itself will. The binomial coefficient is defined by n n! k!(n-k)! 0≤k≤n (6.1.6) #include 1.00 float bico(int n,int k) Returns the binomial coefficient (as a floating-point number. to any float factln(int n); return floor(0.5+exp(factln(n)-factln(k)-factln(n-k))); The floor function cleans up roundoff error for smaller values of n and k. server computer, (North America make one paper University Press. THE ART which uses 是 Programs float factln(int n) strictly proh Returns In(n!). float gammln(float xx); void nrerror(char error._text[☐)： static float a[101]; A static array is automatically initialized to zero. to dir if (n<0)nrerror("Negative factorial in routine factin"); 1881892 OF SCIENTIFIC COMPUTING(ISBN if (n <1)return 0.0; if (n <100)return a[n]a[n](a[n]=gammln(n+1.0));In range of table. else return gammln(n+1.0); Out of range of table. If your problem requires a series of related binomial coefficients,a good idea Numerical Recipes 10-621 43106 is to use recurrence relations,for example (outside )=t+(=(份+() North Software. (6.1.7) n）n-kn k+1=+k Finally,turning away from the combinatorial functions with integer valued arguments,we come to the beta function, B(a,w)=B(u,2)=t-1(1-t)"-1dt (6.1.8)
6.1 Gamma, Beta, and Related Functions 215 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). A useful point is that factrl will be exact for the smaller values of n, since floating-point multiplies on small integers are exact on all computers. This exactness will not hold if we turn to the logarithm of the factorials. For binomial coefficients, however, we must do exactly this, since the individual factorials in a binomial coefficient will overflow long before the coefficient itself will. The binomial coefficient is defined by n k = n! k!(n − k)! 0 ≤ k ≤ n (6.1.6) #include float bico(int n, int k) Returns the binomial coefficient n k as a floating-point number. { float factln(int n); return floor(0.5+exp(factln(n)-factln(k)-factln(n-k))); The floor function cleans up roundoff error for smaller values of n and k. } which uses float factln(int n) Returns ln(n!). { float gammln(float xx); void nrerror(char error_text[]); static float a[101]; A static array is automatically initialized to zero. if (n < 0) nrerror("Negative factorial in routine factln"); if (n <= 1) return 0.0; if (n <= 100) return a[n] ? a[n] : (a[n]=gammln(n+1.0)); In range of table. else return gammln(n+1.0); Out of range of table. } If your problem requires a series of related binomial coefficients, a good idea is to use recurrence relations, for example n + 1 k = n + 1 n − k + 1n k = n k + n k − 1 n k + 1 = n − k k + 1 n k (6.1.7) Finally, turning away from the combinatorial functions with integer valued arguments, we come to the beta function, B(z,w) = B(w, z) = 1 0 t z−1(1 − t) w−1dt (6.1.8)
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 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 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)