#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 a_0 + a_1x^1 +   + a_nx^n, and given a complex value x, this utine 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--){         f=Cadd(Cmul(*x,f),d);         d=Cadd(Cmul(*x,d),b);         b=Cadd(Cmul(*x,b),a[j]);         err=Cabs(b)+abx*err;      }      err *= EPSS;      /* Estimate the round off 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(exp(log(1+abx)),Complex(cos((float)iter),sin((float)iter)))));      x1=Csub(*x,dx);      if (x->r == x1.r && x->i == x1.i) return;     /* Converged. */      if (iter % MT) *x=x1;      else *x=Csub(*x,RCmul(frac[iter/MT],dx));      /* Every so often we take a fractional step, to break any limit cycle (itself a rare occurence). */   }   nrerror("too many iterations in laguer");   /* Very unusual - can occur only for complex roots. Try a different guess for the root. */   return;}