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 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). 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.