c6-11 (779512), страница 2
Текст из файла (страница 2)
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) < 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 (7significant digits).
Since the error scales as 6n , we see that 0.0025 will yield double precision(16 significant digits) and require at most two or three more iterations. Since the coefficientsof the sixth-order truncation error are different for the other elliptic functions, these values forthe error tolerance should be changed to 0.04 and 0.0012 in the algorithm for RC , and 0.05 and0.0015 for RJ and RD . As well as being an algorithm in its own right for certain combinationsof 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 thecomputation. We have chosen conservative values, corresponding to a machine minimumof 3 × 10−39 and a machine maximum of 1.7 × 1038 .
You can always extend the range ofadmissible argument values by using the homogeneity relations (6.11.22), below.#include <math.h>#include "nrutil.h"#define ERRTOL 0.05Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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)6.11 Elliptic Integrals and Jacobian Elliptic FunctionsTINY 1.0e-25BIG 4.5e21C1 (3.0/14.0)C2 (1.0/6.0)C3 (9.0/22.0)C4 (3.0/26.0)C5 (0.25*C3)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 thenegative 2/3 power of the machine overflow limit. BIG must be at most 0.1 × ERRTOL timesthe 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) < 0.0 || FMIN(x+y,z) < TINY || FMAX(FMAX(x,y),z) > 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 <math.h>#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)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).#define#define#define#define#define#define#define#define265266Chapter 6.Special Functionsif (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(x+y,x+z),FMIN(y+z,fabs(p))) < TINY|| FMAX(FMAX(x,y),FMAX(z,fabs(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(fabs(delx),fabs(dely)),FMAX(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;}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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 benonnegative, and at most one can be zero. p must be nonzero. If p < 0, the Cauchy principalvalue is returned. TINY must be at least twice the cube root of the machine underflow limit,BIG at most one fifth the cube root of the machine overflow limit.{float rc(float x, float y);float rf(float x, float y, float z);float a,alamb,alpha,ans,ave,b,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee,fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt;6.11 Elliptic Integrals and Jacobian Elliptic Functions267float rc(float x, float y)Computes Carlson’s degenerate elliptic integral, RC (x, y).
x must be nonnegative and y mustbe nonzero. If y < 0, the Cauchy principal value is returned. TINY must be at least 5 timesthe machine underflow limit, BIG at most one fifth the machine maximum overflow limit.{float alamb,ave,s,w,xt,yt;if (x < 0.0 || y == 0.0 || (x+fabs(y)) < TINY || (x+fabs(y)) > BIG ||(y<-COMP1 && x > 0.0 && x < COMP2))nrerror("invalid arguments in rc");if (y > 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 valueswith the programs given above. It is a simple matter to transform back and forth.The Legendre elliptic integral of the 1st kind is defined asZF (φ, k) ≡0φpdθ1 − k 2 sin2 θ(6.11.17)The complete elliptic integral of the 1st kind is given byK(k) ≡ F (π/2, k)(6.11.18)In terms of RF ,F (φ, k) = sin φRF (cos2 φ, 1 − k 2 sin2 φ, 1)K(k) = RF (0, 1 − k 2 , 1)(6.11.19)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).#include <math.h>#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)268Chapter 6.Special FunctionsThe Legendre elliptic integral of the 2nd kind and the complete elliptic integral ofthe 2nd kind are given byZφE(φ, k) ≡p1 − k 2 sin2 θ dθ0(6.11.20)− 13 k 2 sin3 φRD (cos2 φ, 1 − k 2 sin2 φ, 1)E(k) ≡ E(π/2, k) = RF (0, 1 − k 2 , 1) − 13 k 2 RD (0, 1 − k 2 , 1)Finally, the Legendre elliptic integral of the 3rd kind isZφΠ(φ, n, k) ≡0dθp(1 + n sin θ) 1 − k 2 sin2 θ2= sin φRF (cos2 φ, 1 − k 2 sin2 φ, 1)(6.11.21)− 13 n sin3 φRJ (cos2 φ, 1 − k 2 sin2 φ, 1, 1 + n sin2 φ)(Note that this sign convention for n is opposite that of Abramowitz and Stegun [12],and that their sin α is our k.)#include <math.h>#include "nrutil.h"float ellf(float phi, float ak)Legendre elliptic integral of the 1st kind F (φ, k), evaluated using Carlson’s function RF .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
