176 Chapter 5. Evaluation of Functions Rational Functions You evaluate a rational function like P(z)p0+p1x+…+PxH R()=Q,=90+9nz++9c (5.3.4) in the obvious way,namely as two separate polynomials followed by a divide.As a matter of convention one usually chooses go =1,obtained by dividing numerator and denominator by any other go.It is often convenient to have both sets of 81 coefficients stored in a single array,and to have a standard function available for doing the evaluation: double ratval(double x,double cof[,int mm,int kk) Given mm,kk,and cof [0..mm+kk],evaluate and return the rational function (cof [O]+ Cam ICAL cof[]+..+cof [mm]xmm)/(1+cof [mm+1]x+...+cof [mm+kk]xkk). int j; RECIPES double sumd,sumn; Note precision!Change to float if desired for (sumn=cof [mm],j=mm-1;j>=0;j--)sumn=sumn*x+cof [j]; for (sumd=0.0,j=mm+kk;j>=mm+1;j--)sumd=(sumd+cof [j])*x; computer, return sumn/(1.0+sumd); Americ Press. 9 Programs CITED REFERENCES AND FURTHER READING: SCIENTIFIC Acton,F.S.1970,Numerica/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.Seminumerical Algorithms,2nd ed.,vol.2 of The Art of Computer Programming (Reading,MA:Addison-Wesley),84.6.[3] Fike,C.T.1968,Computer Evaluation of Mathematical Functions(Englewood Cliffs,NJ:Prentice- Hall),Chapter 4. 10-621 Winograd,S.1970,Communications on Pure and Applied Mathematics,vol.23,pp.165-179.[4] Kronsjo,L.1987,Algorithms:Their Complexity and Efficiency,2nd ed.(New York:Wiley).[5] uction, Numerical Recipes 43106 (outside North Software. 5.4 Complex Arithmetic As we mentioned in $1.2,the lack of built-in complex arithmetic in c is a nuisance for numerical work.Even in languages like FORTRAN that have complex data types,it is disconcertingly common to encounter complex operations that produce overflows or underflows when both the complex operands and 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
176 Chapter 5. Evaluation of 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). 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 convention one usually chooses q 0 = 1, obtained by dividing numerator 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 function available for doing the evaluation: double ratval(double x, double cof[], int mm, int kk) Given mm, kk, and cof[0..mm+kk], evaluate and return the rational function (cof[0] + cof[1]x + ··· + cof[mm]xmm)/(1 + cof[mm+1]x + ··· + cof[mm+kk]xkk). { int j; double sumd,sumn; Note precision! Change to float if desired. for (sumn=cof[mm],j=mm-1;j>=0;j--) sumn=sumn*x+cof[j]; for (sumd=0.0,j=mm+kk;j>=mm+1;j--) sumd=(sumd+cof[j])*x; return sumn/(1.0+sumd); } CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical 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, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), §4.6. [3] Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: PrenticeHall), Chapter 4. Winograd, S. 1970, Communications on Pure and Applied Mathematics, vol. 23, pp. 165–179. [4] Kronsj¨o, L. 1987, Algorithms: Their Complexity and Efficiency, 2nd ed. (New York: Wiley). [5] 5.4 Complex Arithmetic As we mentioned in §1.2, the lack of built-in complex arithmetic in C is a nuisance for numerical work. Even in languages like FORTRAN that have complex data types, it is disconcertingly common to encounter complex operations that produce overflows or underflows when both the complex operands and 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
5.4 Complex Arithmetic 177 Actually,complex arithmetic is not quite trivial.Addition and subtraction are done in the obvious way,performing the operation separately on the real and imaginary parts of the operands.Multiplication can also be done in the obvious way, with 4 multiplications,one addition,and one subtraction, (a+ib)(c+id)=(ac-bd)+i(bc+ad) (5.4.1) (the addition before the i doesn't count;it just separates the real and imaginary parts notationally).But it is sometimes faster to multiply via (a +ib)(c+id)=(ac-bd)+il(a+b)(c+d)-ac-bd (5.4.2) which has only three multiplications (ac,bd,(a+b)(c+d)).plus two additions and 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 形 overflow even when the final result is representable,this happens only when the final answer is on the edge of representability.Not so for the complex modulus,if you 令 are misguided enough to compute it as la+ibl Va2+62 (5.4.3) Press. (bad!) whose intermediate result will overflow if either a or b is as large as the square Programs root of the largest representable number(e.g.,1019 as compared to 1038).The right way to do the calculation is SCIENTIFIC 6 |a+b= fla1+(b/a严la≥bl 1Ib1+(a/b)2 (5.4.4) lal b Complex division should use a similar trick to prevent avoidable overflows, 1920 COMPUTING (ISBN underflow,or loss of precision, [a+b(dlcl+b-a(dlc】 装 Numerical 10521 a ib c+d(d/c) d≥ (5.4.5) 43108 c+id [a(c/d)++b(cld-a可 ldl ldl Recipes c(c/d)+d (outside Of course you should calculate repeated subexpressions,like c/d or d/c,only once. North Software. 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 0 c=d=0 1+vV1+(d/c)2 1c≥ld 0三 2 (5.4.6) lc/d+1+(c/d)2 lel la 2
5.4 Complex Arithmetic 177 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). Actually, complex arithmetic is not quite trivial. Addition and subtraction are done in the obvious way, performing the operation separately on the real and imaginary parts of the operands. Multiplication can also be done in the obvious way, with 4 multiplications, one addition, and one subtraction, (a + ib)(c + id)=(ac − bd) + i(bc + ad) (5.4.1) (the addition before the i doesn’t count; it just separates the real and imaginary parts 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) which has only three multiplications (ac, bd, (a + b)(c + d)), plus two additions and 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 overflow even when the final result is representable, this happens only when the final answer is on the edge of representability. Not so for the complex modulus, if you are misguided enough to compute it as |a + ib| = a2 + b2 (bad!) (5.4.3) whose intermediate result will overflow if either a or b is as large as the square root of the largest representable number (e.g., 10 19 as compared to 1038). The right way to do the calculation is |a + ib| = |a| 1+(b/a)2 |a|≥|b| |b| 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) Of course you should calculate repeated subexpressions, like c/d or d/c, only once. 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 |c| 1 + 1+(d/c)2 2 |c|≥|d| |d| |c/d| + 1+(c/d)2 2 |c| < |d| (5.4.6)
178 Chapter 5. Evaluation of Functions Then the answer is 0 w=0 d 2 w≠0，c≥0 Vc+id= 4 (5.4.7) iw w≠0，c<0,d≥0 2w d叫 iw w≠0，c<0,d<0 2w Routines implementing these algorithms are listed in Appendix C. nted for CITED REFERENCES AND FURTHER READING: Midy,P,and Yakovlev,Y.1991,Mathematics and Computers in Simulation,vol.33,pp.33-49. Knuth,D.E.1981.Seminumerical Algorithms,2nd ed.,vol.2 of The Art of Computer Programming (Reading,MA:Addison-Wesley)[see solutions to exercises 4.2.1.16 and 4.6.4.41]. (North University RECIPES Press. THE 5.5 Recurrence Relations and Clenshaw's Recurrence Formula Programs Many useful functions satisfy recurrence relations,e.g., SCIENTIFIC (m+1)Pn+1(x)=(2m+1)xPn(x)-nPn-1(x) (5.5.1) 6 +()=2n ()-n-1() (5.5.2) nEn+i(z)=e-x-IEn(I) (5.5.3) cosne =2 cos0 cos(n-1)0-cos(n-2)0 (5.5.4) Numerical sin ne =2cos0 sin(n-1)0-sin(n-2)0 (5.5.5) where the first three functions are Legendre polynomials,Bessel functions of the first (outside kind,and exponential integrals,respectively.(For notation see [1].)These relations are useful for extending computational methods from two successive values of n to North other values,either larger or smaller. Equations(5.5.4)and(5.5.5)motivate us to say a few words about trigonometric functions.If your program's running time is dominated by evaluating trigonometric functions,you are probably doing something wrong.Trig functions whose arguments form a linear sequence 0 =0o+no,n =0,1,2,...,are efficiently calculated by the following recurrence, cos(0+6)=cos0-[a cos0+Bsin] (5.5.6) sin(0+6)=sin0-a sin-3 cos
178 Chapter 5. Evaluation of 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). Then the answer is √ c + id = 0 w = 0 w + i d 2w w = 0, c ≥ 0 |d| 2w + iw w = 0, c < 0, d ≥ 0 |d| 2w − iw w = 0, c < 0, d < 0 (5.4.7) Routines implementing these algorithms are listed in Appendix C. CITED REFERENCES AND FURTHER READING: Midy, P., and Yakovlev, Y. 1991, Mathematics and Computers in Simulation, vol. 33, pp. 33–49. Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming (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) = (2n + 1)xPn(x) − nPn−1(x) (5.5.1) Jn+1(x) = 2n x Jn(x) − Jn−1(x) (5.5.2) nEn+1(x) = e−x − xEn(x) (5.5.3) cos nθ = 2 cos θ cos(n − 1)θ − cos(n − 2)θ (5.5.4) sin nθ = 2 cos θ sin(n − 1)θ − sin(n − 2)θ (5.5.5) where the first three functions are Legendre polynomials, Bessel functions of the first kind, and exponential integrals, respectively. (For notation see [1].) These relations are useful for extending computational methods from two successive values of n to other values, either larger or smaller. Equations (5.5.4) and (5.5.5) motivate us to say a few words about trigonometric functions. If your program’s running time is dominated by evaluating trigonometric functions, you are probably doing something wrong. Trig functions whose arguments form a linear sequence θ = θ0 + nδ, n = 0, 1, 2,... , are efficiently calculated by the following recurrence, cos(θ + δ) = cos θ − [α cos θ + β sin θ] sin(θ + δ) = sin θ − [α sin θ − β cos θ] (5.5.6)