f5-3

5.3PolynomialsandRationalFunctions 167Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Thompson,I.J.,andBarnett,A.R.1986, JournalofComputationalPhysics ,vol.64,pp.490–509. [5] Lentz, W.J. 1976, Applied Optics , vol. 15, pp. 668–671. [6] Jones, W.B. 1973, in Pad´e Approximants and Their Applications , P.R. Graves-Morris, ed. (Lon- don: Academic Press), p. 125. [7] 5.3 Polynomials and Rational Functions A polynomial of degree N−1is represented numerically as a stored array of coefficients, c(j)with j=1 ,...,N. We will always take c(1)to be the constant term in the polynomial, c(N)the coefficientof xN−1; but of course other conventions are possible. There are two kinds of manipulations that you can do with a polynomial: numerical manipulations (such as evaluation), where you are given the numerical value of its argument, or algebraic manipulations, where you want totransformthe coefficientarrayin someway withoutchoosinganyparticular argument. Let’s start with the numerical. We assume that youknow enough neverto evaluatea polynomialthis way: p=c(1)+c(2)*x+c(3)*x**2+c(4)*x**3+c(5)*x**4 Come the (computer) revolution, all persons found guilty of such criminal behavior will be summarily executed, and their programs won’t be! It is a matter of taste, however, whether to write p=c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*c(5)))) or p=(((c(5)*x+c(4))*x+c(3))*x+c(2))*x+c(1) If the number of coefficients is a large number n, one writes p=c(n) do11j=n-1,1,-1 p=p*x+c(j) enddo 11 Another useful trick is for evaluating a polynomial P(x)and its derivative dP(x)/dxsimultaneously: p=c(n) dp=0. do11j=n-1,1,-1 dp=dp*x+pp=p*x+c(j) enddo 11 which returns the polynomial as pand its derivative as dp. The above trick, which is basically synthetic division [1,2], generalizes to the evaluation of the polynomial and nd-1of its derivatives simultaneously: 168 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE ddpoly(c,nc,x,pd,nd) INTEGER nc,nd REAL x,c(nc),pd(nd) Given the coefficients of a polynomial of degree nc-1as an array c(1:nc) withc(1)being the constant term, and given a value x, and given a value nd>1, this routine returns the polynomial evaluated at xaspd(1) andnd-1derivatives as pd(2:nd) . INTEGER i,j,nndREAL constpd(1)=c(nc) do 11j=2,nd pd(j)=0. enddo 11 do13i=nc-1,1,-1 nnd=min(nd,nc+1-i) do12j=nnd,2,-1 pd(j)=pd(j)*x+pd(j-1) enddo 12 pd(1)=pd(1)*x+c(i) enddo 13 const=2. After the first derivative, factorial constants come in. do14i=3,nd pd(i)=const*pd(i)const=const*i enddo 14 return END As a curiosity, you might be interested to know that polynomials of degree n> 3can be evaluated in fewerthan nmultiplications, at least if you are willing to precompute some auxiliary coefficients and, in some cases, do an extra addition. For example, the polynomial P(x)= a0+a1x+a2x2+a3x3+a4x4(5.3.1 ) where a4>0, can be evaluatedwith 3 multiplicationsand5 additionsas follows: P(x)=[ ( Ax +B)2+Ax +C][(Ax +B)2+D]+E (5.3.2 ) where A, B, C, D, and Eare to be precomputed by A=(a4)1/4 B=a3−A3 4A3 D=3B2+8B3+a1A−2a2B A2 C=a2 A2−2B−6B2−D E=a0−B4−B2(C+D)−CD(5.3.3 ) Fifth degree polynomials can be evaluated in 4 multiplies and 5 adds; sixth degree polynomials can be evaluated in 4 multiplies and 7 adds; if any of this strikes you as interesting, consult references [3-5]. The subject has something of the same entertaining, if impractical, flavor as that of fast matrix multiplication, discussed in§2.11. 5.3PolynomialsandRationalFunctions 169Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Turnnowtoalgebraicmanipulations. Youmultiplyapolynomialofdegree n−1 (arrayof length n)by a monomialfactor x−abya bit of codelike thefollowing, c(n+1)=c(n) do11j=n,2,-1 c(j)=c(j-1)-c(j)*a enddo 11 c(1)=-c(1)*a Likewise,youdivideapolynomialofdegree n−1byamonomialfactor x−a (synthetic division again) using rem=c(n) c(n)=0. do11i=n-1,1,-1 swap=c(i)c(i)=remrem=swap+rem*a enddo 11 which leaves youwith a new polynomialarray and a numericalremainder rem. Multiplication of two general polynomials involves straightforward summing of the products, each involving one coefficient from each polynomial. Division oftwogeneralpolynomials,whileitcanbedoneawkwardlyinthefashiontaughtusing pencilandpaper,issusceptibletoagooddealofstreamlining. Witnessthefollowing routine based on the algorithm in [3]. SUBROUTINE poldiv(u,n,v,nv,q,r) INTEGER n,nvREAL q(n),r(n),u(n),v(nv) Given the ncoefficients of a polynomial in u(1:n) ,a n dt h e nvcoefficients of another polynomial in v(1:nv) , divide the polynomial uby the polynomial v(“u”/“v”)giving a quotient polynomial whose coefficients are returned in q(1:n-nv+1) , and a remainder polynomial whose coefficients are returned in r(1:nv-1) . The arrays qandrare dimen- sioned with lengths n, but the elements r(nv) ...r(n)andq(n-nv+2) ...q(n)will be returned as zero. INTEGER j,kdo 11j=1,n r(j)=u(j) q(j)=0. enddo 11 do13k=n-nv,0,-1 q(k+1)=r(nv+k)/v(nv)do 12j=nv+k-1,k+1,-1 r(j)=r(j)-q(k+1)*v(j-k) enddo 12 enddo 13 do14j=nv,n r(j)=0. enddo 14 returnEND 170 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Rational Functions You evaluate a rational function like R(x)=Pµ(x) Qν(x)=p0+p1x+··· +pµxµ q0+q1x+··· +qνxν(5.3.4 ) in the obvious way, namely as two separate polynomials followed by a divide. As a matter of conventionone usually chooses q0=1, obtainedby dividingnumerator and denominator by any other q0. It is often convenient to have both sets of coefficients stored in a single array, and to have a standard subroutine available for doing the evaluation: FUNCTION ratval(x,cof,mm,kk) INTEGER kk,mmDOUBLE PRECISION ratval,x,cof(mm+kk+1) Note precision! Change to REALif desired. Given mm,kk,a n d cof(1:mm+kk+1) , evaluate and return the rational function (cof(1) + cof(2)x +··· +cof(mm+1)xmm)/(1 + cof(mm+2)x +··· +cof(mm+kk+1)xkk). INTEGER j DOUBLE PRECISION sumd,sumnsumn=cof(mm+1) do 11j=mm,1,-1 sumn=sumn*x+cof(j) enddo 11 sumd=0.d0 do12j=mm+kk+1,mm+2,-1 sumd=(sumd+cof(j))*x enddo 12 ratval=sumn/(1.d0+sumd)returnEND CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 183, 190. [1] Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics , 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), pp. 361–363. [2] Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), §4.6. [3] Fike,C.T.1968, ComputerEvaluationofMathematicalFunctions (EnglewoodCliffs,NJ:Prentice- Hall), Chapter 4. Winograd,S.1970, CommunicationsonPureandAppliedMathematics ,vol.23,pp.165–179.[4] Kronsj¨o, L. 1987, Algorithms: TheirComplexity and Efficiency , 2nd ed. (NewYork: Wiley).[5] 5.4ComplexArithmetic 171Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).5.4 Complex Arithmetic Since FORTRAN has the built-in data type COMPLEX, you can generally let the compiler and intrinsic function library take care of complex arithmetic for you. Generally, but not always. For a program with only a small number of complex operations, you may want to code these yourself, in-line. Or, you may find that yourcompilerisnotuptosnuff: Itisdisconcertinglycommontoencountercomplex operations that produce overflows or underflows when both the complex operandsand the complex result are perfectly representable. This occurs, we think, because software companies assign inexperienced programmers to what they believe to be the perfectly trivial task of implementing complex arithmetic. Actually, complex arithmetic is not quitetrivial. Addition and subtraction are done in the obvious way, performing the operation separately on the real and imaginarypartsoftheoperands. Multiplicationcanalsobedoneintheobviousway, with 4 multiplications, one addition, and one subtraction, (a+ib)(c+id)=( ac−bd)+ i(bc+ad)( 5.4.1 ) (theadditionbeforethe idoesn’tcount;itjustseparatestherealandimaginaryparts notationally). But it is sometimes faster to multiply via (a+ib)(c+id)=( ac−bd)+ i[(a+b)(c+d)−ac−bd]( 5.4.2 ) whichhasonlythreemultiplications( ac,bd,(a+b)(c+d)),plustwoadditionsand three subtractions. The total operations count is higher by two, but multiplication is a slow operation on some machines. While it is true that intermediate results in equations (5.4.1) and (5.4.2) can overflowevenwhenthefinalresultisrepresentable,thishappensonlywhenthefinal answer is on the edge of representability. Not so for the complex modulus, if you or your compiler are misguided enough to compute it as |a+ib|=/radicalbig a2+b2 (bad!) (5.4.3 ) whose intermediate result will overflow if either aorbis as large as the square root of the largest representablenumber(e.g., 1019as comparedto 1038). The right way to do the calculation is |a+ib|=/braceleftbigg |a|/radicalbig 1+( b/a )2|a|≥| b| |b|/radicalbig 1+( a/b )2|a|<|b|(5.4.4 ) Complex division should use a similar trick to prevent avoidable overflows, underflow, or loss of precision, a+ib c+id=  [a+b(d/c )] + i[b−a(d/c )] c+d(d/c )|c|≥| d| [a(c/d )+ b]+i[b(c/d )−a] c(c/d )+ d|c|<|d|(5.4.5 )