f5-4

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 ) 172 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).Ofcourseyoushouldcalculaterepeatedsubexpressions,like c/dord/c,onlyonce. Complex square root is even more complicated, since we must both guard intermediate results, and also enforce a chosen branch cut (here taken to be the negative real axis). To take the square root of c+id, first compute w≡  0 c=d=0 /radicalbig |c|/radicalBigg 1+/radicalbig 1+( d/c )2 2|c|≥|d| /radicalbig |d|/radicalBigg |c/d|+/radicalbig 1+( c/d )2 2|c|<|d|(5.4.6 ) Then the answer is √ c+id=  0 w=0 w+i/parenleftbiggd 2w/parenrightbigg w/negationslash=0,c≥0 |d| 2w+iw w /negationslash=0,c< 0,d≥0 |d| 2w−iw w /negationslash=0,c< 0,d< 0(5.4.7 ) CITED REFERENCES AND FURTHER READING: Midy, P., andYakovlev, Y. 1991, Mathematics and Computers inSimulation , vol. 33, pp. 33–49. Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley) [see solutions to exercises 4.2.1.16 and 4.6.4.41]. 5.5 Recurrence Relations and Clenshaw’s Recurrence Formula Many useful functions satisfy recurrence relations, e.g., (n+1 )Pn+1(x)=( 2 n+1 )xP n(x)−nP n−1(x) (5.5.1) Jn+1(x)=2n xJn(x)−Jn−1(x) (5.5.2) nE n+1(x)=e−x−xE n(x) (5.5.3) cosnθ =2 c o s θcos(n−1)θ−cos(n−2)θ (5.5.4) sinnθ =2 c o s θsin(n−1)θ−sin(n−2)θ (5.5.5) wherethefirstthreefunctionsareLegendrepolynomials,Besselfunctionsofthefirst kind, and exponential integrals, respectively. (For notation see [1].) These relations