9.5 Roots of Polynomials 373 The method operates iteratively:For a trial value z.a is calculated by equation (9.5.11).Then x-a becomes the next trial value.This continues until a is sufficiently small. The following routine implements the Laguerre method to find one root of a given polynomial of degree m,whose coefficients can be complex.As usual,the first coefficient a[o]is the constant term,while a [m]is the coefficient of the highest power of x.The routine implements a simplified version of an elegant stopping criterion due to Adams [5].which neatly balances the desire to achieve full machine accuracy.on the one hand,with the danger of iterating forever in the presence of roundoff error,on the other. #include <math.h> 19881992 #include "complex.h" #include "nrutil.h" 1.800 #define EPSS 1.0e-7 #define MR 8 #define MT 10 from NUMERICAL RECIPESI #define MAXIT (MT*MR) Here EPSS is the estimated fractional roundoff error.We try to break (rare)limit cycles with MR different fractional values,once every MT steps,for MAXIT total allowed iterations. server (Nort void laguer(fcomplex a[],int m,fcomplex *x,int *its) THE Given the degreemand the m1complex coefficients a[.]of the polynomial America computer and given a complex value x,this routine improves x by Laguerre's method until it converges, ART within the achievable roundoff limit,to a root of the given polynomial.The number of iterations taken is returned as its. 9 Programs int iter,j; float abx,abp,abm,err; fcomplex dx,x1,b,d,f,g,h,sq,gp,gm,g2; stat1cf1 oat frac[R+1]={0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0J: Fractions used to break a limit cycle. for (iter=1;iter<=MAXIT;iter++){ Loop over iterations up to allowed maximum. 米its=1ter; OF SCIENTIFIC COMPUTING(ISBN b=a[m]; 19881982 err=Cabs(b); d=f=Complex(0.0,0.0) v@cam abx=Cabs(*x); 10-621 for（j=m-1;j>=0;j--)( Efficient computation of the polynomial and f=Cadd(Cmul(*x,f),d); its first two derivatives.f stores P/2. Numerical Recipes -43108 d=Cadd(Cmul(*x,d),b); b=Cadd(Cmul(*x,b),a[j]) err-Cabs(b)+abxerr; (outside err米=EPSS North Software. Estimate of roundoff error in evaluating polynomial. if (Cabs(b)<=err)return; We are on the root. g=Cdiv(d,b); The generic case:use Laguerre's formula Ame ying of g2=Cmul(g,g); h=Csub(g2,RCmul(2.0,Cdiv(f,b))); sq=Csqrt(RCmul((float)(m-1),Csub(RCmul((float)m,h),g2))); gp=Cadd(g,sq); gm=Csub(g,sq): abp=Cabs(gp); abm=Cabs(gm); if (abp abm)gp=gm; dx=((FMAX(abp,abm)>0.0 Cdiv(Complex((float)m,0.0),gp) RCmul(1+abx,Complex(cos((float)iter),sin((float)iter))))); x1=Csub(*x,dx);9.5 Roots of Polynomials 373 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 machine￾readable 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). The method operates iteratively: For a trial value x, a is calculated by equation (9.5.11). Then x − a becomes the next trial value. This continues until a is sufficiently small. The following routine implements the Laguerre method to find one root of a given polynomial of degree m, whose coefficients can be complex. As usual, the first coefficient a[0] is the constant term, while a[m] is the coefficient of the highest power of x. The routine implements a simplified version of an elegant stopping criterion due to Adams [5], which neatly balances the desire to achieve full machine accuracy, on the one hand, with the danger of iterating forever in the presence of roundoff error, on the other. #include <math.h> #include "complex.h" #include "nrutil.h" #define EPSS 1.0e-7 #define MR 8 #define MT 10 #define MAXIT (MT*MR) Here EPSS is the estimated fractional roundoff error. We try to break (rare) limit cycles with MR different fractional values, once every MT steps, for MAXIT total allowed iterations. void laguer(fcomplex a[], int m, fcomplex *x, int *its) Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial m i=0 a[i]xi, and given a complex value x, this routine improves x by Laguerre’s method until it converges, within the achievable roundoff limit, to a root of the given polynomial. The number of iterations taken is returned as its. { int iter,j; float abx,abp,abm,err; fcomplex dx,x1,b,d,f,g,h,sq,gp,gm,g2; static float frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0}; Fractions used to break a limit cycle. for (iter=1;iter<=MAXIT;iter++) { Loop over iterations up to allowed maximum. *its=iter; b=a[m]; err=Cabs(b); d=f=Complex(0.0,0.0); abx=Cabs(*x); for (j=m-1;j>=0;j--) { Efficient computation of the polynomial and its first two derivatives. f stores P

f=Cadd(Cmul(*x,f),d); /2. d=Cadd(Cmul(*x,d),b); b=Cadd(Cmul(*x,b),a[j]); err=Cabs(b)+abx*err; } err *= EPSS; Estimate of roundoff error in evaluating polynomial. if (Cabs(b) <= err) return; We are on the root. g=Cdiv(d,b); The generic case: use Laguerre’s formula. g2=Cmul(g,g); h=Csub(g2,RCmul(2.0,Cdiv(f,b))); sq=Csqrt(RCmul((float) (m-1),Csub(RCmul((float) m,h),g2))); gp=Cadd(g,sq); gm=Csub(g,sq); abp=Cabs(gp); abm=Cabs(gm); if (abp < abm) gp=gm; dx=((FMAX(abp,abm) > 0.0 ? Cdiv(Complex((float) m,0.0),gp) : RCmul(1+abx,Complex(cos((float)iter),sin((float)iter))))); x1=Csub(*x,dx);