258 Chapter 6. Special Functions #include <math.h> #include "complex.h" #define EPS 6.0e-8 Relative error,or absolute error near a zero of Ci(). #define EULER 0.57721566 Euler's constant y. #define MAXIT 100 Maximum number of iterations allowed. #define PIBY2 1.5707963 π/2. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define TMIN 2.0 Dividing line between using the series and continued frac- #define TRUE 1 tion. #define ONE Complex(1.0,0.0) void cisi(float x,float *ci,float *si) Computes the cosine and sine integrals Ci(z)and Si(x).Ci(0)is returned as a large negative number and no error message is generated.For x<0 the routine returns Ci(-x)and you must supply the -i yourself. f (including this one) granted for 18881992 void nrerror(char error_text []) int i,k,odd; 100 float a,err,fact,sign,sum,sumc,sums,t,term fcomplex h,b,c,d,del; /Cambridge t=fabs(x); from NUMERICAL RECIPES IN 1f(t=0.0) Special case. *g1=0.0; to make *c1=-1.0/FPMIN; return; server computer, (North America UnN电.t THE if (t>TMIN){ Evaluate continued fraction by modified one paper ART b=Complex(1.0,t); Lentz's method ($5.2). c=Complex(1.0/FPMIN,0.0); 是 d=h=Cdiv(ONE,b); Programs for (i=2;i<=MAXIT;i++){ a=-(1-1)*(1-1): b=Cadd(b,Complex(2.0,0.0)); st st d=Cdiv(ONE,Cadd(RCmul(a,d),b)); Denominators cannot be zero. to dir Copyright (C) c=Cadd(b,Cdiv(Complex(a,0.0),c)) del-Cmul(c,d): h=Cmul(h,del); if (fabs(del.r-1.0)+fabs(del.i)<EPS)break; ectcustser OF SCIENTIFIC COMPUTING(ISBN if (i>MAXIT)nrerror("cf failed in cisi"); h=Cmul(Complex(cos(t),-sin(t)),h); tci =-h.r; 10-621 米si=PIBY2+h.1: else Evaluate both series simultaneously. 、海 if (t sqrt(FPMIN)){ Special case:avoid failure of convergence 1988-1992 by Numerical Recipes -43108 sumc=0.0; test because of underflow. sumset e1se (outside sum=sums=sumc=0.0; sign=fact=1.0; North Software. odd=TRUE; for (k=1;k<=MAXIT;k++){ Amer fact *=t/k; term=fact/k; sum +sign*term; err=term/fabs(sum) if (odd){ s1g即=-s1g即： sums=sum; sum=sumc else sumc=sum; sum=sums258 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). #include <math.h> #include "complex.h" #define EPS 6.0e-8 Relative error, or absolute error near a zero of Ci(x). #define EULER 0.57721566 Euler’s constant γ. #define MAXIT 100 Maximum number of iterations allowed. #define PIBY2 1.5707963 π/2. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define TMIN 2.0 Dividing line between using the series and continued frac- #define TRUE 1 tion. #define ONE Complex(1.0,0.0) void cisi(float x, float *ci, float *si) Computes the cosine and sine integrals Ci(x) and Si(x). Ci(0) is returned as a large negative number and no error message is generated. For x < 0 the routine returns Ci(−x) and you must supply the −iπ yourself. { void nrerror(char error_text[]); int i,k,odd; float a,err,fact,sign,sum,sumc,sums,t,term; fcomplex h,b,c,d,del; t=fabs(x); if (t == 0.0) { Special case. *si=0.0; *ci = -1.0/FPMIN; return; } if (t > TMIN) { Evaluate continued fraction by modified b=Complex(1.0,t); Lentz’s method (§5.2). c=Complex(1.0/FPMIN,0.0); d=h=Cdiv(ONE,b); for (i=2;i<=MAXIT;i++) { a = -(i-1)*(i-1); b=Cadd(b,Complex(2.0,0.0)); d=Cdiv(ONE,Cadd(RCmul(a,d),b)); Denominators cannot be zero. c=Cadd(b,Cdiv(Complex(a,0.0),c)); del=Cmul(c,d); h=Cmul(h,del); if (fabs(del.r-1.0)+fabs(del.i) < EPS) break; } if (i > MAXIT) nrerror("cf failed in cisi"); h=Cmul(Complex(cos(t),-sin(t)),h); *ci = -h.r; *si=PIBY2+h.i; } else { Evaluate both series simultaneously. if (t < sqrt(FPMIN)) { Special case: avoid failure of convergence sumc=0.0; test because of underflow. sums=t; } else { sum=sums=sumc=0.0; sign=fact=1.0; odd=TRUE; for (k=1;k<=MAXIT;k++) { fact *= t/k; term=fact/k; sum += sign*term; err=term/fabs(sum); if (odd) { sign = -sign; sums=sum; sum=sumc; } else { sumc=sum; sum=sums;