222 Chapter 6.Special Functions CITED REFERENCES AND FURTHER READING: 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 Dover Publications,New York),Chapters 6,7,and 26. Pearson,K.(ed.)1951,Tables of the Incomplete Gamma Function (Cambridge:Cambridge University Press). 6.3 Exponential Integrals The standard definition of the exponential integral is 鱼 En()= d tn x>0,n=0,1,. (6.3.1) The function defined by the principal value of the integral ⊙ 令 =-= x>0 (6.3.2) Press. is also called an exponential integral.Note that Ei(-)is related to -E1(x)by analytic continuation. 9 Program The function En(x)is a special case of the incomplete gamma function En(x)=r"-1I(1-n,x) (6.3.3) OF SCIENTIFIC( We can therefore use a similar strategy for evaluating it.The continued fraction- just equation(6.2.6)rewritten-converges for all z >0: En(x)=e-x 1n1n+12 COMPUTING (ISBN x+1+x+1++… (6.3.4) 1920 Numerical 10-521 We use it in its more rapidly converging even form, 43108 1 1.n 2(n+1) 、 En(x）=e- (6.3.5) (outside Recipes x+n-x+n+2-x+n+4- The continued fraction only really converges fast enough to be useful for z1. Software. For 0<z1,we can use the series representation E.(回= a2-lnx+(n]、 (-x)m (6.3.6) m-n+1)m! 020 The quantity (n)here is the digamma function,given for integer arguments by -1 (1)=-7, =+ (6.3.7)
222 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). 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), Chapters 6, 7, and 26. Pearson, K. (ed.) 1951, Tables of the Incomplete Gamma Function (Cambridge: Cambridge University Press). 6.3 Exponential Integrals The standard definition of the exponential integral is En(x) = ∞ 1 e−xt tn dt, x > 0, n = 0, 1,... (6.3.1) The function defined by the principal value of the integral Ei(x) = − ∞ −x e−t t dt = x −∞ et t dt, x > 0 (6.3.2) is also called an exponential integral. Note that Ei(−x) is related to −E1(x) by analytic continuation. The function En(x) is a special case of the incomplete gamma function En(x) = xn−1Γ(1 − n, x) (6.3.3) We can therefore use a similar strategy for evaluating it. The continued fraction — just equation (6.2.6) rewritten — converges for all x > 0: En(x) = e−x 1 x + n 1 + 1 x + n + 1 1 + 2 x + ··· (6.3.4) We use it in its more rapidly converging even form, En(x) = e−x 1 x + n − 1 · n x + n + 2 − 2(n + 1) x + n + 4 − ··· (6.3.5) The continued fraction only really converges fast enough to be useful for x >∼ 1. For 0 < x <∼ 1, we can use the series representation En(x) = (−x)n−1 (n − 1)! [− ln x + ψ(n)] − ∞ m=0 m=n−1 (−x)m (m − n + 1)m! (6.3.6) The quantity ψ(n) here is the digamma function, given for integer arguments by ψ(1) = −γ, ψ(n) = −γ + n −1 m=1 1 m (6.3.7)
6.3 Exponential Integrals 223 where y=0.5772156649...is Euler's constant.We evaluate the expression(6.3.6) in order of ascending powers of x: 1 En(x)=- (-x)n-2 1-n) -2-m)1+-n1②-…+ (-1)(n-21 a-血x+j--y+-n+ 1·n +2n++ (6.3.8) The first square bracket is omitted when n=1.This method of evaluation has the 83g advantage that for large n the series converges before reaching the term containing (n).Accordingly,one needs an algorithm for evaluating (n)only for small n, n s20-40.We use equation (6.3.7),although a table look-up would improve efficiency slightly. Amos [1]presents a careful discussion of the truncation error in evaluating equation(6.3.8),and gives a fairly elaborate termination criterion.We have found that simply stopping when the last term added is smaller than the required tolerance (Nort server works about as well. Two special cases have to be handled separately: Americ computer, University Press. THE ART e-z Eo(x)= (6.3.9) 9 Program En(0)= n-1 n>1 The routine expint allows fast evaluation of En()to any accuracy EPS within the reach of your machine's word length for floating-point numbers. The to dir only modification required for increased accuracy is to supply Euler's constant with enough significant digits.Wrench [2]can provide you with the first 328 digits if OF SCIENTIFIC COMPUTING(ISBN 1988-19920 necessary! #include 10-521 #define MAXIT 100 Maximum allowed number of iterations. #define EULER 0.5772156649 Euler's constant #define FPMIN 1.0e-30 Close to smallest representable floating-point number. Numerical Recipes 43108 #define EPS 1.0e-7 Desired relative error,not smaller than the machine pre- cision. float expint(int n,float x) (outside 膜 Evaluates the exponential integral En(z). Software. void nrerror(char error_text[]); int i,ii,nm1; float a,b,c,d,del,fact,h,psi,ansi nm1=n-1; if(n<0Ilx<0.0Il(x==0.02数(n=0I川n==1))) nrerror("bad arguments in expint"); else if (n ==0)ans=exp(-x)/x; Special case. else if(x==0.0)ans=1.0/nm1; Another special case. else
6.3 Exponential Integrals 223 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). where γ = 0.5772156649 ... is Euler’s constant. We evaluate the expression (6.3.6) in order of ascending powers of x: En(x) = − 1 (1 − n) − x (2 − n) · 1 + x2 (3 − n)(1 · 2) −··· + (−x)n−2 (−1)(n − 2)! + (−x)n−1 (n − 1)! [− ln x + ψ(n)] − (−x)n 1 · n! + (−x)n+1 2 · (n + 1)! + ··· (6.3.8) The first square bracket is omitted when n = 1. This method of evaluation has the advantage that for large n the series converges before reaching the term containing ψ(n). Accordingly, one needs an algorithm for evaluating ψ(n) only for small n, n 1 (6.3.9) The routine expint allows fast evaluation of En(x) to any accuracy EPS within the reach of your machine’s word length for floating-point numbers. The only modification required for increased accuracy is to supply Euler’s constant with enough significant digits. Wrench [2] can provide you with the first 328 digits if necessary! #include #define MAXIT 100 Maximum allowed number of iterations. #define EULER 0.5772156649 Euler’s constant γ. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define EPS 1.0e-7 Desired relative error, not smaller than the machine precision. float expint(int n, float x) Evaluates the exponential integral En(x). { void nrerror(char error_text[]); int i,ii,nm1; float a,b,c,d,del,fact,h,psi,ans; nm1=n-1; if (n < 0 || x < 0.0 || (x==0.0 && (n==0 || n==1))) nrerror("bad arguments in expint"); else { if (n == 0) ans=exp(-x)/x; Special case. else { if (x == 0.0) ans=1.0/nm1; Another special case. else {
224 Chapter 6.Special Functions 1f(x>1.0)[ Lentz's algorithm ($5.2). bex+n; c=1.0/FPMIN; d=1.0/b: hed: for (i=1;i<=MAXIT;i++) a=-i*(nm1+i); b+=2.0; d=1.0/(a*d+b); Denominators cannot be zero. c=b+a/c; del=c*d; h*=del: if (fabs(del-1.0)<EPS){ ans=h*exp(-x); return ans; 83 鱼 granted for 19881992 nrerror("continued fraction failed in expint"); 1.800 else Evaluate series. ans (nm1!=0 1.0/nm1 -log(x)-EULER); Set first term. fact=1.0; Cambridge for (i=1;i<=MAXIT;i++) from NUMERICAL RECIPES IN fact -x/i; if (i !nm1)del =-fact/(i-nm1); else psi =-EULER; Compute (n). for(1i=1;1i<=nm1;ii+)psi+=1.0/11; server computer, (North America UnN电.t THE del=fact*(-log(x)+psi); 是 ART ans +del; if (fabs(del)<fabs(ans)*EPS)return ans; Programs nrerror("series failed in expint"); strictly prohibited. nail to dir return ans; 2 198819920 OF SCIENTIFIC COMPUTING(ISBN A good algorithm for evaluating Ei is to use the power series for small x and the asymptotic series for large z.The power series is Numerical Recipes 10-6211 43108 Ei()=y+lnx+1.五+2：2+ (6.3.10) (outside where is Euler's constant.The asymptotic expansion is North Software. Ei(z) +++ 1 + (6.3.11) The lower limit for the use of the asymptotic expansion is approximately |In EPS, where EPS is the required relative error
224 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). if (x > 1.0) { Lentz’s algorithm (§5.2). b=x+n; c=1.0/FPMIN; d=1.0/b; h=d; for (i=1;i<=MAXIT;i++) { a = -i*(nm1+i); b += 2.0; d=1.0/(a*d+b); Denominators cannot be zero. c=b+a/c; del=c*d; h *= del; if (fabs(del-1.0) < EPS) { ans=h*exp(-x); return ans; } } nrerror("continued fraction failed in expint"); } else { Evaluate series. ans = (nm1!=0 ? 1.0/nm1 : -log(x)-EULER); Set first term. fact=1.0; for (i=1;i<=MAXIT;i++) { fact *= -x/i; if (i != nm1) del = -fact/(i-nm1); else { psi = -EULER; Compute ψ(n). for (ii=1;ii<=nm1;ii++) psi += 1.0/ii; del=fact*(-log(x)+psi); } ans += del; if (fabs(del) < fabs(ans)*EPS) return ans; } nrerror("series failed in expint"); } } } } return ans; } A good algorithm for evaluating Ei is to use the power series for small x and the asymptotic series for large x. The power series is Ei(x) = γ + ln x + x 1 · 1! + x2 2 · 2! + ··· (6.3.10) where γ is Euler’s constant. The asymptotic expansion is Ei(x) ∼ ex x 1 + 1! x + 2! x2 + ··· (6.3.11) The lower limit for the use of the asymptotic expansion is approximately | ln EPS|, where EPS is the required relative error.
6.3 Exponential Integrals 225 #include #define EULER 0.57721566 Euler's constant y. #define MAXIT 100 Maximum number of iterations allowed. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define EPS 6.0e-8 Relative error,or absolute error near the zero of Ei at x=0.3725. float ei(float x) Computes the exponential integral Ei(r)for x >0. void nrerror(char error._text[☐)； int k; float fact,prev,sum,term; if (x <0.0)nrerror("Bad argument in ei"); if (x FPMIN)return log(x)+EULER; Special case:avoid failure of convergence if (x <=-log(EPS)){ test because of underflow. 19881992 gum=0.0; Use power series. fact=1.0; 1-800 for (k=1;k<=MAXIT;k++){ fact *x/k; term=fact/k; to any n NUMERICAL RECIPES sum +term; if (term EPS*sum)break; if (k MAXIT)nrerror("Series failed in ei"); (Nort return sum+log(x)+EULER; else Use asymptotic series. America server computer, 令 Press. THE sum=0.0; Start with second term. ART term=1.0; for (k=1;k<=MAXIT;k++){ prev=term; Programs term *k/x; if (term EPS)break; Since final sum is greater than one,term itself approximates the relative error. SCIENTIFIC if (term prev)sum +term; Still converging:add new term. else sum -prev; Diverging:subtract previous term and break; exit. 1920 COMPUTING(ISBN return exp(x)*(1.0+sum)/x; Numerical Recipes 43108 CITED REFERENCES AND FURTHER READING: Stegun,I.A.,and Zucker,R.1974,Journal of Research of the National Bureau of Standards. (outside 膜 vol.78B.pp.199-216;1976,op.ct,vol.80B.pp.291-311. North oftware. Amos D.E.1980,ACM Transactions on Mathematical Software,vol.6,pp.365-377 [1]:also vol.6,Pp.420-428. 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 Dover Publications,New York),Chapter 5. Wrench J.W.1952.Mathematical Tables and Other Aids to Computation,vol.6,p.255.[2]
6.3 Exponential Integrals 225 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). #include #define EULER 0.57721566 Euler’s constant γ. #define MAXIT 100 Maximum number of iterations allowed. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define EPS 6.0e-8 Relative error, or absolute error near the zero of Ei at x = 0.3725. float ei(float x) Computes the exponential integral Ei(x) for x > 0. { void nrerror(char error_text[]); int k; float fact,prev,sum,term; if (x MAXIT) nrerror("Series failed in ei"); return sum+log(x)+EULER; } else { Use asymptotic series. sum=0.0; Start with second term. term=1.0; for (k=1;k<=MAXIT;k++) { prev=term; term *= k/x; if (term < EPS) break; Since final sum is greater than one, term itself approximates the relative error. if (term < prev) sum += term; Still converging: add new term. else { sum -= prev; Diverging: subtract previous term and break; exit. } } return exp(x)*(1.0+sum)/x; } } CITED REFERENCES AND FURTHER READING: Stegun, I.A., and Zucker, R. 1974, Journal of Research of the National Bureau of Standards, vol. 78B, pp. 199–216; 1976, op. cit., vol. 80B, pp. 291–311. Amos D.E. 1980, ACM Transactions on Mathematical Software, vol. 6, pp. 365–377 [1]; also vol. 6, pp. 420–428. 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 5. Wrench J.W. 1952, Mathematical Tables and Other Aids to Computation, vol. 6, p. 255. [2]
226 Chapter 6.Special Functions TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT (0.5,5.0) (8.010.0) /(1.0,3.0) Permission is 6 (0.5,0.5) .com or call 1-800-872- granted for (including this one) internet to any server computer, t users to make one paper 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN C: (5.0,0.5) 0 (North America THE 0 .6 Programs Figure 6.4.1.The incomplete beta function (,b)for five different pairs of (a,b).Notice that the pairs (0.5,5.0)and(5.0,0.5)are symmetrically related as indicated in equation (6.4.3). 6.4 Incomplete Beta Function,Student's Distribution,F-Distribution,Cumulative 1788-1982 ART OF SCIENTIFIC COMPUTING (ISBN 0-521 Binomial Distribution The incomplete beta function is defined by Numerical Recipes -43108-5 I(a,b)= Ba,0=1 B(a,b) B(a,b)Jo ta-1(1-t)b-1dt (a,b>0) (6.4.1) (outside Software. It has the limiting values North Io(a,b)=0 11(a,b)=1 (6.4.2) visit website f machine and the symmetry relation Iz(a,b)=1-11-x(b,a) (6.4.3) If a and b are both rather greater than one,then I(a,b)rises from"near-zero"to "near-unity"quite sharply at about x =a/(a+b).Figure 6.4.1 plots the function for several pairs (a,b)
226 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). 0 (5.0,0.5) (0.5,0.5) (8.0,10.0) (1.0,3.0) (0.5,5.0) .2 .4 .6 1 .8 0 .2 .4 .6 .8 1 incomplete beta function Ix(a,b) x Figure 6.4.1. The incomplete beta function Ix(a, b) for five different pairs of (a, b). Notice that the pairs (0.5, 5.0) and (5.0, 0.5) are symmetrically related as indicated in equation (6.4.3). 6.4 Incomplete Beta Function, Student’s Distribution, F-Distribution, Cumulative Binomial Distribution The incomplete beta function is defined by Ix(a, b) ≡ Bx(a, b) B(a, b) ≡ 1 B(a, b) x 0 t a−1(1 − t) b−1dt (a, b > 0) (6.4.1) It has the limiting values I0(a, b)=0 I1(a, b)=1 (6.4.2) and the symmetry relation Ix(a, b)=1 − I1−x(b, a) (6.4.3) If a and b are both rather greater than one, then I x(a, b) rises from “near-zero” to “near-unity” quite sharply at about x = a/(a + b). Figure 6.4.1 plots the function for several pairs (a, b)
