muller.c

Polynomial Root Finder is a reliable and fast C program (+ Matlab gateway) for finding all roots of

/**** suppress overflow *****/
void suppress_overflow(int nred)
/*int nred;        the highest exponent of the deflated polynomial        */
{
     int           kiter;            /* internal iteration counter        */
     unsigned char loop;             /* loop = FALSE => terminate loop    */
     double        help;             /* help variable                     */

     kiter = 0;                      /* reset iteration counter           */
     do { 
          loop=FALSE;                /* initial estimation: no overflow   */
          help = Cabs(x2);           /* help = |x2|                       */
          if (help>1. && fabs(nred*log10(help))>BOUND6) {
               kiter++;              /* if |x2|>1 and |x2|^nred>10^BOUND6 */ 
               if (kiter<KITERMAX) { /* then halve the distance between   */
                    h2=RCmul(.5,h2); /* new and old x2                    */
                    q2=RCmul(.5,q2); 
                    x2=Csub(x2,h2);
                    loop=TRUE;
               } else 
                    kiter=0;  
          }
     } while(loop);
}










/***** check of too big function values *****/
void too_big_functionvalues(double *f2absq)
/*double *f2absq;                                   f2absq=|f2|^2          */
{
     if ((fabs(f2.r)+fabs(f2.i))>BOUND4)         /* limit |f2|^2, when     */
          *f2absq = fabs(f2.r)+fabs(f2.i);       /* |f2.r|+|f2.i|>BOUND4   */
     else                                  
          *f2absq = (f2.r)*(f2.r)+(f2.i)*(f2.i); /* |f2|^2 = f2.r^2+f2.i^2 */
}










/***** Muller's modification to improve convergence *****/
void convergence_check(int *overflow,double f1absq,double f2absq,
                       double epsilon)
/*double f1absq,       f1absq = |f1|^2                          */
/*       f2absq,       f2absq = |f2|^2                          */
/*       epsilon;      bound for |q2|                           */
/*int    *overflow;    *overflow = TRUE  => overflow occures    */ 
/*                     *overflow = FALSE => no overflow occures */
{
     if ((f2absq>(CONVERGENCE*f1absq)) && (Cabs(q2)>epsilon) && 
     (iter<ITERMAX)) {
          q2 = RCmul(.5,q2); /* in case of overflow:            */
          h2 = RCmul(.5,h2); /* halve q2 and h2; compute new x2 */
          x2 = Csub(x2,h2);
          *overflow = TRUE;
     }
}










/***** compute P(x2) and make some checks *****/
void compute_function(dcomplex *pred,int nred,double f1absq,
                      double *f2absq,double epsilon)
/*dcomplex *pred;      coefficient vector of the deflated polynomial   */    
/*int      nred;       the highest exponent of the deflated polynomial */
/*double   f1absq,     f1absq = |f1|^2                                 */
/*         *f2absq,    f2absq = |f2|^2                                 */
/*         epsilon;    bound for |q2|                                  */
{
     int    overflow;    /* overflow = TRUE  => overflow occures       */ 
                         /* overflow = FALSE => no overflow occures    */

     do {
         overflow = FALSE; /* initial estimation: no overflow          */ 

                   /* suppress overflow                                */
         suppress_overflow(nred); 

                   /* calculate new value => result in f2              */
         fdvalue(pred,nred,&f2,&f2,x2,FALSE);

                   /* check of too big function values                 */
         too_big_functionvalues(f2absq);

                   /* increase iterationcounter                        */
         iter++;

                   /* Muller's modification to improve convergence     */
         convergence_check(&overflow,f1absq,*f2absq,epsilon);
     } while (overflow);
}










/***** is the new x2 the best approximation? *****/
void check_x_value(dcomplex *xb,double *f2absqb,int *rootd,
                   double f1absq,double f2absq,double epsilon,
                   int *noise)
/*dcomplex *xb;        best x-value                                    */
/*double   *f2absqb,   f2absqb |P(xb)|^2                               */
/*         f1absq,     f1absq = |f1|^2                                 */
/*         f2absq,     f2absq = |f2|^2                                 */
/*         epsilon;    bound for |q2|                                  */
/*int      *rootd,     *rootd = TRUE  => root determined               */
/*                     *rootd = FALSE => no root determined            */
/*         *noise;     noisecounter                                    */
{
     if ((f2absq<=(BOUND1*f1absq)) && (f2absq>=(BOUND2*f1absq))) {
                                  /* function-value changes slowly     */
          if (Cabs(h2)<BOUND3) {  /* if |h[2]| is small enough =>      */
              q2 = RCmul(2.,q2);  /* double q2 and h[2]                */
              h2 = RCmul(2.,h2);     
          } else {                /* otherwise: |q2| = 1 and           */
                                  /*            h[2] = h[2]*q2         */
              q2 = Complex(cos(iter),sin(iter));
              h2 = Cmul(h2,q2);
          }
     } else if (f2absq<*f2absqb) {
          *f2absqb = f2absq;      /* the new function value is the     */
          *xb      = x2;          /* best approximation                */
          *noise   = 0;           /* reset noise counter               */
          if ((sqrt(f2absq)<epsilon) && 
          (Cabs(Cdiv(Csub(x2,x1),x2))<epsilon))
               *rootd = TRUE;     /* root determined                   */
     }
}










/***** check, if determined root is good enough. *****/
void root_check(dcomplex *pred,int nred,double f2absqb,int *seconditer,
                int *rootd,int *noise,dcomplex xb)
/*dcomplex *pred,        coefficient vector of the deflated polynomial      */
/*         xb;           best x-value                                       */
/*int      nred,         the highest exponent of the deflated polynomial    */
/*         *noise,       noisecounter                                       */
/*         *rootd,       *rootd = TRUE  => root determined                  */
/*                       *rootd = FALSE => no root determined               */
/*         *seconditer;  *seconditer = TRUE  => start second iteration with */
/*                                              new initial estimations     */
/*                       *seconditer = FALSE => end routine                 */
/*double   f2absqb;      f2absqb |P(xb)|^2                                  */
{

     dcomplex df;     /* df=P'(x0)                                          */

     if ((*seconditer==1) && (f2absqb>0)) { 
          fdvalue(pred,nred,&f2,&df,xb,TRUE); /* f2=P(x0), df=P'(x0)        */
         if (Cabs(f2)/(Cabs(df)*Cabs(xb))>BOUND7) {
              /* start second iteration with new initial estimations        */
/*              x0 = Complex(-1./sqrt(2),1./sqrt(2)); 
              x1 = Complex(1./sqrt(2),-1./sqrt(2)); 
              x2 = Complex(-1./sqrt(2),-1./sqrt(2)); */
/*ml, 12-21-94: former initial values: */
              x0 = Complex(1.,0.);                 
              x1 = Complex(-1.,0.);                
              x2 = Complex(0.,0.);       /*   */
              fdvalue(pred,nred,&f0,&df,x0,FALSE); /* f0 =  P(x0)           */
              fdvalue(pred,nred,&f1,&df,x1,FALSE); /* f1 =  P(x1)           */
              fdvalue(pred,nred,&f2,&df,x2,FALSE); /* f2 =  P(x2)           */
              iter = 0;                /* reset iteration counter           */
              (*seconditer)++;         /* increase seconditer               */
              *rootd = FALSE;          /* no root determined                */
              *noise = 0;              /* reset noise counter               */
          }
     }
}