264 Chapter 6. Special Functions #include <math.h> #include "nrutil.h" #define ERRTOL 0.08 #define TINY 1.5e-38 #define BIG 3.0e37 #define THIRD (1.0/3.0) #define C1(1.0/24.0) #define C2 0.1 #define C3(3.0/44.0) #define C4(1.0/14.0) float rf(float x,float y,float z) Computes Carlson's elliptic integral of the first kind,RF(r,y,z).y,and z must be nonneg- ative,and at most one can be zero.TINY must be at least 5 times the machine underflow limit, BIG at most one fifth the machine overflow limit. 83g r granted for float alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt; 1.800 (including this one) 19881992 if (FMIN(FMIN(x,y),z)<0.0 II FMIN(FMIN(x+y,x+z),y+z)<TINY II FMAX(FMAX(x,y),2)>BIG) 872 nrerror("invalid arguments in rf"); Cambridge xt=x; yt=y; 乙t=Z： do sqrtx=sqrt(xt); (North America to any server computer, sqrty=sgrt(yt); to make one paper UnN电.t THE sqrtz=sgrt(zt); ART alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); copyfor their ave=THIRD*(xt+yt+zt); delx=(ave-xt)/ave; st st Programs dely=(ave-yt)/ave; delz=(ave-zt)/ave; to din while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz))>ERRTOL); e2-delx*dely-delz+delz; e3=delx*dely*delz; OF SCIENTIFIC COMPUTING(ISBN return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave); ectcustser 18881920 v@cam A value of 0.08 for the error tolerance parameter is adequate for single precision (7 Numerical Recipes 10-621 43108 significant digits).Since the error scales as en,we see that 0.0025 will yield double precision (16 significant digits)and require at most two or three more iterations.Since the coefficients of the sixth-order truncation error are different for the other elliptic functions,these values for (outside 膜 the error tolerance should be changed to 0.04 and 0.0012 in the algorithm for Rc,and 0.05 and 0.0015 for Ry and Rp.As well as being an algorithm in its own right for certain combinations North Software. of elementary functions,the algorithm for Rc is used repeatedly in the computation of Ry. The c implementations test the input arguments against two machine-dependent con- stants,TINY and BIG,to ensure that there will be no underflow or overflow during the computation.We have chosen conservative values,corresponding to a machine minimum of 3 x 10-39 and a machine maximum of 1.7 x 1038.You can always extend the range of admissible argument values by using the homogeneity relations (6.11.22),below. #include <math.h> #include "nrutil.h" #define ERRTOL 0.05264 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 "nrutil.h" #define ERRTOL 0.08 #define TINY 1.5e-38 #define BIG 3.0e37 #define THIRD (1.0/3.0) #define C1 (1.0/24.0) #define C2 0.1 #define C3 (3.0/44.0) #define C4 (1.0/14.0) float rf(float x, float y, float z) Computes Carlson’s elliptic integral of the first kind, RF (x, y, z). x, y, and z must be nonneg￾ative, and at most one can be zero. TINY must be at least 5 times the machine underflow limit, BIG at most one fifth the machine overflow limit. { float alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt; if (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(x+y,x+z),y+z) < TINY || FMAX(FMAX(x,y),z) > BIG) nrerror("invalid arguments in rf"); xt=x; yt=y; zt=z; do { sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); ave=THIRD*(xt+yt+zt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; } while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz)) > ERRTOL); e2=delx*dely-delz*delz; e3=delx*dely*delz; return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave); } A value of 0.08 for the error tolerance parameter is adequate for single precision (7 significant digits). Since the error scales as  6 n, we see that 0.0025 will yield double precision (16 significant digits) and require at most two or three more iterations. Since the coefficients of the sixth-order truncation error are different for the other elliptic functions, these values for the error tolerance should be changed to 0.04 and 0.0012 in the algorithm for RC , and 0.05 and 0.0015 for RJ and RD. As well as being an algorithm in its own right for certain combinations of elementary functions, the algorithm for RC is used repeatedly in the computation of RJ . The C implementations test the input arguments against two machine-dependent con￾stants, TINY and BIG, to ensure that there will be no underflow or overflow during the computation. We have chosen conservative values, corresponding to a machine minimum of 3 × 10−39 and a machine maximum of 1.7 × 1038. You can always extend the range of admissible argument values by using the homogeneity relations (6.11.22), below. #include <math.h> #include "nrutil.h" #define ERRTOL 0.05