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 <math.h> 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. else6.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 machine￾readable 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 <∼ 20 – 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 works about as well. Two special cases have to be handled separately: E0(x) = e−x x En(0) = 1 n − 1 , 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 <math.h> #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 pre￾cision. 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 {