muller.c

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

/******************************************************************/
/*                                                                */
/* file:          muller.c                                        */
/*                                                                */
/* main function: muller()                                        */
/*                                                                */
/* version:       1.2                                             */
/*                                                                */
/* author:        B. Frenzel                                      */
/*                                                                */
/* date:          Jan 7 1993                                      */ 
/*                                                                */
/* input:         pred[]  coefficient vector of the deflated      */
/*                        polynomial                              */
/*                nred    the highest exponent of the deflated    */
/*                        polynomial                              */  
/*                                                                */
/* return:        xb      determined root                         */
/*                                                                */
/* subroutines:   initialize(),root_of_parabola(),                */
/*                iteration_equation(), compute_function(),       */
/*                check_x_value(), root_check()                   */
/*                                                                */
/* description:                                                   */
/* muller() determines the roots of a polynomial with complex     */
/* coefficients with the Muller method; these roots are the       */
/* initial estimations for the following Newton method            */
/*                                                                */
/* Copyright:                                                     */
/* Lehrstuhl fuer Nachrichtentechnik Erlangen                     */
/* Cauerstr. 7, 91054 Erlangen, FRG, 1993                         */
/* e-mail: [e-mail]                       */
/*                                                                */
/******************************************************************/
#define  MULLER
#include "header.h"











/***** main routine of Mueller's method *****/
dcomplex muller(dcomplex *pred,int nred)
/*dcomplex *pred;        coefficient vector of the deflated polynomial      */
/*int      nred;         the highest exponent of the deflated polynomial    */
{
     double   f1absq=FVALUE,  /* f1absq=|f1|^2                              */
              f2absq=FVALUE,  /* f2absq=|f2|^2                              */
              f2absqb=FVALUE, /* f2absqb=|P(xb)|^2                          */
              h2abs,          /* h2abs=|h2|                                 */
              epsilon;        /* bound for |q2|                             */
     int      seconditer=0,   /* second iteration, when root is too bad     */
              noise=0,        /* noise counter                              */
              rootd=FALSE;    /* rootd = TRUE  => root determined           */
                              /* rootd = FALSE => no root determined        */
     dcomplex xb;             /* best x-value                               */

                              /* initializing routine                       */
     initialize(pred,&xb,&epsilon); 

     fdvalue(pred,nred,&f0,&f0,x0,FALSE);   /* compute exact function value */
     fdvalue(pred,nred,&f1,&f1,x1,FALSE);   /* oct-29-1993 ml               */
     fdvalue(pred,nred,&f2,&f2,x2,FALSE);


do {                          /* loop for possible second iteration         */
     do {                     /* main iteration loop                        */
                              /* calculate the roots of the parabola        */
          root_of_parabola();
  
                              /* store values for the next iteration        */
          x0 = x1;  
          x1 = x2; 
          h2abs = Cabs(h2);   /* distance between x2 and x1                 */
          
                              /* main iteration-equation                    */
          iteration_equation(&h2abs);

                              /* store values for the next iteration        */ 
          f0 = f1;
          f1 = f2;
          f1absq = f2absq;

                              /* compute P(x2) and make some checks         */
          compute_function(pred,nred,f1absq,&f2absq,epsilon);

          /* printf("betrag %10.5e  %4.2d  %4.2d\n",f2absq,iter,seconditer);  */

                              /* is the new x-value the best approximation? */
          check_x_value(&xb,&f2absqb,&rootd,f1absq,f2absq,
                        epsilon,&noise);

                              /* increase noise counter                     */
         if (fabs((Cabs(xb)-Cabs(x2))/Cabs(xb))<NOISESTART)
              noise++;
     } while ((iter<=ITERMAX) && (!rootd) && (noise<=NOISEMAX));

     seconditer++;            /* increase seconditer                        */

                              /* check, if determined root is good enough   */
     root_check(pred,nred,f2absqb,&seconditer,&rootd,&noise,xb); 
} while (seconditer==2);

     return xb;               /* return best x value                        */
}










/***** initializing routine *****/
void initialize(dcomplex *pred,dcomplex *xb,double *epsilon)
/*dcomplex *pred,     coefficient vector of the deflated polynomial */    
/*         *xb;       best x-value                                  */
/*double   *epsilon;  bound for |q2|                                */
{
     /* initial estimations for x0,...,x2 and its values            */
     /* ml, 12-21-94 changed                                        */

     x0 = Complex(0.,0.);                 /* x0 = 0 + j*1           */ 
     x1 = Complex(-1./sqrt(2),-1./sqrt(2));                /* x1 = 0 - j*1           */
     x2 = Complex(1./sqrt(2),1./sqrt(2)); /* x2 = (1 + j*1)/sqrt(2) */

     h1 = Csub(x1,x0);                         /* h1 = x1 - x0      */
     h2 = Csub(x2,x1);                         /* h2 = x2 - x1      */
     q2 = Cdiv(h2,h1);                         /* q2 = h2/h1        */

     *xb      = x2;    /* best initial x-value = zero   */
     *epsilon = FACTOR*DBL_EPSILON;/* accuracy for determined root  */ 
     iter     = 0;                 /* reset iteration counter       */
}










/***** root_of_parabola() calculate smaller root of Muller's parabola *****/
void root_of_parabola(void)
{
     dcomplex A2,B2,C2,  /* variables to get q2                */
              discr,     /* discriminante                      */
              N1,N2;     /* denominators of q2                 */

                     /* A2 = q2(f2 - (1+q2)f1 + f0q2)          */
                     /* B2 = q2[q2(f0-f1) + 2(f2-f1)] + (f2-f1)*/
                     /* C2 = (1+q2)f[2]                        */
     A2   = Cmul(q2,Csub(Cadd(f2,Cmul(q2,f0)),
                         Cmul(f1,RCadd(1.,q2))));
     B2   = Cadd(Csub(f2,f1),Cmul(q2,Cadd(Cmul(q2,
                         Csub(f0,f1)),RCmul(2.,Csub(f2,f1)))));
     C2   = Cmul(f2,RCadd(1.,q2));
                     /* discr = B2^2 - 4A2C2                   */
     discr = Csub(Cmul(B2,B2),RCmul(4.,Cmul(A2,C2)));
                     /* denominators of q2                     */
     N1 = Csub(B2,Csqrt(discr));  
     N2 = Cadd(B2,Csqrt(discr));  
                 /* choose denominater with largest modulus    */
     if (Cabs(N1)>Cabs(N2) && Cabs(N1)>DBL_EPSILON)
          q2 = Cdiv(RCmul(-2.,C2),N1);  
     else if (Cabs(N2)>DBL_EPSILON)
          q2 = Cdiv(RCmul(-2.,C2),N2);  
     else 
          q2 = Complex(cos(iter),sin(iter));  
}










/***** main iteration equation: x2 = h2*q2 + x2 *****/
void iteration_equation(double *h2abs)
/*double *h2abs;                  Absolute value of the old distance        */
{
     double h2absnew,          /* Absolute value of the new h2              */
            help;              /* help variable                             */

     h2 = Cmul(h2,q2);         
     h2absnew = Cabs(h2);      /* distance between old and new x2           */

     if (h2absnew > (*h2abs*MAXDIST)) { /* maximum relative change          */ 
          help = MAXDIST/h2absnew;
          h2 = RCmul(help,h2);
          q2 = RCmul(help,q2); 
     } 

     *h2abs = h2absnew; /* actualize old distance for next iteration        */

     x2 = Cadd(x2,h2);
}