256 Chapter 6. Special Functions #include <math.h> #include "complex.h" #define EPS 6.0e-8 #define MAXIT 100 #define FPMIN 1.0e-30 #define XMIN 1.5 #define PI 3.1415927 #define PIBY2 (PI/2.0) Here EPS is the relative error:MAXIT is the maximum number of iterations allowed;FPMIN is a number near the smallest representable floating-point number:XMIN is the dividing line between using the series and continued fraction #define TRUE 1 #define ONE Complex(1.0,0.0) void frenel(float x,float *s,float *c) 83g Computes the Fresnel integrals S(x)and C(z)for all real x. granted for 19881992 void nrerror(char error_text[]); 1600 (including this one) int k,n,odd; float a,ax,fact,pix2,sign,sum,sumc,sums,term,test; fcomplex b,cc,d,h,del,cs; 7423 from NUMERICAL RECIPES IN ax=fabs(x); if (ax sqrt(FPMIN)){ Special case:avoid failure of convergence *8=0.0; test because of underflow. 米C=aX: to any server computer, else if (ax <XMIN){ Evaluate both series simultaneously. (North America sum=sums=0.0; one paper /Cambridge University Press.Programs THE ART sumc=ax; s1gn=1.0: only). fact=PIBY2*ax*ax; odd=TRUE; send copy for their term=ax; n=3; for (k=1;k<=MAXIT;k++){ email term *fact/k; sum sign*term/n; test=fabs(sum)*EPS; if (odd){ sign =-sign; sums=sum; OF SCIENTIFIC COMPUTING (ISBN 0-521 sum=sumc; else sumc=sum: sum=sums; @cambridge.org 1988-1992 by Numerical Recipes if (term test)break; -431085 odd=!odd; n+=2; if (k MAXIT)nrerror("series failed in frenel"); (outside North Software. *s=sums; *c=sumc; else Evaluate continued fraction by modified pix2-PI*ax*ax; Lentz's method ($5.2). b=Complex(1.0,-pix2); cc=Complex(1.0/FPMIN,0.0) d=h=Cdiv(ONE,b); n=-1; for (k=2;k<=MAXIT;k++){ n+=2; a=-n*(n+1); b=Cadd(b,Complex(4.0,0.0)); d=Cdiv(ONE,Cadd(RCmul(a,d),b)); Denominators cannot be zero256 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 #define MAXIT 100 #define FPMIN 1.0e-30 #define XMIN 1.5 #define PI 3.1415927 #define PIBY2 (PI/2.0) Here EPS is the relative error; MAXIT is the maximum number of iterations allowed; FPMIN is a number near the smallest representable floating-point number; XMIN is the dividing line between using the series and continued fraction. #define TRUE 1 #define ONE Complex(1.0,0.0) void frenel(float x, float *s, float *c) Computes the Fresnel integrals S(x) and C(x) for all real x. { void nrerror(char error_text[]); int k,n,odd; float a,ax,fact,pix2,sign,sum,sumc,sums,term,test; fcomplex b,cc,d,h,del,cs; ax=fabs(x); if (ax < sqrt(FPMIN)) { Special case: avoid failure of convergence *s=0.0; test because of underflow. *c=ax; } else if (ax <= XMIN) { Evaluate both series simultaneously. sum=sums=0.0; sumc=ax; sign=1.0; fact=PIBY2*ax*ax; odd=TRUE; term=ax; n=3; for (k=1;k<=MAXIT;k++) { term *= fact/k; sum += sign*term/n; test=fabs(sum)*EPS; if (odd) { sign = -sign; sums=sum; sum=sumc; } else { sumc=sum; sum=sums; } if (term < test) break; odd=!odd; n += 2; } if (k > MAXIT) nrerror("series failed in frenel"); *s=sums; *c=sumc; } else { Evaluate continued fraction by modified pix2=PI*ax*ax; Lentz’s method (§5.2). b=Complex(1.0,-pix2); cc=Complex(1.0/FPMIN,0.0); d=h=Cdiv(ONE,b); n = -1; for (k=2;k<=MAXIT;k++) { n += 2; a = -n*(n+1); b=Cadd(b,Complex(4.0,0.0)); d=Cdiv(ONE,Cadd(RCmul(a,d),b)); Denominators cannot be zero