6.3 Exponential Integrals 225 #include <math.h> #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 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). #include <math.h> #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 <= 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. sum=0.0; Use power series. fact=1.0; for (k=1;k<=MAXIT;k++) { fact *= x/k; term=fact/k; sum += term; if (term < EPS*sum) break; } if (k > 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 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]