5.2 Evaluation of Continued Fractions 169 into equation (5.1.11),and then setting 2=1. Sometimes you will want to compute a function from a series representation even when the computation is not efficient.For example,you may be using the values obtained to fit the function to an approximating form that you will use subsequently (cf.$5.8).If you are summing very large numbers of slowly convergent terms,pay attention to roundoff errors!In floating-point representation it is more accurate to sum a list of numbers in the order starting with the smallest one,rather than starting with the largest one.It is even better to group terms pairwise,then in pairs of pairs, etc.,so that all additions involve operands of comparable magnitude. 81 CITED REFERENCES AND FURTHER READING: Goodwin,E.T.(ed.)1961,Modern Computing Methods,2nd ed.(New York:Philosophical Li- brary),Chapter 13 [van Wijngaarden's transformations].[1] Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). ICAL Chapter 3. Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55(Washington:National Bureau of Standards;reprinted 1968 by RECIPES Dover Publications,New York),$3.6. Mathews,J.,and Walker,R.L.1970,Mathematical Methods of Physics,2nd ed.(Reading,MA: W.A.Benjamin/Addison-Wesley),82.3.[2] g 令 Press. 5.2 Evaluation of Continued Fractions IENTIFIC Continued fractions are often powerful ways of evaluating functions that occur in scientific applications.A continued fraction looks like this: f(z)=6o+ (5.2.1) MPUTING 02 b1+ 1920 b2+ (ISBN 44 b3+ b4++ Numerica 10.621 Printers prefer to write this as 是 43108 a1 a2 a3 a4 a5 Recipes f()=b如+1+2+b朗+1+5+ (5.2.2) (outside In either(5.2.1)or(5.2.2),the a's and b's can themselves be functions of x,usually Software. 首 linear or quadratic monomials at worst (i.e.,constants times x or times z2).For example,the continued fraction representation of the tangent function is tamr=1-3-5-7- (5.2.3) Continued fractions frequently converge much more rapidly than power series expansions,and in a much larger domain in the complex plane (not necessarily including the domain of convergence of the series,however).Sometimes the continued fraction converges best where the series does worst,although this is not
5.2 Evaluation of Continued Fractions 169 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). into equation (5.1.11), and then setting z = 1. Sometimes you will want to compute a function from a series representation even when the computation is not efficient. For example, you may be using the values obtained to fit the function to an approximating form that you will use subsequently (cf. §5.8). If you are summing very large numbers of slowly convergent terms, pay attention to roundoff errors! In floating-point representation it is more accurate to sum a list of numbers in the order starting with the smallest one, rather than starting with the largest one. It is even better to group terms pairwise, then in pairs of pairs, etc., so that all additions involve operands of comparable magnitude. CITED REFERENCES AND FURTHER READING: Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library), Chapter 13 [van Wijngaarden’s transformations]. [1] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Chapter 3. Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), §3.6. Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics, 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), §2.3. [2] 5.2 Evaluation of Continued Fractions Continued fractions are often powerful ways of evaluating functions that occur in scientific applications. A continued fraction looks like this: f(x) = b0 + a1 b1 + a2 b2+ a3 b3+ a4 b4+ a5 b5+··· (5.2.1) Printers prefer to write this as f(x) = b0 + a1 b1 + a2 b2 + a3 b3 + a4 b4 + a5 b5 + ··· (5.2.2) In either (5.2.1) or (5.2.2), the a’s and b’s can themselves be functions of x, usually linear or quadratic monomials at worst (i.e., constants times x or times x2). For example, the continued fraction representation of the tangent function is tan x = x 1 − x2 3 − x2 5 − x2 7 − ··· (5.2.3) Continued fractions frequently converge much more rapidly than power series expansions, and in a much larger domain in the complex plane (not necessarily including the domain of convergence of the series, however). Sometimes the continued fraction converges best where the series does worst, although this is not
170 Chapter 5.Evaluation of Functions a general rule.Blanch [1]gives a good review of the most useful convergence tests for continued fractions. There are standard techniques,including the important quotient-difference algo- rithm,for going back and forth between continued fraction approximations,power series approximations,and rational function approximations.Consult Acton [2]for an introduction to this subject,and Fike [3]for further details and references. How do you tell how far to go when evaluating a continued fraction?Unlike a series,you can't just evaluate equation(5.2.1)from left to right,stopping when the change is small.Written in the form of(5.2.1),the only way to evaluate the continued fraction is from right to left,first (blindly!)guessing how far out to start.This is not the right way. The right way is to use a result that relates continued fractions to rational 菲 approximations,and that gives a means of evaluating(5.2.1)or(5.2.2)from left 超日支 to right.Let fn denote the result of evaluating (5.2.2)with coefficients through ICAL an and bn.Then fnBn An (5.2.4) RECIPES 9 where A,and B,are given by the following recurrence: A-1=1B-1≡0 A0≡b0B0≡1 Aj=bjAj-1+ajAj-2 Bj=bj Bj-1+ajBj-2 j=1,2,,n s9&公 (5.2.5) This method was invented by J.Wallis in 1655(!),and is discussed in his Arithmetica Infinitorum [41.You can easily prove it by induction. In practice,this algorithm has some unattractive features:The recurrence(5.2.5) frequently generates very large or very small values for the partial numerators and denominators A;and B;.There is thus the danger of overflow or underflow of the floating-point representation.However,the recurrence(5.2.5)is linear in the A's and 多名a Numerica 10621 B's.At any point you can rescale the currently saved two levels of the recurrence, e.g.,divide Aj,Bj,Aj-1,and Bj-1 all by Bj.This incidentally makes Aj=fj and is convenient for testing whether you have gone far enough:See if fi and fj-1 from the last iteration are as close as you would like them to be.(If B;happens to %9N be zero,which can happen,just skip the renormalization for this cycle.A fancier level of optimization is to renormalize only when an overflow is imminent,saving the unnecessary divides.All this complicates the program logic.) Two newer algorithms have been proposed for evaluating continued fractions. Steed's method does not use Ai and Bj explicitly,but only the ratio Dj=Bj-1/Bj One calculates Dj and Afi=fi-fj-1 recursively using D1=1/(b+aD5-1) (5.2.6) △f1=(b5D3-1)△f5-1 (5.2.7) Steed's method(see,e.g.,[51)avoids the need for rescaling of intermediate results. However,for certain continued fractions you can occasionally run into a situation
170 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). a general rule. Blanch [1] gives a good review of the most useful convergence tests for continued fractions. There are standard techniques, including the important quotient-difference algorithm, for going back and forth between continued fraction approximations, power series approximations, and rational function approximations. Consult Acton [2] for an introduction to this subject, and Fike [3] for further details and references. How do you tell how far to go when evaluating a continued fraction? Unlike a series, you can’t just evaluate equation (5.2.1) from left to right, stopping when the change is small. Written in the form of (5.2.1), the only way to evaluate the continued fraction is from right to left, first (blindly!) guessing how far out to start. This is not the right way. The right way is to use a result that relates continued fractions to rational approximations, and that gives a means of evaluating (5.2.1) or (5.2.2) from left to right. Let fn denote the result of evaluating (5.2.2) with coefficients through an and bn. Then fn = An Bn (5.2.4) where An and Bn are given by the following recurrence: A−1 ≡ 1 B−1 ≡ 0 A0 ≡ b0 B0 ≡ 1 Aj = bjAj−1 + ajAj−2 Bj = bjBj−1 + ajBj−2 j = 1, 2,...,n (5.2.5) This method was invented by J. Wallis in 1655 (!), and is discussed in his Arithmetica Infinitorum [4]. You can easily prove it by induction. In practice, this algorithm has some unattractive features: The recurrence (5.2.5) frequently generates very large or very small values for the partial numerators and denominators Aj and Bj . There is thus the danger of overflow or underflow of the floating-point representation. However, the recurrence (5.2.5) is linear in the A’s and B’s. At any point you can rescale the currently saved two levels of the recurrence, e.g., divide Aj , Bj , Aj−1, and Bj−1 all by Bj . This incidentally makes Aj = fj and is convenient for testing whether you have gone far enough: See if f j and fj−1 from the last iteration are as close as you would like them to be. (If Bj happens to be zero, which can happen, just skip the renormalization for this cycle. A fancier level of optimization is to renormalize only when an overflow is imminent, saving the unnecessary divides. All this complicates the program logic.) Two newer algorithms have been proposed for evaluating continued fractions. Steed’s method does not use Aj and Bj explicitly, but only the ratio Dj = Bj−1/Bj . One calculates Dj and ∆fj = fj − fj−1 recursively using Dj = 1/(bj + ajDj−1) (5.2.6) ∆fj = (bjDj − 1)∆fj−1 (5.2.7) Steed’s method (see, e.g., [5]) avoids the need for rescaling of intermediate results. However, for certain continued fractions you can occasionally run into a situation
5.2 Evaluation of Continued Fractions 171 where the denominator in (5.2.6)approaches zero,so that Di and Afi are very large.The next Afj+will typically cancel this large change,but with loss of accuracy in the numerical running sum of the fi's.It is awkward to program around this,so Steed's method can be recommended only for cases where you know in advance that no denominator can vanish.We will use it for a special purpose in the routine bessik (86.7). The best general method for evaluating continued fractions seems to be the modified Lentz's method [61.The need for rescaling intermediate results is avoided by using both the ratios Cj=Aj/Aj-1;Dj=Bj-1/Bj (5.2.8) and calculating f;by f方=fi-1C5D5 (5.2.9) 满 From equation(5.2.5),one easily shows that the ratios satisfy the recurrence relations D=1/(b+aD5-1),C3=b+a/Ci-1 (5.2.10 令 In this algorithm there is the danger that the denominator in the expression for Dj, or the quantity C;itself,might approach zero.Either of these conditions invalidates (5.2.10).However,Thompson and Barnett [5]show how to modify Lentz's algorithm to fix this:Just shift the offending term by a small amount,e.g.,10-30.If you work through a cycle of the algorithm with this prescription,you will see that fj+ 9」 is accurately calculated. In detail,the modified Lentz's algorithm is this: Set fo bo;if bo =0 set fo tiny ·Set Co=fo: to dir 。Set Do=0. ·Forj=1,2, OF SCIENTIFIC COMPUTING (ISBN Set Dj=bj ajDj-1. If Di=0,set Di tiny. Set Cj=bj aj/Cj-1. If Cj=0 set Cj=tiny. Set Di=1/Di v@cambridge.org 1988-1992 by Numerical Recipes Set△，=CiDj 12-:621-43106-50 Set fi=fi-1△i Ifl△j-l<eps then exit. (outside Here eps is your floating-point precision,say 10-7or 10-15.The parameter tiny North Software. should be less than typical values of eps say 10-30. The above algorithm assumes that you can terminate the evaluation of the continued fraction when f;-fi-1 is sufficiently small.This is usually the case, but by no means guaranteed.Jones [7]gives a list of theorems that can be used to justify this termination criterion for various kinds of continued fractions. There is at present no rigorous analysis oferror propagation in Lentz's algorithm. However,empirical tests suggest that it is at least as good as other methods
5.2 Evaluation of Continued Fractions 171 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). where the denominator in (5.2.6) approaches zero, so that D j and ∆fj are very large. The next ∆fj+1 will typically cancel this large change, but with loss of accuracy in the numerical running sum of the fj’s. It is awkward to program around this, so Steed’s method can be recommended only for cases where you know in advance that no denominator can vanish. We will use it for a special purpose in the routine bessik (§6.7). The best general method for evaluating continued fractions seems to be the modified Lentz’s method [6]. The need for rescaling intermediate results is avoided by using both the ratios Cj = Aj/Aj−1, Dj = Bj−1/Bj (5.2.8) and calculating fj by fj = fj−1CjDj (5.2.9) From equation (5.2.5), one easily shows that the ratios satisfy the recurrence relations Dj = 1/(bj + ajDj−1), Cj = bj + aj/Cj−1 (5.2.10) In this algorithm there is the danger that the denominator in the expression for D j , or the quantity Cj itself, might approach zero. Either of these conditions invalidates (5.2.10). However, Thompson and Barnett [5] show how to modify Lentz’s algorithm to fix this: Just shift the offending term by a small amount, e.g., 10 −30. If you work through a cycle of the algorithm with this prescription, you will see that f j+1 is accurately calculated. In detail, the modified Lentz’s algorithm is this: • Set f0 = b0; if b0 = 0 set f0 = tiny. • Set C0 = f0. • Set D0 = 0. • For j = 1, 2,... Set Dj = bj + ajDj−1. If Dj = 0, set Dj = tiny. Set Cj = bj + aj/Cj−1. If Cj = 0 set Cj = tiny. Set Dj = 1/Dj. Set ∆j = CjDj. Set fj = fj−1∆j . If |∆j − 1| < eps then exit. Here eps is your floating-point precision, say 10−7 or 10−15. The parameter tiny should be less than typical values of eps|bj|, say 10−30. The above algorithm assumes that you can terminate the evaluation of the continued fraction when |fj − fj−1| is sufficiently small. This is usually the case, but by no means guaranteed. Jones [7] gives a list of theorems that can be used to justify this termination criterion for various kinds of continued fractions. There is at present no rigorous analysis of error propagationin Lentz’s algorithm. However, empirical tests suggest that it is at least as good as other methods
172 Chapter 5.Evaluation of Functions Manipulating Continued Fractions Several important properties of continued fractions can be used to rewrite them in forms that can speed up numerical computation.An equivalence transformation an→an,bn→λbn:an+1一λan+1 (5.2.11) leaves the value of a continued fraction unchanged.By a suitable choice of the scale factor A you can often simplify the form of the a's and the b's.Of course,you can carry out successive equivalence transformations,possibly with different A's,on successive terms of the continued fraction. The even and odd parts of a continued fraction are continued fractions whose 菲 successive convergents are f2n and f2n+1,respectively.Their main use is that they converge twice as fast as the original continued fraction,and so if their terms are not much more complicated than the terms in the original there can be a big savings in computation.The formula for the even part of(5.2.2)is RECIPESI feven=dodd+ C1 C2 (5.2.12) where in terms of intermediate variables Press. b1 (5.2.13) 三兰。 an=bnbn-1 n≥2 we have d0=b0,c=a1,d1=1+a2 Cn=-a2n-1a2m-2,dn=1+a2m-1+a2m: n≥2 (5.2.14) You can find the similar formula for the odd part in the review by Blanch [1].Often x a combination of the transformations(5.2.14)and(5.2.11)is used to get the best form for numerical work. Numerical We will make frequent use of continued fractions in the next chapter. 经”> -431 CITED REFERENCES AND FURTHER READING: (outside Recipes Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- North matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York),83.10. Blanch,G.1964,S/AM Review,vol.6,pp.383-421.[1] Acton,F.S.1970,Numerica/Methods That Work,1990,corrected edition (Washington:Mathe- matical Association of America),Chapter 11.[2] Cuyt,A.,and Wuytack,L.1987,Nonlinear Methods in Numerical Analysis (Amsterdam:North- Holland),Chapter 1. Fike,C.T.1968,Computer Evaluation of Mathematical Functions(Englewood Cliffs,NJ:Prentice- Hall,ss8.2,10.4,and10.5.[3] Wallis,J.1695,in Opera Mathematica,vol.1,p.355.Oxoniae e Theatro Shedoniano.Reprinted by Georg Olms Verlag,Hildeshein,New York(1972).[4]
172 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). Manipulating Continued Fractions Several important properties of continued fractions can be used to rewrite them in forms that can speed up numerical computation. An equivalence transformation an → λan, bn → λbn, an+1 → λan+1 (5.2.11) leaves the value of a continued fraction unchanged. By a suitable choice of the scale factor λ you can often simplify the form of the a’s and the b’s. Of course, you can carry out successive equivalence transformations, possibly with different λ’s, on successive terms of the continued fraction. The even and odd parts of a continued fraction are continued fractions whose successive convergents are f2n and f2n+1, respectively. Their main use is that they converge twice as fast as the original continued fraction, and so if their terms are not much more complicated than the terms in the original there can be a big savings in computation. The formula for the even part of (5.2.2) is feven = d0 + c1 d1 + c2 d2 + ··· (5.2.12) where in terms of intermediate variables α1 = a1 b1 αn = an bnbn−1 , n ≥ 2 (5.2.13) we have d0 = b0, c1 = α1, d1 =1+ α2 cn = −α2n−1α2n−2, dn =1+ α2n−1 + α2n, n ≥ 2 (5.2.14) You can find the similar formula for the odd part in the review by Blanch [1]. Often a combination of the transformations (5.2.14) and (5.2.11) is used to get the best form for numerical work. We will make frequent use of continued fractions in the next chapter. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), §3.10. Blanch, G. 1964, SIAM Review, vol. 6, pp. 383–421. [1] Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 11. [2] Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: NorthHolland), Chapter 1. Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: PrenticeHall), §§8.2, 10.4, and 10.5. [3] Wallis, J. 1695, in Opera Mathematica, vol. 1, p. 355, Oxoniae e Theatro Shedoniano. Reprinted by Georg Olms Verlag, Hildeshein, New York (1972). [4]
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];}