6.11 Elliptic Integrals and Jacobian Elliptic Functions 261 Other methods for computing Dawson's integral are also known [2.31. CITED REFERENCES AND FURTHER READING: Rybicki,G.B.1989,Computers in Physics,vol.3,no.2,pp.85-87.[1] Cody,W.J.,Pociorek,K.A.,and Thatcher,H.C.1970,Mathematics of Computation,vol.24, pp.171-178.2] McCabe,J.H.1974,Mathematics of Computation,vol.28.pp.811-816.[3] 6.11 Elliptic Integrals and Jacobian Elliptic Functions Elliptic integrals occur in many applications,because any integral of the form R(t,s）dt (6.11.1) where R is a rational function of t and s,and s is the square root of a cubic or 美e物s3间】 America quartic polynomial in t,can be evaluated in terms of elliptic integrals.Standard references [1]describe how to carry out the reduction,which was originally done by Legendre.Legendre showed that only three basic elliptic integrals are required. 9 The simplest of these is 心量二 OF SCIENTIFIC dt (6.11.2) (a1+b1t)(a2+b2t)(a3 +b3t)(a4+bat) 6 where we have written the quartic s2 in factored form.In standard integral tables [2], one of the limits of integration is always a zero of the quartic,while the other limit lies closer than the next zero,so that there is no singularity within the interval.To evaluate I1,we simply break the interval [y,into subintervals,each of which either 10-6211 begins or ends on a singularity.The tables,therefore,need only distinguish the eight cases in which each of the four zeros(ordered according to size)appears as the upper or lower limit of integration.In addition,when one of the b's in (6.11.2)tends to Numerical Recipes 431985 zero,the quartic reduces to a cubic,with the largest or smallest singularity moving to too:this leads to eight more cases (actually just special cases of the first eight) (outside The sixteen cases in total are then usually tabulated in terms of Legendre's standard North elliptic integral of the Ist kind,which we will define below.By a change of the variable of integration t,the zeros of the quartic are mapped to standard locations on the real axis.Then only two dimensionless parameters are needed to tabulate Legendre's integral.However,the symmetry of the original integral(6.11.2)under permutation of the roots is concealed in Legendre's notation.We will get back to Legendre's notation below.But first,here is a better way: Carlson [3]has given a new definition of a standard elliptic integral of the first kind, RF,y,司=互J0 1 dt (6.11.3) /(t+x)(t+y)(t+z)
6.11 Elliptic Integrals and Jacobian Elliptic Functions 261 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 machinereadable 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). Other methods for computing Dawson’s integral are also known [2,3]. CITED REFERENCES AND FURTHER READING: Rybicki, G.B. 1989, Computers in Physics, vol. 3, no. 2, pp. 85–87. [1] Cody, W.J., Pociorek, K.A., and Thatcher, H.C. 1970, Mathematics of Computation, vol. 24, pp. 171–178. [2] McCabe, J.H. 1974, Mathematics of Computation, vol. 28, pp. 811–816. [3] 6.11 Elliptic Integrals and Jacobian Elliptic Functions Elliptic integrals occur in many applications, because any integral of the form R(t, s) dt (6.11.1) where R is a rational function of t and s, and s is the square root of a cubic or quartic polynomial in t, can be evaluated in terms of elliptic integrals. Standard references [1] describe how to carry out the reduction, which was originally done by Legendre. Legendre showed that only three basic elliptic integrals are required. The simplest of these is I1 = x y dt (a1 + b1t)(a2 + b2t)(a3 + b3t)(a4 + b4t) (6.11.2) where we have written the quartic s2 in factored form. In standard integral tables [2], one of the limits of integration is always a zero of the quartic, while the other limit lies closer than the next zero, so that there is no singularity within the interval. To evaluate I1, we simply break the interval [y, x] into subintervals, each of which either begins or ends on a singularity. The tables, therefore, need only distinguish the eight cases in which each of the four zeros (ordered according to size) appears as the upper or lower limit of integration. In addition, when one of the b’s in (6.11.2) tends to zero, the quartic reduces to a cubic, with the largest or smallest singularity moving to ±∞; this leads to eight more cases (actually just special cases of the first eight). The sixteen cases in total are then usually tabulated in terms of Legendre’s standard elliptic integral of the 1st kind, which we will define below. By a change of the variable of integration t, the zeros of the quartic are mapped to standard locations on the real axis. Then only two dimensionless parameters are needed to tabulate Legendre’s integral. However, the symmetry of the original integral (6.11.2) under permutation of the roots is concealed in Legendre’s notation. We will get back to Legendre’s notation below. But first, here is a better way: Carlson [3] has given a new definition of a standard elliptic integral of the first kind, RF (x, y, z) = 1 2 ∞ 0 dt (t + x)(t + y)(t + z) (6.11.3)
262 Chapter 6.Special Functions where y,and z are nonnegative and at most one is zero.By standardizing the range of integration,he retains permutation symmetry for the zeros.(Weierstrass'canonical form also has this property.)Carlson first shows that when x or y is a zero of the quartic in (6.11.2),the integral Ii can be written in terms of RF in a form that is symmetric under permutation of the remaining three zeros.In the general case when neither nor y is a zero,two such RF functions can be combined into a single one by an addition theorem, leading to the fundamental formula h=2Rr(Ui,Ui,Ui4)） (6.11.4) where Uig =(XiXjYkYm+YiYjXkXm)/(x-y) (6.11.5) X:=(a:+b)412,Y.=(as+b)2 (6.11.6) and i,j,k,m is any permutation of 1,2,3,4.A short-cut in evaluating these expressions is Uis =Ui2-(a1ba-a4b1)(a2b3-a3b2) (6.11.7) Ui=U12-(a1b3-a3b1)(a2b4-a4b2) 的可 The U's correspond to the three ways of pairing the four zeros,and I is thus manifestly RECIPES symmetric under permutation of the zeros.Equation (6.11.4)therefore reproduces all sixteen 2 cases when one limit is a zero,and also includes the cases when neither limit is a zero Thus Carlson's function allows arbitrary ranges of integration and arbitrary positions of the branch points of the integrand relative to the interval of integration.To handle elliptic integrals of the second and third kind,Carlson defines the standard integral of the third kind as Rj(z,y,zp)= dt 20(t+p√t+x)t+)t+z) (6.11.8) 9 which is symmetric in z,y,and z.The degenerate case when two arguments are equal is denoted RD(x,y,2)=RJ(x,y,z,z) (6.11.9) and is symmetric in x and y.The function Ro replaces Legendre's integral of the second kind.The degenerate form of Rr is denoted Rc(,y)=RF(,y,y) (6.11.10) 寺2g》 It embraces logarithmic,inverse circular,and inverse hyperbolic functions. Carlson [4-7]gives integral tables in terms of the exponents of the linear factors of the quartic in (6.11.1).For example,the integral where the exponents are ( Numerica 10621 can be expressed as a single integral in terms of Rp;it accounts for 144 separate cases in Gradshteyn and Ryzhik[2]! 431 Refer to Carlson's papers [3-7]for some of the practical details in reducing elliptic Recipes integrals to his standard forms,such as handling complex conjugate zeros. Turn now to the numerical evaluation of elliptic integrals.The traditional methods [8] are Gauss or Landen transformations.Descending transformations decrease the modulus North k of the Legendre integrals towards zero,increasing transformations increase it towards unity.In these limits the functions have simple analytic expressions.While these methods converge quadratically and are quite satisfactory for integrals of the first and second kinds, they generally lead to loss of significant figures in certain regimes for integrals of the third kind.Carlson's algorithms [9.101,by contrast,provide a unified method for all three kinds with no significant cancellations. The key ingredient in these algorithms is the duplication theorem: RF(x,y,z)=2RF(x+入，y+入，z+入) =RF(巴+入，+A2+A (6.11.11) 4,4,4
262 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 machinereadable 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). where x, y, and z are nonnegative and at most one is zero. By standardizing the range of integration, he retains permutation symmetry for the zeros. (Weierstrass’ canonical form also has this property.) Carlson first shows that when x or y is a zero of the quartic in (6.11.2), the integral I1 can be written in terms of RF in a form that is symmetric under permutation of the remaining three zeros. In the general case when neither x nor y is a zero, two such RF functions can be combined into a single one by an addition theorem, leading to the fundamental formula I1 = 2RF (U2 12, U2 13, U2 14) (6.11.4) where Uij = (XiXjYkYm + YiYjXkXm)/(x − y) (6.11.5) Xi = (ai + bix) 1/2 , Yi = (ai + biy) 1/2 (6.11.6) and i, j, k, m is any permutation of 1, 2, 3, 4. A short-cut in evaluating these expressions is U2 13 = U2 12 − (a1b4 − a4b1)(a2b3 − a3b2) U2 14 = U2 12 − (a1b3 − a3b1)(a2b4 − a4b2) (6.11.7) The U’s correspond to the three ways of pairing the four zeros, and I1 is thus manifestly symmetric under permutation of the zeros. Equation (6.11.4) therefore reproduces all sixteen cases when one limit is a zero, and also includes the cases when neither limit is a zero. Thus Carlson’s function allows arbitrary ranges of integration and arbitrary positions of the branch points of the integrand relative to the interval of integration. To handle elliptic integrals of the second and third kind, Carlson defines the standard integral of the third kind as RJ (x, y, z, p) = 3 2 ∞ 0 dt (t + p) (t + x)(t + y)(t + z) (6.11.8) which is symmetric in x, y, and z. The degenerate case when two arguments are equal is denoted RD(x, y, z) = RJ (x, y, z, z) (6.11.9) and is symmetric in x and y. The function RD replaces Legendre’s integral of the second kind. The degenerate form of RF is denoted RC (x, y) = RF (x, y, y) (6.11.10) It embraces logarithmic, inverse circular, and inverse hyperbolic functions. Carlson [4-7] gives integral tables in terms of the exponents of the linear factors of the quartic in (6.11.1). For example, the integral where the exponents are (1 2 , 1 2 ,−1 2 ,−3 2 ) can be expressed as a single integral in terms of RD; it accounts for 144 separate cases in Gradshteyn and Ryzhik [2]! Refer to Carlson’s papers [3-7] for some of the practical details in reducing elliptic integrals to his standard forms, such as handling complex conjugate zeros. Turn now to the numerical evaluation of elliptic integrals. The traditional methods [8] are Gauss or Landen transformations. Descending transformations decrease the modulus k of the Legendre integrals towards zero, increasing transformations increase it towards unity. In these limits the functions have simple analytic expressions. While these methods converge quadratically and are quite satisfactory for integrals of the first and second kinds, they generally lead to loss of significant figures in certain regimes for integrals of the third kind. Carlson’s algorithms [9,10], by contrast, provide a unified method for all three kinds with no significant cancellations. The key ingredient in these algorithms is the duplication theorem: RF (x, y, z)=2RF (x + λ, y + λ, z + λ) = RF x + λ 4 , y + λ 4 , z + λ 4 (6.11.11)
6.11 Elliptic Integrals and Jacobian Elliptic Functions 263 where 入=(x)/2+(ez/P+(g2)2n (6.11.12) This theorem can be proved by a simple change of variable of integration[11].Equation (6.11.11)is iterated until the arguments of Rr are nearly equal.For equal arguments we have RF(e,x,x)=x-1/2 (6.11.13) When the arguments are close enough,the function is evaluated from a fixed Taylor expansion about(6.11.13)through fifth-order terms.While the iterative part of the algorithm is only linearly convergent,the error ultimately decreases by a factor of 4=4096 for each iteration. Typically only two or three iterations are required,perhaps six or seven if the initial values 81 of the arguments have huge ratios.We list the algorithm for Rr here,and refer you to Carlson's paper [9]for the other cases. a Stage 1:For n 0,1,2,...compute 4n=(xn十yn+zn)/3 Xn =1-(En/un),Yn=1-(yn/un),Zn=1-(zn/un) RECIPES En max(Xn,Yn,2n) (North 令 If en<tol go to Stage 2;else compute America computer, make one paper University Press. THE ART n=(EnUn)112+(Enzn)112+(Unzn)1/2 Programs xn+1=(xn+入n)/4,yn+1=(ym+入n)/4,zm+1=(zn+入n)/4 ictly proh and repeat this stage to dir Stage 2:Compute E2 XnYn-Zn;Es XnYnZn 1881992 OF SCIENTIFIC COMPUTING(ISBN RF=(1-六E2+六E+六E-是E2E3)/(un)2 In some applications the argument p in Ry or the argument y in Rc is negative,and the Numerical 10-521 Cauchy principal value of the integral is required.This is easily handled by using the formulas Recipes 43108 RJ(x,y之，p)= [(Y-y)RJ(E,y,2,7)-3RF(,y,z)+3Rc(Ez/y,PY/y)]/(y-p) (outside (6.11.14) Software. where Y三y+2-g-) (6.11.15) y-P is positive if p is negative,and 1/2 Rc(z,y)= Rc(-y,-y) (6.11.16) T-y The Cauchy principal value of Ry has a zero at some value of p<0,so (6.11.14)will give some loss of significant figures near the zero
6.11 Elliptic Integrals and Jacobian Elliptic Functions 263 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 machinereadable 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). where λ = (xy) 1/2 + (xz) 1/2 + (yz) 1/2 (6.11.12) This theorem can be proved by a simple change of variable of integration [11]. Equation (6.11.11) is iterated until the arguments of RF are nearly equal. For equal arguments we have RF (x, x, x) = x−1/2 (6.11.13) When the arguments are close enough, the function is evaluated from a fixed Taylor expansion about (6.11.13) through fifth-order terms. While the iterative part of the algorithm is only linearly convergent, the error ultimately decreases by a factor of 46 = 4096 for each iteration. Typically only two or three iterations are required, perhaps six or seven if the initial values of the arguments have huge ratios. We list the algorithm for RF here, and refer you to Carlson’s paper [9] for the other cases. Stage 1: For n = 0, 1, 2,... compute µn = (xn + yn + zn)/3 Xn = 1 − (xn/µn), Yn = 1 − (yn/µn), Zn = 1 − (zn/µn) n = max(|Xn|, |Yn|, |Zn|) If n < tol go to Stage 2; else compute λn = (xnyn) 1/2 + (xnzn) 1/2 + (ynzn) 1/2 xn+1 = (xn + λn)/4, yn+1 = (yn + λn)/4, zn+1 = (zn + λn)/4 and repeat this stage. Stage 2: Compute E2 = XnYn − Z2 n, E3 = XnYnZn RF = (1 − 1 10E2 + 1 14E3 + 1 24E2 2 − 3 44E2E3)/(µn) 1/2 In some applications the argument p in RJ or the argument y in RC is negative, and the Cauchy principal value of the integral is required. This is easily handled by using the formulas RJ (x, y,z, p) = [(γ − y)RJ(x, y, z, γ) − 3RF (x, y, z)+3RC (xz/y, pγ/y)] /(y − p) (6.11.14) where γ ≡ y + (z − y)(y − x) y − p (6.11.15) is positive if p is negative, and RC (x, y) = x x − y 1/2 RC(x − y,−y) (6.11.16) The Cauchy principal value of RJ has a zero at some value of p < 0, so (6.11.14) will give some loss of significant figures near the zero.
264 Chapter 6. Special Functions #include #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)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 #include "nrutil.h" #define ERRTOL 0.05
264 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 machinereadable 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 #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 nonnegative, 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) 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 constants, 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 #include "nrutil.h" #define ERRTOL 0.05
6.11 Elliptic Integrals and Jacobian Elliptic Functions 265 #define TINY 1.0e-25 #define BIG 4.5e21 #define C1 (3.0/14.0) #define C2 (1.0/6.0) 姓daf1na3 (9.0/22.0) #def1nec4(3.0/26.0) #define C5 (0.25*C3) #define C6 (1.5*C4) float rd(float x,float y,float z) Computes Carlson's elliptic integral of the second kind,Rp(z,y,2).x and y must be non- negative,and at most one can be zero.2 must be positive.TINY must be at least twice the http://www.n Permission is read able files negative 2/3 power of the machine overflow limit.BIG must be at most 0.1 x ERRTOL times the negative 2/3 power of the machine underflow limit. 83g float alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty (including this one) granted fori sqrtz,sum,xt,yt,zt; if (FMIN(x,y)BIG) 11-800-872 internet nrerror("invalid arguments in rd"); xt=x; yt-y; from NUMERICAL RECIPES IN 2t=2 sum=0.0; fac=1.0; do 7423 (North America sqrtx=sgrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sgrtx*(sqrty+sqrtz)+sqrty.sqrtz; to any server computer,is strictly prohibited. t users to make one paper copy for their THE sum +fac/(sqrtz*(zt+alamb)); fac=0.25*fac; xt=0.25*(xt+alamb) yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); own ave=0.2*(xt+yt+3.0*zt) delx=(ave-xt)/ave; dely=(ave-yt)/ave; only),or send email to directcustsen delz=(ave-zt)/ave; while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz))>ERRTOL); ART OF SCIENTIFIC COMPUTING(ISBN 0-521 ea=delx*dely; 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes eb=delz*delz; ecsea-eb; ed=ea-6.0*eb; eesedtectec; return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee) @cambridge.org personal use.Further reproduction,or +delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave)): 1-431085 Software. #include #include "nrutil.h" (outside North America) #define ERRTOL 0.05 #define TINY 2.5e-13 ying of machine #define BIG 9.0e11 #define C1 (3.0/14.0) #def1neC2(1.0/3.0) #define C3 (3.0/22.0) #define C4 (3.0/26.0) #define C5 (0.75*C3) #define C6 (1.5米C4 #define C7 (0.5*C2) #define C8 (C3+C3)
6.11 Elliptic Integrals and Jacobian Elliptic Functions 265 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 machinereadable 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). #define TINY 1.0e-25 #define BIG 4.5e21 #define C1 (3.0/14.0) #define C2 (1.0/6.0) #define C3 (9.0/22.0) #define C4 (3.0/26.0) #define C5 (0.25*C3) #define C6 (1.5*C4) float rd(float x, float y, float z) Computes Carlson’s elliptic integral of the second kind, RD(x, y, z). x and y must be nonnegative, and at most one can be zero. z must be positive. TINY must be at least twice the negative 2/3 power of the machine overflow limit. BIG must be at most 0.1 × ERRTOL times the negative 2/3 power of the machine underflow limit. { float alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty, sqrtz,sum,xt,yt,zt; if (FMIN(x,y) BIG) nrerror("invalid arguments in rd"); xt=x; yt=y; zt=z; sum=0.0; fac=1.0; do { sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; sum += fac/(sqrtz*(zt+alamb)); fac=0.25*fac; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); ave=0.2*(xt+yt+3.0*zt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; } while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz)) > ERRTOL); ea=delx*dely; eb=delz*delz; ec=ea-eb; ed=ea-6.0*eb; ee=ed+ec+ec; return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee) +delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave)); } #include #include "nrutil.h" #define ERRTOL 0.05 #define TINY 2.5e-13 #define BIG 9.0e11 #define C1 (3.0/14.0) #define C2 (1.0/3.0) #define C3 (3.0/22.0) #define C4 (3.0/26.0) #define C5 (0.75*C3) #define C6 (1.5*C4) #define C7 (0.5*C2) #define C8 (C3+C3)
266 Chapter 6. Special Functions float rj(float x,float y,float z,float p) Computes Carlson's elliptic integral of the third kind,RJ(,y,2,p).x.y,and 2 must be nonnegative,and at most one can be zero.p must be nonzero.If pBIG) nrerror("invalid arguments in rj"); 83 gm=0,0 fac=1.0; 1f(p>0.0) 11-800-872 (including this one) xt=xi internet yt=y; zt=z; pt=pi else xt=FMIN(FMIN(x,y),z) zt=FMAX(FMAX(x,y),z) yt=x+y+z-xt-zt; -7423 (North America to any server computer,is strictly prohibited. from NUMERICAL RECIPES IN C: THE a=1.0/(yt-p); b=a*(zt-yt)*(yt-xt); pt=yt+b; rho=xt*zt/yt; tau=p*pt/yt; rcx=rc(rho,tau); do sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sgrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; alpha-SQR(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz); beta=pt*SQR(pt+alamb); sum +fac*rc(alpha,beta); fac=0.25*fac; xt=0.25*(xt+alamb); only),orsend email to directcustserv@cambridge.org 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) yt=0.25*(yt+alamb): zt=0.25*(zt+alamb) pt=0.25*(pt+alamb); personal use.Further reproduction,or ave=0.2*(xt+yt+zt+pt+pt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; delp=(ave-pt)/ave: (outside North America) un Software. while (FMAX(FMAX(FMAX(fabs(delx),fabs(dely)), fabs(delz)),fabs(delp))>ERRTOL); ea=delx*(dely+delz)+dely*delz: eb=delx*dely*delz; ying of machine ec=delp*delp; ed=ea-3.0*ec; ee=eb+2.0*delp*(ea-ec); ans=3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4)) +delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sgrt(ave)); if (p <=0.0)ans=a*(b*ans+3.0*(rcx-rf(xt,yt,zt))); return ans;
266 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 machinereadable 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). float rj(float x, float y, float z, float p) Computes Carlson’s elliptic integral of the third kind, RJ (x, y, z, p). x, y, and z must be nonnegative, and at most one can be zero. p must be nonzero. If p BIG) nrerror("invalid arguments in rj"); sum=0.0; fac=1.0; if (p > 0.0) { xt=x; yt=y; zt=z; pt=p; } else { xt=FMIN(FMIN(x,y),z); zt=FMAX(FMAX(x,y),z); yt=x+y+z-xt-zt; a=1.0/(yt-p); b=a*(zt-yt)*(yt-xt); pt=yt+b; rho=xt*zt/yt; tau=p*pt/yt; rcx=rc(rho,tau); } do { sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; alpha=SQR(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz); beta=pt*SQR(pt+alamb); sum += fac*rc(alpha,beta); fac=0.25*fac; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); pt=0.25*(pt+alamb); ave=0.2*(xt+yt+zt+pt+pt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; delp=(ave-pt)/ave; } while (FMAX(FMAX(FMAX(fabs(delx),fabs(dely)), fabs(delz)),fabs(delp)) > ERRTOL); ea=delx*(dely+delz)+dely*delz; eb=delx*dely*delz; ec=delp*delp; ed=ea-3.0*ec; ee=eb+2.0*delp*(ea-ec); ans=3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4)) +delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave)); if (p <= 0.0) ans=a*(b*ans+3.0*(rcx-rf(xt,yt,zt))); return ans; }
6.11 Elliptic Integrals and Jacobian Elliptic Functions 267 #include #include "nrutil.h" #define ERRTOL 0.04 #define TINY 1.69e-38 #define SQRTNY 1.3e-19 #define BIG 3.e37 #define TNBG (TINYBIG) #define COMP1 (2.236/SQRTNY) #define COMP2 (TNBG*TNBG/25.0) #define THIRD (1.0/3.0) #define C1 0.3 #define C2(1.0/7.0) http://www.nr. #define C3 0.375 #define C4(9.0/22.0) 83 float rc(float x,float y) 鱼 granted for 19881992 Computes Carlson's degenerate elliptic integral,Rc(x,y).x must be nonnegative and y must be nonzero.If y 0,the Cauchy principal value is returned.TINY must be at least 5 times 11800 the machine underflow limit,BIG at most one fifth the machine maximum overflow limit. float alamb,ave,s,w,xt,yt; Cambridge from NUMERICAL RECIPESI if (x 0.0 II y =0.0 II (x+fabs(y))BIG II (y0.0&x0.0) xt=x; yt=y; server computer, (North America to make one paper UnN电.t THE w=1.0; ART 】e1se{ xt=x-yi 是 yt =-y; Programs w=sqrt(x)/sqrt(xt); send do alamb=2.0*sgrt(xt)*sqrt(yt)+yt; email st st Copyright (C) xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); ave=THIRD*(xt+yt+yt); s=(yt-ave)/ave; OF SCIENTIFIC COMPUTING(ISBN while (fabs(s)>ERRTOL); return w*(1.0+s*s*(C1+s*(C2+s*(C3+s*C4))))/sgrt(ave); v@cam 10-621 At times you may want to express your answer in Legendre's notation.Alter- natively,you may be given results in that notation and need to compute their values 1988-1992 by Numerical Recipes 43108 with the programs given above.It is a simple matter to transform back and forth. The Legendre elliptic integral of the Ist kind is defined as (outside do North Software. F(中，k)≡ (6.11.17) /1-k2sin20 The complete elliptic integral of the Ist kind is given by K(k)≡F(π/2，) (6.11.18) In terms of RF. F(,k)=sinoRg(cos2p,1-k2sin2中，1) (6.11.19) K(k)=RF(0,1-k2,1)
6.11 Elliptic Integrals and Jacobian Elliptic Functions 267 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 machinereadable 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 #include "nrutil.h" #define ERRTOL 0.04 #define TINY 1.69e-38 #define SQRTNY 1.3e-19 #define BIG 3.e37 #define TNBG (TINY*BIG) #define COMP1 (2.236/SQRTNY) #define COMP2 (TNBG*TNBG/25.0) #define THIRD (1.0/3.0) #define C1 0.3 #define C2 (1.0/7.0) #define C3 0.375 #define C4 (9.0/22.0) float rc(float x, float y) Computes Carlson’s degenerate elliptic integral, RC(x, y). x must be nonnegative and y must be nonzero. If y BIG || (y 0.0 && x 0.0) { xt=x; yt=y; w=1.0; } else { xt=x-y; yt = -y; w=sqrt(x)/sqrt(xt); } do { alamb=2.0*sqrt(xt)*sqrt(yt)+yt; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); ave=THIRD*(xt+yt+yt); s=(yt-ave)/ave; } while (fabs(s) > ERRTOL); return w*(1.0+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave); } At times you may want to express your answer in Legendre’s notation. Alternatively, you may be given results in that notation and need to compute their values with the programs given above. It is a simple matter to transform back and forth. The Legendre elliptic integral of the 1st kind is defined as F(φ, k) ≡ φ 0 dθ 1 − k2 sin2 θ (6.11.17) The complete elliptic integral of the 1st kind is given by K(k) ≡ F(π/2, k) (6.11.18) In terms of RF , F(φ, k) = sin φRF (cos2 φ, 1 − k2 sin2 φ, 1) K(k) = RF (0, 1 − k2, 1) (6.11.19)
268 Chapter 6.Special Functions The Legendre elliptic integral of the 2nd kind and the complete elliptic integral of the 2nd kind are given by E(o,k)V1-k2sin20 d0 sin oRr(cos2 ,1-k2 sin2 o,1) (6.11.20) -k2sin3oRp(cos2,1-k2 sin2 1) E(k)≡E(π/2，k)=RF(0,1-k2,1)-k2RD(0,1-k2,1) Finally,the Legendre elliptic integral of the 3rd kind is d Π(，n,k)≡ Jo (1+nsin20)V1-k2 sin20 from NUMERICAL RECIPESI 18881892 sin oRF(cos2,1-k2 sin2,1) (6.11.21) (Nort server 令 -in sin3 oRJ(cos2,1-k2 sin2 0,1,1+nsin2o) America Press. THE (Note that this sign convention for n is opposite that of Abramowitz and Stegun [121, ART and that their sina is our k. Progra #include #include "nrutil.h" float ellf(float phi,float ak) Legendre elliptic integral of the 1st kind F(,k),evaluated using Carlson's function Rp.The argument ranges are0≤中≤r/2,0≤ksinφ≤1. float rf(float x,float y,float z); OF SCIENTIFIC COMPUTING(ISBN float s; s=sin(phi); return s*rf(SQR(cos(phi)),(1.0-s*ak)*(1.0+s*ak),1.0); v@cambri 1988-1992 by Numerical Recipes 10-521 43108-5 #include #include "nrutil.h" float elle(float phi,float ak) (outside North Software. Legendre elliptic integral of the 2nd kind E(,)evaluated using Carlson's functions Rp and Rp.The argument ranges are0≤φ≤π/2,0≤ksin≤1. float rd(float x,float y,float z); float rf(float x,float y,float z); float cc,q,s; s=sin(phi); cc=SQR(cos(phi)); q=(1.0-s*ak)*(1.0+s*ak); return s*(rf(cc,q,1.0)-(SQR(s*ak))*rd(cc,q,1.0)/3.0);
268 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 machinereadable 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). The Legendre elliptic integral of the 2nd kind and the complete elliptic integral of the 2nd kind are given by E(φ, k) ≡ φ 0 1 − k2 sin2 θ dθ = sin φRF (cos2 φ, 1 − k2 sin2 φ, 1) − 1 3 k2 sin3 φRD(cos2 φ, 1 − k2 sin2 φ, 1) E(k) ≡ E(π/2, k) = RF (0, 1 − k2, 1) − 1 3 k2RD(0, 1 − k2, 1) (6.11.20) Finally, the Legendre elliptic integral of the 3rd kind is Π(φ, n, k) ≡ φ 0 dθ (1 + n sin2 θ) 1 − k2 sin2 θ = sin φRF (cos2 φ, 1 − k2 sin2 φ, 1) − 1 3n sin3 φRJ (cos2 φ, 1 − k2 sin2 φ, 1, 1 + n sin2 φ) (6.11.21) (Note that this sign convention for n is opposite that of Abramowitz and Stegun [12], and that their sinα is our k.) #include #include "nrutil.h" float ellf(float phi, float ak) Legendre elliptic integral of the 1st kind F(φ, k), evaluated using Carlson’s function RF . The argument ranges are 0 ≤ φ ≤ π/2, 0 ≤ k sin φ ≤ 1. { float rf(float x, float y, float z); float s; s=sin(phi); return s*rf(SQR(cos(phi)),(1.0-s*ak)*(1.0+s*ak),1.0); } #include #include "nrutil.h" float elle(float phi, float ak) Legendre elliptic integral of the 2nd kind E(φ, k), evaluated using Carlson’s functions RD and RF . The argument ranges are 0 ≤ φ ≤ π/2, 0 ≤ k sin φ ≤ 1. { float rd(float x, float y, float z); float rf(float x, float y, float z); float cc,q,s; s=sin(phi); cc=SQR(cos(phi)); q=(1.0-s*ak)*(1.0+s*ak); return s*(rf(cc,q,1.0)-(SQR(s*ak))*rd(cc,q,1.0)/3.0); }
6.11 Elliptic Integrals and Jacobian Elliptic Functions 269 #include #include "nrutil.h" float ellpi(float phi,float en,float ak) Legendre elliptic integral of the 3rd kind II(,n,k),evaluated using Carlson's functions RJ and RF.(Note that the sign convention on n is opposite that of Abramowitz and Stegun.The ranges of and k are0≤p≤r/2,0≤ksin≤l. float rf(float x,float y,float z); float rj(float x,float y,float z,float p); float cc,enss,q,si s=sin(phi); enss=en*s*s; cc=SQR(cos(phi)); 83g q=(1.0-s*ak)*(1.0+s*ak); granted for return s*(rf(cc,q,1.0)-enss*rj(cc,q,1.0,1.0+enss)/3.0); 1-.200 872 Carlson's functions are homogeneous of degree-and-3,so from NUMERICAL RECIPES I 19881992 Rr(Ax,λy,Xz)=入-1/2RF(x,y,z) (Nort server (6.11.22) R(Ax,入y,入z,Ap)=入-3/2R(z,y,z,p) America computer University Press. THE Thus to express a Carlson function in Legendre's notation,permute the first three arguments into ascending order,use homogeneity to scale the third argument to be Program 1,and then use equations (6.11.19)-(6.11.21). OF SCIENTIFIC( Jacobian Elliptic Functions 6 The Jacobian elliptic function sn is defined as follows:instead of considering the elliptic integral COMPUTING (ISBN u(,k)≡u=F(⊙，) (6.11.23) 1920 有 Numerical 10521 consider the inverse function 43108 y=sin中=sn(u,k) (6.11.24) Equivalently, (outside Recipes Software. sn North dy (6.11.25) 1-2y2)1-k22y2) Whenk=0,sn is just sin.The functions cn and dn are defined by the relations sn2 +cn2 =1,k2sn2 dn2 =1 (6.11.26) The routine given below actually takes me=2=1-k2 as an input parameter. It also computes all three functions sn,cn,and dn since computing all three is no harder than computing any one of them.For a description of the method,see [81
6.11 Elliptic Integrals and Jacobian Elliptic Functions 269 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 machinereadable 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 #include "nrutil.h" float ellpi(float phi, float en, float ak) Legendre elliptic integral of the 3rd kind Π(φ, n, k), evaluated using Carlson’s functions RJ and RF . (Note that the sign convention on n is opposite that of Abramowitz and Stegun.) The ranges of φ and k are 0 ≤ φ ≤ π/2, 0 ≤ k sin φ ≤ 1. { float rf(float x, float y, float z); float rj(float x, float y, float z, float p); float cc,enss,q,s; s=sin(phi); enss=en*s*s; cc=SQR(cos(phi)); q=(1.0-s*ak)*(1.0+s*ak); return s*(rf(cc,q,1.0)-enss*rj(cc,q,1.0,1.0+enss)/3.0); } Carlson’s functions are homogeneous of degree − 1 2 and −3 2 , so RF (λx, λy, λz) = λ−1/2RF (x, y, z) RJ (λx, λy, λz, λp) = λ−3/2RJ (x, y, z, p) (6.11.22) Thus to express a Carlson function in Legendre’s notation, permute the first three arguments into ascending order, use homogeneity to scale the third argument to be 1, and then use equations (6.11.19)–(6.11.21). Jacobian Elliptic Functions The Jacobian elliptic function sn is defined as follows: instead of considering the elliptic integral u(y, k) ≡ u = F(φ, k) (6.11.23) consider the inverse function y = sin φ = sn(u, k) (6.11.24) Equivalently, u = sn 0 dy (1 − y2)(1 − k2y2) (6.11.25) When k = 0, sn is just sin. The functions cn and dn are defined by the relations sn2 + cn2 = 1, k2sn2 + dn2 =1 (6.11.26) The routine given below actually takes mc ≡ k2 c = 1 − k2 as an input parameter. It also computes all three functions sn, cn, and dn since computing all three is no harder than computing any one of them. For a description of the method, see [8].
270 Chapter 6. Special Functions #include #define CA 0.0003 The accuracy is the square of CA void sncndn(float uu,float emmc,float *sn,float *cn,float *dn) Returns the Jacobian elliptic functions sn(u,ke),cn(u,ke),and dn(u,ke).Here uu =u,while emmc =2 float a,b,c,d,emc,ui float em[14],en[14]; int i,ii,l,bo; emc=emmc; Permission is Sample page usuu; if (emc){ bo=(emc=1;ii--)[ b=em[ii]; Copyright (C)1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) c*=(*dn); *dn=(en[ii]+a)/(b+a): a=c/b; a=1.0/sqrt(c*c+1.0); *sn=(*sn>=0.0？a:-a) *cn=c*(*sn); if (bo){ a=(*dn); only),orsend email to directcustserv@cambridge.org(outside North America). *dn=(*cn); Software. *cn=a; *sn /d; else ying of machine *cn=1.0/cosh(u); *dn=(*cn); 来gn=tanh (u)
270 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 machinereadable 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 #define CA 0.0003 The accuracy is the square of CA. void sncndn(float uu, float emmc, float *sn, float *cn, float *dn) Returns the Jacobian elliptic functions sn(u, kc), cn(u, kc), and dn(u, kc). Here uu = u, while emmc = k2 c . { float a,b,c,d,emc,u; float em[14],en[14]; int i,ii,l,bo; emc=emmc; u=uu; if (emc) { bo=(emc =1;ii--) { b=em[ii]; a *= c; c *= (*dn); *dn=(en[ii]+a)/(b+a); a=c/b; } a=1.0/sqrt(c*c+1.0); *sn=(*sn >= 0.0 ? a : -a); *cn=c*(*sn); } if (bo) { a=(*dn); *dn=(*cn); *cn=a; *sn /= d; } } else { *cn=1.0/cosh(u); *dn=(*cn); *sn=tanh(u); } }