/**** 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               */          }     }}