/******************************************************************//*                                                                *//* 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: int@nt.e-technik.uni-erlangen.de                       *//*                                                                *//******************************************************************/#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);}