5.3 Polynomials and Rational Functions 173 Thompson,I.J.,and Barnett,A.R.1986,Journal of Computationa/Physics,vol.64,pp.490-509. Lentz,W.J.1976,Applied Optics,vol.15,pp.668-671.[6] Jones,W.B.1973,in Pade Approximants and Their Applications,PR.Graves-Morris,ed.(Lon- don:Academic Press),p.125.[7] 5.3 Polynomials and Rational Functions A polynomial of degree N is represented numerically as a stored array of 83 coefficients,c[j]with j=0,...,N.We will always take c [o]to be the constant 鱼 18881892 term in the polynomial,c[N]the coefficient of N;but of course other conventions are possible.There are two kinds of manipulations that you can do with a polynomial: 100 numerical manipulations(such as evaluation),where you are given the numerical value of its argument,or algebraic manipulations,where you want to transform from NUMERICAL RECIPES I the coefficient array in some way without choosing any particular argument.Let's start with the numerical (Nort server We assume that you know enough never to evaluate a polynomial this way: America UnN电.t p=C[0]+c[1]*x+c[2]*x*x+c[3]*X*x*x+c[4]*x*X*x*x make one paper computer, ART or (even worse!). 9 Programs p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0); 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 OF SCIENTIFIC COMPUTING(ISBN p=c[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*c[4])); or v@cambri 10-621 p=((c[4]*x+c[3])*x+c[2])x+c[1])*x+c[0]; 1988-1992 by Numerical Recipes 43108 If the number of coefficients c[o..n]is large,one writes 68i-450J-）Ppxtct: Software. or (outside North America). p=c[j=n]; while (j>0)p=p*x+c[--j]; Another useful trick is for evaluating a polynomial P(x)and its derivative dP(x)/dx simultaneously: p=c[n]; dp=0.0; for(j=n-1;j>=0;j--){dp=dp*x+p;p=p*x+c[j];}
5.3 Polynomials and Rational Functions 173 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). Thompson, I.J., and Barnett, A.R. 1986, Journal of Computational Physics, 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. (London: Academic Press), p. 125. [7] 5.3 Polynomials and Rational Functions A polynomial of degree N is represented numerically as a stored array of coefficients, c[j] with j= 0,...,N. We will always take c[0] to be the constant term in the polynomial, c[N] the coefficient of xN ; 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 to transform the coefficient array in some way without choosing any particular argument. Let’s start with the numerical. We assume that you know enough never to evaluate a polynomial this way: p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x; or (even worse!), p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0); 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[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*c[4]))); or p=(((c[4]*x+c[3])*x+c[2])*x+c[1])*x+c[0]; If the number of coefficients c[0..n] is large, one writes p=c[n]; for(j=n-1;j>=0;j--) p=p*x+c[j]; or p=c[j=n]; while (j>0) p=p*x+c[--j]; Another useful trick is for evaluating a polynomial P(x) and its derivative dP(x)/dx simultaneously: p=c[n]; dp=0.0; for(j=n-1;j>=0;j--) {dp=dp*x+p; p=p*x+c[j];}
174 Chapter 5. Evaluation of Functions or p=c[j=n]; dp=0.0; while (j>0){dp=dp*x+p;p=p*x+c[--j];} which yields the polynomial as p and its derivative as dp. The above trick,which is basically synthetic division [1,2],generalizes to the evaluation of the polynomial and nd of its derivatives simultaneously: void ddpoly(float c[],int nc,float x,float pd[],int nd) Given the nc+1 coefficients of a polynomial of degree nc as an array c[0..nc]with c[o] being the constant term,and given a value x,and given a value nd>1,this routine returns the 9 polynomial evaluated at x as pd [o]and nd derivatives as pd [1..nd]. int nnd,j,i; float cnst=1.0; pd[o]=c[nc]; for(j=1;j=0;1--){ nnd=(nd (nc-i)?nd nc-i); for (j=nnd;j>=1;j--) (Nort pd[j]=pd[j]*x+pd[j-1]; pd [o]=pd[o]*x+c[i]; America computer, make one paper Press. THE ART for(i=2;13 can be evaluated in fewer than n multiplications,at least if you are willing to precompute some auxiliary coefficients and,in some cases,do an extra addition. For example,the polynomial P(r)=ao +aix a2x2+a3x3+aax (5.3.1) Numerica 10621 where a>0,can be evaluated with 3 multiplications and 5 additions as follows: 431 P(z)=[(Ax+B)2+Ax+C][(Ax+B)2+D]+E (5.3.2) (outside Recipes where A,B,C,D,and E are to be precomputed by North Software. A=(a4)1/A B=3-A3 4A3 D=3B2+8B3+11-2a2B (.3.3) 2 C=器-2B-6B-D E=ao-B4-B2(C+D)-CD
174 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). or p=c[j=n]; dp=0.0; while (j>0) {dp=dp*x+p; p=p*x+c[--j];} which yields the polynomial as p and its derivative as dp. The above trick, which is basically synthetic division [1,2], generalizes to the evaluation of the polynomial and nd of its derivatives simultaneously: void ddpoly(float c[], int nc, float x, float pd[], int nd) Given the nc+1 coefficients of a polynomial of degree nc as an array c[0..nc] with c[0] being the constant term, and given a value x, and given a value nd>1, this routine returns the polynomial evaluated at x as pd[0] and nd derivatives as pd[1..nd]. { int nnd,j,i; float cnst=1.0; pd[0]=c[nc]; for (j=1;j=0;i--) { nnd=(nd =1;j--) pd[j]=pd[j]*x+pd[j-1]; pd[0]=pd[0]*x+c[i]; } for (i=2;i 3 can be evaluated in fewer than n multiplications, 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 evaluated with 3 multiplications and 5 additions as follows: P(x) = [(Ax + B) 2 + Ax + C][(Ax + B) 2 + D] + E (5.3.2) where A, B, C, D, and E are 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)
5.3 Polynomials and Rational Functions 175 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-51.The subject has something of the same entertaining,if impractical,flavor as that of fast matrix multiplication,discussed ins2.11. Turn now to algebraic manipulations.You multiply a polynomial of degree n-1 (array of range [o..n-1])by a monomial factor x-a by a bit of code like the following, c[n]=c[n-1]; 19881992 for(j=n-1;j>=1;j-)c[0]=c[j-1]-c[j]*a; c[0]*=(-a); 拿 1-800 Likewise,you divide a polynomial of degree n by a monomial factor x-a from NUMERICAL RECIPESI (synthetic division again)using rem=c[n]; c[n]=0.0: (North America 州bMe se UnN电.t THE for(i=n-1;1>=0;1--){ to make one paper swap=c[i]; ART c[i]=rem; rem=swap+rem*a; strictly proh Programs which leaves you with a new polynomial array and a numerical remainder rem Multiplication of two general polynomials involves straightforward summing to dir of the products,each involving one coefficient from each polynomial.Division of two general polynomials,while it can be done awkwardly in the fashion taught using OF SCIENTIFIC COMPUTING(ISBN pencil and paper,is susceptible to a good deal of streamlining.Witness the following routine based on the algorithm in [3]. 10-621 void poldiv(float u[],int n,float v[],int nv,float q[],float r[]) Given the n+1 coefficients of a polynomial of degree n in u[0..n],and the nv+1 coefficients 1988-1992 by Numerical Recipes of another polynomial of degree nv in v[0..nv],divide the polynomial u by the polynomial 43106 v("u"/"v")giving a quotient polynomial whose coefficients are returned in q[0..n],and a remainder polynomial whose coefficients are returned in r[0..n].The elements r[nv..n] and q[n-nv+1..n]are returned as zero. (outside Software. int k,j; for (j=0;j=0;k--){ q[k]=r[nv+k]/v[nv]; for (j=nv+k-1;j>=k;j--)r[j]-=q[k]*v[j-k]; for (j=nv;j<=n;j++)r[j]=0.0;
5.3 Polynomials and Rational Functions 175 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). 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. Turn now to algebraic manipulations. You multiply a polynomial of degree n − 1 (array of range [0..n-1]) by a monomial factor x − a by a bit of code like the following, c[n]=c[n-1]; for (j=n-1;j>=1;j--) c[j]=c[j-1]-c[j]*a; c[0] *= (-a); Likewise, you divide a polynomial of degree n by a monomial factor x − a (synthetic division again) using rem=c[n]; c[n]=0.0; for(i=n-1;i>=0;i--) { swap=c[i]; c[i]=rem; rem=swap+rem*a; } which leaves you with a new polynomial array and a numerical remainder rem. Multiplication of two general polynomials involves straightforward summing of the products, each involving one coefficient from each polynomial. Division of two general polynomials, while it can be done awkwardly in the fashion taught using pencil and paper, is susceptible to a good deal of streamlining. Witness the following routine based on the algorithm in [3]. void poldiv(float u[], int n, float v[], int nv, float q[], float r[]) Given the n+1 coefficients of a polynomial of degree n in u[0..n], and the nv+1 coefficients of another polynomial of degree nv in v[0..nv], divide the polynomial u by the polynomial v (“u”/“v”)giving a quotient polynomial whose coefficients are returned in q[0..n], and a remainder polynomial whose coefficients are returned in r[0..n]. The elements r[nv..n] and q[n-nv+1..n] are returned as zero. { int k,j; for (j=0;j=0;k--) { q[k]=r[nv+k]/v[nv]; for (j=nv+k-1;j>=k;j--) r[j] -= q[k]*v[j-k]; } for (j=nv;j<=n;j++) r[j]=0.0; }
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