//----------------------------------------------------------------------------
//  EgalTech 2013-2014
//----------------------------------------------------------------------------
// File : JenkinsTraub.cpp          Data : 08.01.14        Versione : 1.5a1
// Contenuto : Implementazione calcolo degli zeri di polinomi a coefficienti
//             reali o complessi con il metodo di Jenkins e Traub.
//             Rpoly deriva da TOMS493. Cpoly deriva da TOMS419.
//
//
// Modifiche : 08.01.14 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "JenkinsTraub.h"

//--------------------------- Constants --------------------------------------
#define RADIX_DOUBLE          2
#define EPS_DOUBLE            2.22e-16
#define SMALL_DOUBLE          2.3e-308
#define BIG_DOUBLE            1.7e+308


//--------------------------- Class Rpoly --------------------------------------
//------------------------------------------------------------------------------
// IN:  op - double precision vector of coefficients in order of decreasing powers.
//      degree - integer degree of polynomial
// OUT: zeror,zeroi - output double precision vectors of the real and imaginary parts of the zeros.
// RET: -1 if leading coefficient is zero, otherwise number of roots found. 
//------------------------------------------------------------------------------
int
Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) 
{
   bool   bZerOk ;
   bool   bScale ;
   bool   bContinue ;
   int    cnt, nz, i, j, jj, l, nm1 ;
   double t, aa, bb, cc, factor, rot ;
   double lo, max, min, xx, yy, cosr, sinr, xxx, x, sc, bnd ;
   double xm, ff, df, dx, infin, smalno, base ;
   double temp[POLY_MAXDEG+1] ;
   double pt[POLY_MAXDEG+1] ;


  // The following statements set machine constants.
   base = RADIX_DOUBLE ;
   eta = EPS_DOUBLE ;
   infin = BIG_DOUBLE ;
   smalno = SMALL_DOUBLE ;
   are = eta ;
   mre = eta ;
   lo = smalno / eta ;

  //  Initialization of constants for shift rotation.
   xx = sqrt( 0.5) ;
   yy = - xx ;
   rot = 94.0 ;
   rot *= 0.017453293 ;
   cosr = cos( rot) ;
   sinr = sin( rot) ;
   n = degree ;

  // Inizializzo numero iterazioni
   itercnt = 0 ;

  // Algorithm fails if the leading coefficient is zero.
   if ( op[0] == 0.0)
      return -1 ;

  // Remove the zeros at the origin, if any.
   while ( op[n] == 0.0) {
      j = degree - n ;
      zeror[j] = 0.0 ;
      zeroi[j] = 0.0 ;
      n -- ;
   }
   if ( n < 1)
      return degree ;

  // Make a copy of the coefficients.
   for ( i = 0 ; i <= n ; i++)
      p[i] = op[i] ;

  //  Start the algorithm for one zero.
   bContinue = true ;
   while ( bContinue) {    

     // Calculate the final zero 
      if ( n == 1) {
         zeror[degree-1] = - p[1] / p[0] ;
         zeroi[degree-1] = 0.0 ;
         n -= 1 ;
         return ( degree - n) ;
      }

     // Calculate a pair of zeros.
      if ( n == 2) {
         quad( p[0], p[1], p[2], &zeror[degree-2], &zeroi[degree-2],
                                 &zeror[degree-1], &zeroi[degree-1]) ;
         n -= 2 ;
         return ( degree - n) ;
      }

     // Find largest and smallest moduli of coefficients.
      max = 0.0 ;
      min = infin ;
      for ( i = 0 ; i <= n ; i++) {
         x = abs( p[i]) ;
         if ( x > max)
            max = x ;
         if ( x != 0.0  &&  x < min)
            min = x ;
      }
     // Scale if there are large or very small coefficients.
     // Computes a scale factor to multiply the coefficients of the
     // polynomial. The scaling si done to avoid overflow and to
     // avoid undetected underflow interfering with the convergence
     // criterion. The factor is a power of the base.
      bScale = true ;
      sc = lo / min ;
      if ( sc > 1.0  &&  ( infin / sc) < max)
         bScale = false ;
      if ( sc <= 1.0) {
         if ( max < 10.0)
            bScale = false ;
         if ( sc == 0.0)
            sc = smalno ;
      }
      if ( bScale) {
        // Scale polynomial.
         l = (int)( log( sc) / log( base) + 0.5) ;
         factor = pow( base * 1.0, l) ;
         if ( factor != 1.0) {
            for ( i = 0 ; i <= n ; i ++) 
               p[i] = factor * p[i] ;
         }
      }

     // Compute lower bound on moduli of roots.
      for ( i = 0 ; i <= n ; i++) {
         pt[i] = abs( p[i]) ;
      }
      pt[n] = - pt[n] ;
     // Compute upper estimate of bound.
      x = exp( ( log( - pt[n]) - log( pt[0])) / (double) n) ;
     // If Newton step at the origin is better, use it.
      if ( pt[n-1] != 0.0) {
         xm = - pt[n] / pt[n-1] ;
         if ( xm < x)
            x = xm ;
      }

     // Chop the interval (0,x) until ff <= 0
      while ( true) {
         xm = x * 0.1 ;
         ff = pt[0] ;
         for ( i = 1 ; i <= n ; i ++) 
            ff = ff * xm + pt[i] ;
         if ( ff <= 0.0)
            break ;
         x = xm ;
      }

     // Do Newton iteration until x converges to two decimal places.
      dx = x ;
      while ( abs( dx / x) > 0.005) {
         ff = pt[0] ;
         df = ff ;
         for ( i = 1 ; i < n ; i++) { 
            ff = ff * x + pt[i] ;
            df = df * x + ff ;
         }
         ff = ff * x + pt[n] ;
         dx = ff / df ;
         x -= dx ;
         itercnt ++ ;
      }
      bnd = x ;

     // Compute the derivative as the initial k polynomial and do 5 steps with no shift.
      nm1 = n - 1 ;
      for ( i = 1 ; i < n ; i++)
         k[i] = (double)( n - i) * p[i] / (double) n ;
      k[0] = p[0] ;
      aa = p[n] ;
      bb = p[n-1] ;
      bZerOk = ( k[n-1] == 0) ;
      for ( jj = 0 ; jj < 5 ; jj ++) {
         itercnt ++ ;
         cc = k[n-1] ;
        // Use a scaled form of recurrence if value of k at 0 is nonzero.
         if ( ! bZerOk) {
            t = - aa / cc ;
            for ( i = 0 ; i < nm1 ; i ++) {
               j = n - i - 1 ;
               k[j] = t * k[j-1] + p[j] ;
            }
            k[0] = p[0] ;
            bZerOk = ( abs( k[n-1]) <= abs( bb) * eta * 10.0) ;
         }
         else {
           // Use unscaled form of recurrence.
            for ( i = 0 ; i < nm1 ; i ++) {
               j = n - i - 1 ;
               k[j] = k[j-1] ;
            }
            k[0] = 0.0 ;
            bZerOk = ( k[n-1] == 0.0) ;
         }
      }

     // Save k for restarts with new shifts.
      for ( i = 0 ; i < n ; i ++) 
         temp[i] = k[i] ;

     // Loop to select the quadratic corresponding to each new shift.
      bContinue = false ;
      for ( cnt = 0 ; cnt < 20 ; cnt ++) {
        //  Quadratic corresponds to a double shift to a non-real point and its complex conjugate.
        //  The point has modulus bnd and amplitude rotated by 94 degrees from the previous shift.
         xxx = cosr * xx - sinr * yy ;
         yy = sinr * xx + cosr * yy ;
         xx = xxx ;
         sr = bnd * xx ;
         si = bnd * yy ;
         u = -2.0 * sr ;
         v = bnd ;
         fxshfr( 20 * ( cnt + 1), &nz);
        // The second stage jumps directly to one of the third stage iterations and returns here if successful.
        // Deflate the polynomial, store the zero or zeros and return to the main algorithm.
         if ( nz != 0) {
            j = degree - n ;
            zeror[j] = szr ;
            zeroi[j] = szi ;
            n -= nz ;
            for ( i = 0 ; i <= n ; i ++)
               p[i] = qp[i] ;
            if ( nz != 1) {
               zeror[j+1] = lzr ;
               zeroi[j+1] = lzi ;
            }
           // interrompo il loop e riparto a cercare un nuovo zero
            bContinue = true ;
            break ;
         }
        // If the iteration is unsuccessful another quadratic is chosen after restoring k.
         else {
            for ( i = 0 ; i < n ; i ++)
               k[i] = temp[i] ;
         }
      }
   }

  // Return with failure if no convergence after 20 shifts.
   return ( degree - n) ;
}

//------------------------------------------------------------------------------
//  Computes up to L2 fixed shift k-polynomials,
//  testing for convergence in the linear or quadratic
//  case. Initiates one of the variable shift
//  iterations and returns with the number of zeros found.
//------------------------------------------------------------------------------
void
Rpoly::fxshfr( int l2, int* nz)
{
   bool   bVpass ;
   bool   bSpass ;
   bool   bVtry ;
   bool   bStry ;
   bool   bFflag ;
   bool   bIflag ;
   int    type, i, j ;
   double svu, svv, ui, vi, s ;
   double betas, betav, oss, ovv, ss, vv, ts, tv ;
   double ots, otv, tvv, tss ;
   double svk[POLY_MAXDEG+1] ;


  // Inizializzazioni
   bIflag = true ;
   *nz = 0 ;
   betav = 0.25 ;
   betas = 0.25 ;
   oss = sr ;
   ovv = v ;

  // Evaluate polynomial by synthetic division.
   quadsd( n, &u, &v, p, qp, &a, &b) ;
   calcsc( &type) ;

   for ( j = 0 ; j < l2 ; j ++) {

      bFflag = true ;

     // Calculate next k polynomial and estimate v.
      nextk( &type) ;
      calcsc( &type) ;
      newest( type, &ui, &vi) ;
      vv = vi ;

     // Estimate s.
      ss = 0.0 ;
      if ( k[n-1] != 0.0)
         ss = - p[n] / k[n-1] ;
      tv = 1.0 ;
      ts = 1.0 ;
      if ( j == 0  ||  type == 3) {
         ovv = vv ;
         oss = ss ;
         otv = tv ;
         ots = ts ;
         continue ;
      }

     // Compute relative measures of convergence of s and v sequences.
      if ( vv != 0.0)
         tv = abs( ( vv - ovv) / vv) ;
      if ( ss != 0.0)
         ts = abs( ( ss - oss) / ss) ;
     /*  If decreasing, multiply two most recent convergence measures. */
      tvv = 1.0 ;
      if ( tv < otv)
         tvv = tv * otv ;
      tss = 1.0 ;
      if ( ts < ots)
         tss = ts * ots ;
     // Compare with convergence criteria.
      bVpass = ( tvv < betav) ;
      bSpass = ( tss < betas) ;
      if ( ! ( bSpass  ||  bVpass)) {
         ovv = vv ;
         oss = ss ;
         otv = tv ;
         ots = ts ;
         continue ;
      }

     // At least one sequence has passed the convergence test. Store variables before iterating.
      svu = u;
      svv = v;
      for ( i = 0 ; i < n ; i ++)
         svk[i] = k[i] ;
      s = ss ;

     // Choose iteration according to the fastest converging sequence.
      bVtry = false ;
      bStry = false ;
      while ( true) {
         if ( bFflag  &&  bSpass  &&  ( ! bVpass  ||  tss < tvv))
            ;
         else {
            quadit( &ui, &vi, nz) ;
            if ( *nz > 0)
               return ;
           // Quadratic iteration has failed. Flag that it has been tried and decrease the convergence criterion.
            bIflag = true ;
            bVtry = true ;
            betav *= 0.25;
           // Try linear iteration if it has not been tried and the S sequence is converging.
            if ( bStry  ||  ! bSpass)
               bIflag = false ;
            else {
               for ( i = 0 ; i < n ; i ++)
                  k[i] = svk[i] ;
            }
         }
         bFflag = false ;
         if ( bIflag) {
            realit( s, nz, &bIflag) ;
            if ( *nz > 0)
               return ;
           // Linear iteration has failed. Flag that it has been tried and decrease the convergence criterion.
            bStry = true ;
            betas *= 0.25 ;
           // If linear iteration signals an almost double real zero attempt quadratic iteration.
            if ( bIflag) {
               ui = -( s + s) ;
               vi = s * s ;
               break ;
            }
         }
        // Restore variables
         u = svu;
         v = svv;
         for ( i = 0 ; i < n ; i ++) {
            k[i] = svk[i] ;
         }
        // Try quadratic iteration if it has not been tried and the V sequence is convergin.
         if ( ! bVpass  ||  bVtry)
            break ;
      }

     // Recompute QP and scalar values to continue the second stage.
      quadsd( n, &u, &v, p, qp, &a, &b) ;
      calcsc( &type) ;

     // Salvo valori come precedenti
      ovv = vv ;
      oss = ss ;
      otv = tv ;
      ots = ts ;
    }
}

//------------------------------------------------------------------------------
//  Variable-shift k-polynomial iteration for a
//  quadratic factor converges only if the zeros are
//  equimodular or nearly so.
//  uu, vv - coefficients of starting quadratic.
//  nz - number of zeros found.
//------------------------------------------------------------------------------
void
Rpoly::quadit( double *uu, double *vv, int *nz)
{
   bool   bTried ;
   int    type, i, j ;
   double ui, vi ;
   double mp, omp, ee, relstp, t, zm ;


  // Inizializzazioni
   *nz = 0 ;
   bTried = false ;
   u = *uu ;
   v = *vv ;
   j = 0 ;

  // Main loop.
   while ( true) {
      itercnt ++ ;

      quad( 1.0, u, v, &szr, &szi, &lzr, &lzi) ;

     // Return if roots of the quadratic are real and not
     // close to multiple or nearly equal and of opposite sign.
      if ( abs( abs( szr) - abs( lzr)) > 0.01 * abs( lzr))
         return ;

     // Evaluate polynomial by quadratic synthetic division.
      quadsd( n, &u, &v, p, qp, &a, &b) ;
      mp = abs( a - szr * b) + abs( szi * b) ;
     // Compute a rigorous bound on the rounding error in evaluating p.
      zm = sqrt( abs( v)) ;
      ee = 2.0 * abs( qp[0]) ;
      t = -szr * b ;
      for ( i = 1 ; i < n ; i ++) {
         ee = ee * zm + abs( qp[i]) ;
      }
      ee = ee * zm + abs( a + t) ;
      ee *= (5.0 * mre + 4.0 * are) ;
      ee = ee - ( 5.0 * mre + 2.0 * are) * ( abs( a + t) + abs( b) * zm) ;   
      ee = ee + 2.0 * are * abs( t) ;
     // Iteration has converged sufficiently if the polynomial value is less than 20 times this bound.
      if ( mp <= 20.0 * ee) {
         *nz = 2 ;
         return ;
      }
      j ++ ;
    
     // Stop iteration after 20 steps.
      if ( j > 20)
         return ;

     // A cluster appears to be stalling the convergence.
     // Five fixed shift steps are taken with a u,v close to the cluster.
      if ( j >= 2  &&
           ! ( relstp > 0.01  ||  mp < omp  ||  bTried)) {
         if ( relstp < eta)
            relstp = eta;
         relstp = sqrt( relstp) ;
         u = u - u * relstp ;
         v = v + v * relstp ;
         quadsd( n, &u, &v, p, qp, &a, &b) ;
         for ( i = 0 ; i < 5 ; i ++) {
	         calcsc( &type) ;
	         nextk( &type) ;
         }
         bTried = true ;
         j = 0 ;
      }

     // Salvo valore
      omp = mp ;

     // Calculate next k polynomial and new u and v.
      calcsc( &type) ;
      nextk( &type) ;
      calcsc( &type) ;
      newest( type, &ui, &vi) ;
     // If vi is zero the iteration is not converging.
      if ( vi == 0.0)
         return ;
      relstp = abs( ( vi - v) / vi) ;
      u = ui ;
      v = vi ;
   }
}

//------------------------------------------------------------------------------
//  Variable-shift H polynomial iteration for a real zero.
//  sss - starting iterate
//  nz  - number of zeros found
//  iflag - flag to indicate a pair of zeros near real axis.
//------------------------------------------------------------------------------
void
Rpoly::realit( double sss, int* nz, bool* pbIflag)
{
   int    i, j ;
   double pv, kv, t, s ;
   double ms, mp, omp, ee ;


  // Inizializzazioni
   *nz = 0 ;
   s = sss ;
   *pbIflag = false ;
   j = 0 ;

  // Main loop
   while ( true) {
      itercnt ++ ;
      pv = p[0] ;
     // Evaluate p at s.
      qp[0] = pv ;
      for ( i = 1 ; i <= n ; i++) {
         pv = pv * s + p[i] ;
         qp[i] = pv ;
      }
      mp = abs( pv) ;
     // Compute a rigorous bound on the error in evaluating p.
      ms = abs( s) ;
      ee = ( mre / ( are + mre)) * abs( qp[0]) ;
      for ( i = 1 ; i <= n ; i ++) {
         ee = ee * ms + abs( qp[i]) ;
      }
     // Iteration has converged sufficiently if the polynomial value is less than 20 times this bound.
      if ( mp <= 20.0 * (( are + mre) * ee - mre * mp)) {
         *nz = 1 ;
         szr = s ;
         szi = 0.0 ;
         return ;
      }
      j ++ ;
     // Stop iteration after 10 steps.
      if ( j > 10)
         return ;
     // A cluster of zeros near the real axis has been encountered.
      if ( j >= 2  &&
           ! ( abs( t) > 0.001 * abs( s-t)  ||  mp < omp)) {
        // Return with iflag set to initiate a quadratic iteration.
         *pbIflag = true ;
         return ;
      }

     // Return if the polynomial value has increased significantly.
      omp = mp ;

     // Compute t, the next polynomial, and the new iterate.
      kv = k[0] ;
      qk[0] = kv ;
      for ( i = 1 ; i < n ; i ++) {
         kv = kv*s + k[i] ;
         qk[i] = kv;
      }
      if ( abs( kv) <= abs( k[n-1]) * 10.0 * eta) {
        // Use unscaled form.
         k[0] = 0.0 ;
         for ( i = 1 ; i < n ; i ++) {
            k[i] = qk[i-1] ;
         }
      }
      else {
        // Use the scaled form of the recurrence if the value of k at s is nonzero.
         t = - pv / kv ;
         k[0] = qp[0] ;
         for ( i = 1 ; i < n ; i++) {
            k[i] = t * qk[i-1] + qp[i] ;
         }
      }
      kv = k[0] ;
      for ( i = 1 ; i < n ; i ++) {
         kv = kv * s + k[i] ;
      }
      t = 0.0 ;
      if ( abs( kv) > ( abs( k[n-1] * 10.0 * eta)))
         t = - pv / kv ;
      s += t ;
   }
}

//------------------------------------------------------------------------------
//  This routine calculates scalar quantities used to
//  compute the next k polynomial and new estimates of
//  the quadratic coefficients.
//  type - integer variable set here indicating how the
//  calculations are normalized to avoid overflow.
//------------------------------------------------------------------------------
void
Rpoly::calcsc( int *type)
{
  // Synthetic division of k by the quadratic 1,u,v
   quadsd( n-1, &u, &v, k, qk, &c, &d) ;

  // Type=3 indicates the quadratic is almost a factor of k.
   if ( abs( c) <= abs( k[n-1] * 100.0 * eta)  &&
        abs( d) <= abs( k[n-2] * 100.0 * eta)) {
      *type = 3 ;
      return ;
   }

  // Type=1 indicates that all formulas are divided by c.
   if ( abs( d) < abs( c)) {
      *type = 1 ;
      e = a / c ;
      f = d / c ;
      g = u * e ;
      h = v * b ;
      a3 = a * e + ( h / c + g) * b ;
      a1 = b - a * ( d / c) ;
      a7 = a + g * d + h * f ;
      return ;
   }

  // Type=2 indicates that all formulas are divided by d.
   *type = 2 ;
   e = a / d ;
   f = c / d ;
   g = u * b ;
   h = v * b ;
   a3 = ( a + g) * e + h * ( b / d) ;
   a1 = b * f - a ;
   a7 = ( f + u)*a + h ;
}

//------------------------------------------------------------------------------
//  Computes the next k polynomials using scalars computed in calcsc.
//------------------------------------------------------------------------------
void
Rpoly::nextk( int* type)
{
   double temp ;
   int i ;


   if ( *type == 3) {
    /*  Use unscaled form of the recurrence if type is 3. */
       k[0] = 0.0 ;
       k[1] = 0.0 ;
       for ( i = 2 ; i < n ; i ++) {
           k[i] = qk[i-2] ;
       }
       return ;
   }
   temp = a ;
   if ( *type == 1)
      temp = b ;
   if ( abs( a1) <= abs( temp) * eta * 10.0) {
     // If a1 is nearly zero then use a special form of the recurrence.
       k[0] = 0.0;
       k[1] = -a7*qp[0] ;
       for ( i = 2 ; i < n ; i ++) {
          k[i] = a3 * qk[i-2] - a7 * qp[i-1] ;
       }
       return ;           // HVE return added
   }
   /*  Use scaled form of the recurrence. */
   a7 /= a1 ;
   a3 /= a1 ;
   k[0] = qp[0] ;
   k[1] = qp[1] - a7 * qp[0] ;
   for ( i = 2 ; i < n ; i ++) {
      k[i] = a3 * qk[i-2] - a7 * qp[i-1] + qp[i] ;
   }
}

//------------------------------------------------------------------------------
//  Compute new estimates of the quadratic coefficients using the scalars computed in calcsc.
//------------------------------------------------------------------------------
void
Rpoly::newest( int type, double *uu, double *vv)
{
   double a4,a5,b1,b2,c1,c2,c3,c4,temp;


  /* Use formulas appropriate to setting of type. */
   if ( type == 3) {
     /*  If type=3 the quadratic is zeroed. */
      *uu = 0.0 ;
      *vv = 0.0 ;
      return ;
   }
   if ( type == 2) {
      a4 = ( a + g) * f + h ;
      a5 = ( f + u) * c + v * d ;
   }
   else {
      a4 = a + u * b + h * f ;
      a5 = c + ( u + v * f) * d ;
   }
  /*  Evaluate new quadratic coefficients. */
   b1 = -k[n-1] / p[n] ;
   b2 = -( k[n-2] + b1 * p[n-1]) / p[n] ;
   c1 = v * b2 * a1 ;
   c2 = b1 * a7 ;
   c3 = b1 * b1 * a3 ;
   c4 = c1 - c2 - c3 ;
   temp = a5 + b1 * a4 - c4 ;
   if ( temp == 0.0) {
      *uu = 0.0 ;
      *vv = 0.0 ;
      return ;
   }
   *uu = u - ( u * ( c3 + c2) + v * ( b1 * a1 + b2 * a7)) / temp ;
   *vv = v * ( 1.0 + c4 / temp) ;
   return ;
}

//------------------------------------------------------------------------------
// Divides p by the quadratic 1,u,v placing the quotient in q and the remainder in a,b.
//------------------------------------------------------------------------------
void
Rpoly::quadsd( int nn, double *u, double *v, double *p, double *q,
               double *a, double *b)
{
   int i ;
   double c ;


   *b = p[0] ;
   q[0] = *b ;
   *a = p[1] - (*b) * (*u) ;
   q[1] = *a ;
   for ( i = 2 ; i <= nn ; i++) {
      c = p[i] - (*a) * (*u) - (*b) * (*v) ;
      q[i] = c ;
      *b = *a ;
      *a = c ;
   }
}

//------------------------------------------------------------------------------
//  Calculate the zeros of the quadratic a*z^2 + b1*z + c.
//  The quadratic formula, modified to avoid overflow, is used 
//  to find the larger zero if the zeros are real and both
//  are complex. The smaller real zero is found directly from 
//  the product of the zeros c/a.
//------------------------------------------------------------------------------
void
Rpoly::quad( double a, double b1, double c,
             double* sr, double* si, double* lr, double* li)
{
   double b, d, e ;


   if ( a == 0.0) {         /* less than two roots */
      if ( b1 != 0.0)
         *sr = - c / b1 ;
      else
         *sr = 0.0 ;
      *lr = 0.0 ;
      *si = 0.0 ;
      *li = 0.0 ;
      return;
   }
   if ( c == 0.0) {         /* one real root, one zero root */
      *sr = 0.0 ;
      *lr = - b1 / a ;
      *si = 0.0 ;
      *li = 0.0 ;
      return;
   }
  /* Compute discriminant avoiding overflow. */
   b = b1 / 2.0 ;
   if ( abs( b) < abs( c)) { 
      if ( c < 0.0)
         e = - a ;
      else
         e = a ;
      e = b * ( b / abs( c)) - e ;
      d = sqrt( abs( e)) * sqrt( abs( c)) ;
   }
   else {
      e = 1.0 - ( a / b) *( c / b) ;
      d = sqrt( abs( e)) * abs( b) ;
   }
   if ( e < 0.0) {      /* complex conjugate zeros */
      *sr = - b / a ;
      *lr = *sr ;
      *si = abs( d / a) ;
      *li = - ( *si) ;
   }
   else {
      if ( b >= 0.0)
         d = - d ;       /* real zeros. */
      *lr = ( - b + d) / a ;
      *sr = 0.0 ;
      if ( *lr != 0.0)
         *sr = ( c / *lr) / a ;
      *si = 0.0 ;
      *li = 0.0 ;
   }
}


//--------------------------- Class Cpoly --------------------------------------
//------------------------------------------------------------------------------
// IN:  opr, opi - double precision vector of real and imaginary coefficients in order of decreasing powers.
//      degree - integer degree of polynomial
// OUT: zeror,zeroi - output double precision vectors of the real and imaginary parts of the zeros.
// RET: -1 if leading coefficient is zero, otherwise number of roots found. 
//------------------------------------------------------------------------------
int
Cpoly::Calculate( const double* opr, const double* opi, int degree, double* zeror, double* zeroi) 
{
   bool   bContinue ;
   bool   bConv ;
   int    cnt1, cnt2, idnn2, i ;
   double xx, yy, cosr, sinr, smalno, base, xxx, zr, zi, bnd ;


  // The following statements set machine constants.
   mcon( &eta, &infin, &smalno, &base) ;
   are = eta ;
   mre = 2.0 * sqrt( 2.0 ) * eta ;

  //  Initialization of constants for shift rotation.
   xx = 0.70710678 ;
   yy = -xx ;
   cosr = -0.060756474 ;
   sinr = -0.99756405 ;
   nn = degree ;  

  // Inizializzo numero iterazioni
   itercnt = 0 ;

  // Algorithm fails if the leading coefficient is zero
   if ( opr[0] == 0  &&  opi[0] == 0)
      return - 1 ;

  // Remove the zeros at the origin if any
   while ( opr[nn] == 0  &&  opi[nn] == 0) {
      idnn2 = degree - nn ;
      zeror[idnn2] = 0 ;
      zeroi[idnn2] = 0 ;
      nn -- ;
   }

  // Make a copy of the coefficients
   for ( i = 0 ; i <= nn ; i++) {
      pr[i] = opr[i] ;
      pi[i] = opi[i] ;
      shr[i] = cmod( pr[i], pi[i]) ;
   }

  // Scale the polynomial
   bnd = scale( nn, shr, eta, infin, smalno, base) ;
   if ( bnd != 1)
      for ( i = 0 ; i <= nn ; i++) {
         pr[i] *= bnd ;
         pi[i] *= bnd ;
      }

  // Main loop
   bContinue = true ;
   while ( bContinue) {   
 
      if ( nn <= 1) {
         cdivid( -pr[1], -pi[1], pr[0], pi[0], &zeror[degree-1], &zeroi[degree-1]) ;
         return degree ;       
      }

     // Calculate bnd, alower bound on the modulus of the zeros
      for ( i = 0 ; i <= nn ; i++)
         shr[i] = cmod( pr[i], pi[i]) ;

      cauchy( nn, shr, shi, &bnd) ;
   
     // Outer loop to control 2 Major passes with different sequences of shifts
      bContinue = false ;
      for ( cnt1 = 1 ; cnt1 <= 2  &&  ! bContinue ; cnt1++) {
        // First stage  calculation , no shift
         noshft( 5) ;

        // Inner loop to select a shift
         for ( cnt2 = 1 ; cnt2 <= 9  &&  ! bContinue ; cnt2++) {
           // Shift is chosen with modulus bnd and amplitude rotated by 94 degree from the previous shif
            xxx = cosr * xx - sinr * yy ;
            yy = sinr * xx + cosr * yy ;
            xx = xxx ;
            sr = bnd * xx ;
            si = bnd * yy ;

           // Second stage calculation, fixed shift
            fxshft( 10 * cnt2, &zr, &zi, &bConv) ;
            if ( bConv) {
              // The second stage jumps directly to the third stage ieration
              // If successful the zero is stored and the polynomial deflated
               idnn2 = degree - nn ;
               zeror[idnn2] = zr ;
               zeroi[idnn2] = zi ;
               nn -- ;
               for ( i = 0 ; i <= nn ; i++) {
                  pr[i] = qpr[i] ;
                  pi[i] = qpi[i] ;
               }
               bContinue = true ;
            }
           // If the iteration is unsuccessful another shift is chosen
         }
        // if 9 shifts fail, the outer loop is repeated with another sequence of shifts
      }
   }

  // The zerofinder has failed on two major passes
  // return empty handed with the number of roots found (less than the original degree)
   degree -= nn ;

   return degree ;       
}

//------------------------------------------------------------------------------
// COMPUTES  THE DERIVATIVE  POLYNOMIAL AS THE INITIAL H
// POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS.
//------------------------------------------------------------------------------
void
Cpoly::noshft( const int l1)
{
   int    i, j, jj, n, nm1 ;
   double xni, t1, t2 ;


   n = nn ;
   nm1 = n - 1 ;
   for ( i = 0 ; i < n ; i++) {
      xni = nn - i ;
      hr[i] = xni * pr[i] / n ;
      hi[i] = xni * pi[i] / n ;
   }
   for ( jj = 1 ; jj <= l1 ; jj++) {
      itercnt ++ ;
      if ( cmod( hr[n - 1], hi[n - 1]) > eta * 10 * cmod( pr[n - 1], pi[n - 1])) {
         cdivid( -pr[nn], -pi[nn], hr[n - 1], hi[n - 1], &tr, &ti) ;
         for ( i = 0 ; i < nm1 ; i++) {
            j = nn - i - 1 ;
            t1 = hr[j - 1] ;
            t2 = hi[j - 1] ;
            hr[j] = tr * t1 - ti * t2 + pr[j] ;
            hi[j] = tr * t2 + ti * t1 + pi[j] ;
         }
         hr[0] = pr[0] ;
         hi[0] = pi[0] ;
      }
      else {
         // If the constant term is essentially zero, shift H coefficients
         for ( i = 0 ; i < nm1 ; i++) {
            j = nn - i - 1 ;
            hr[j] = hr[j - 1] ;
            hi[j] = hi[j - 1] ;
         }
         hr[0] = 0 ;
         hi[0] = 0 ;
      }
   }
}

//------------------------------------------------------------------------------
// COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE.
// INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE
// APPROXIMATE ZERO IF SUCCESSFUL.
// L2 - LIMIT OF FIXED SHIFT STEPS
// ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE.
// CONV  - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION
//------------------------------------------------------------------------------
void
Cpoly::fxshft( const int l2, double* zr, double* zi, bool* pbConv)
{
   bool   bBol ;
   bool   bPasd ;
   bool   bTest ;
   int    i, j, n ;
   double otr, oti, svsr, svsi ;


   n = nn ;
   polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ;
   bTest = true ;
   bPasd = false ;

   // Calculate first T = -P(S)/H(S)
   calct( &bBol) ;

   // Main loop for second stage
   for ( j = 1 ; j <= l2 ; j++) {
      itercnt ++ ;

      otr = tr ;
      oti = ti ;

      // Compute the next H Polynomial and new t
      nexth( bBol) ;
      calct( &bBol) ;
      *zr = sr + tr ;
      *zi = si + ti ;

      // Test for convergence unless stage 3 has failed once or this
      // is the last H Polynomial
      if ( ! ( bBol  ||  ! bTest  ||  j == 12)) {
         if ( cmod( tr - otr, ti - oti) < 0.5 * cmod( *zr, *zi)) {
            if ( bPasd) {
               // The weak convergence test has been passwed twice, start the third stage
               // Iteration, after saving the current H polynomial and shift
               for ( i = 0; i < n; i++ ) {
                  shr[i] = hr[i] ;
                  shi[i] = hi[i] ;
               }
               svsr = sr ;
               svsi = si ;
               vrshft( 10, zr, zi, pbConv) ;
               if ( *pbConv)
                  return ;

               //The iteration failed to converge. Turn off testing and restore h,s,pv and T
               bTest = false ;
               for ( i = 0 ; i < n ; i++) {
                  hr[i] = shr[i] ;
                  hi[i] = shi[i] ;
               }
               sr = svsr ;
               si = svsi ;
               polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ;
               calct( &bBol) ;
               continue ;
            }
            bPasd = true ;
         }
         else
            bPasd = false ;
      }
   }

   // Attempt an iteration with final H polynomial from second stage
   vrshft( 10, zr, zi, pbConv) ;
}

//------------------------------------------------------------------------------
// CARRIES OUT THE THIRD STAGE ITERATION.
// L3 - LIMIT OF STEPS IN STAGE 3.
// ZR,ZI   - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE
//           ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT.
// CONV    -  .TRUE. IF ITERATION CONVERGES
//------------------------------------------------------------------------------
void
Cpoly::vrshft( const int l3, double* zr, double* zi, bool* pbConv)
{
   bool   bBol ;
   bool   bFlag ;
   bool   bNext ;
   int    i, j ;
   double mp, ms, omp, relstp, r1, r2, tp ;


   *pbConv = false ;
   bFlag = false ;
   sr = *zr ;
   si = *zi ;

  // Main loop for stage three
   for ( i = 1 ; i <= l3 ; i++) {
      itercnt ++ ;
     // Evaluate P at S and test for convergence
      polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ;
      mp = cmod( pvr, pvi) ;
      ms = cmod( sr, si) ;
      if ( mp <= 20 * errev( nn, qpr, qpi, ms, mp, are, mre)) {
        // Polynomial value is smaller in value than a bound on the error
        // in evaluationg P, terminate the iteration
         *pbConv = true ;
         *zr = sr ;
         *zi = si ;
         return ;
      }
      bNext = false ;
      if ( i != 1) {
         if ( ! ( bFlag  ||  mp < omp  ||  relstp >= 0.05)) {
           // Iteration has stalled. Probably a cluster of zeros. Do 5 fixed 
           // shift steps into the cluster to force one zero to dominate
            tp = relstp ;
            bFlag = true ;
            if ( relstp < eta)
               tp = eta ;
            r1 = sqrt( tp) ;
            r2 = sr * ( 1 + r1 ) - si * r1 ;
            si = sr * r1 + si * ( 1 + r1) ;
            sr = r2 ;
            polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ;
            for ( j = 1 ; j <= 5 ; j++) {
               calct( &bBol) ;
               nexth( bBol) ;
            }
            omp = infin ;
            bNext = true ;
         }
         
        // Exit if polynomial value increase significantly
         if ( ! bNext  &&  mp * 0.1 > omp)
            return ;
      }

     // eventuale salvataggio dato
      if ( ! bNext)
         omp = mp ;

     // Calculate next iterate
      calct( &bBol) ;
      nexth( bBol) ;
      calct( &bBol) ;
      if ( ! bBol) {
         relstp = cmod( tr, ti) / cmod( sr, si) ;
         sr += tr ;
         si += ti ;
      }
   }
}

//------------------------------------------------------------------------------
// COMPUTES  T = -P(S)/H(S).
// bool   - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO.
//------------------------------------------------------------------------------
void
Cpoly::calct( bool* pbBol)
{
   int    n ;
   double hvr, hvi ;


   n = nn ;

  // evaluate h(s)
   polyev( n - 1, sr, si, hr, hi, qhr, qhi, &hvr, &hvi) ;
   *pbBol = ( cmod( hvr, hvi ) <= are * 10 * cmod( hr[n - 1], hi[n - 1])) ;
   if ( ! *pbBol) {
      cdivid( -pvr, -pvi, hvr, hvi, &tr, &ti) ;
      return ;
   }

   tr = 0 ;
   ti = 0 ;
}

//------------------------------------------------------------------------------
// CALCULATES THE NEXT SHIFTED H POLYNOMIAL.
// bool   -  LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO
//------------------------------------------------------------------------------
void
Cpoly::nexth( bool bBol)
{
   int    j, n ;
   double t1, t2 ;


   n = nn ;
   if ( ! bBol) {
      for ( j = 1 ; j < n ; j++) {
         t1 = qhr[j - 1] ;
         t2 = qhi[j - 1] ;
         hr[j] = tr * t1 - ti * t2 + qpr[j] ;
         hi[j] = tr * t2 + ti * t1 + qpi[j] ;
      }
      hr[0] = qpr[0] ;
      hi[0] = qpi[0] ;
      return ;
   }

  // If h[s] is zero replace H with qh
   for ( j = 1 ; j < n ; j++) {
      hr[j] = qhr[j - 1] ;
      hi[j] = qhi[j - 1] ;
   }
   hr[0] = 0 ;
   hi[0] = 0 ;
}

//------------------------------------------------------------------------------
// EVALUATES A POLYNOMIAL  P  AT  S  BY THE HORNER RECURRENCE
// PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV.
//------------------------------------------------------------------------------
void
Cpoly::polyev( const int nn, const double sr, const double si, const double pr[], const double pi[],
               double qr[], double qi[], double *pvr, double *pvi )  
{
   int    i ;
   double t ;


   qr[0] = pr[0] ;
   qi[0] = pi[0] ;
   *pvr = qr[0] ;
   *pvi = qi[0] ;

   for ( i = 1 ; i <= nn ; i++) {
      t = ( *pvr) * sr - ( *pvi) * si + pr[i] ;
      *pvi = ( *pvr) * si + ( *pvi) * sr + pi[i] ;
      *pvr = t ;
      qr[i] = *pvr ;
      qi[i] = *pvi ;
   }
}

//------------------------------------------------------------------------------
// BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE.
// QR,QI - THE PARTIAL SUMS
// MS    -MODULUS OF THE POINT
// MP    -MODULUS OF POLYNOMIAL VALUE
// ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION
//------------------------------------------------------------------------------
double
Cpoly::errev( const int nn, const double qr[], const double qi[], const double ms, const double mp,
              const double are, const double mre )
{
   int    i ;
   double e ;


   e = cmod( qr[0], qi[0]) * mre / ( are + mre) ;
   for ( i = 0 ; i <= nn ; i++)
      e = e * ms + cmod( qr[i], qi[i]) ;

   return ( e * ( are + mre ) - mp * mre) ;
}

//------------------------------------------------------------------------------
// CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A
// POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS.
//------------------------------------------------------------------------------
void
Cpoly::cauchy( const int nn, double pt[], double q[], double* fn_val)
{
   int    i, n ;
   double x, xm, f, dx, df ;


   pt[nn] = -pt[nn] ;

  // Compute upper estimate bound
   n = nn ;
   x = exp( ( log( - pt[n]) - log( pt[0])) / (double) n) ;
   if ( pt[n - 1] != 0) {
      // Newton step at the origin is better, use it
      xm = -pt[nn] / pt[n - 1] ;
      if ( xm < x)
         x = xm ;
   }

  // Chop the interval (0,x) until f < 0
   while ( true) {
      xm = x * 0.1 ;
      f = pt[0] ;
      for ( i = 1 ; i <= nn ; i++)
         f = f * xm + pt[i] ;
      if ( f <= 0)
         break ;
      x = xm ;
   }
   dx = x ;
   
  // Do Newton iteration until x converges to two decimal places
   while ( abs( dx / x ) > 0.005) {
      q[0] = pt[0] ;
      for ( i = 1 ; i <= nn ; i++)
         q[i] = q[i - 1] * x + pt[i] ;
      f = q[nn] ;
      df = q[0] ;
      for ( i = 1 ; i < n ; i++)
         df = df * x + q[i] ;
      dx = f / df ;
      x -= dx ;
      itercnt ++ ;
   }

   *fn_val = x ;
}

//------------------------------------------------------------------------------
// RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE POLYNOMIAL.
// THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID UNDETECTED UNDERFLOW
// INTERFERING WITH THE CONVERGENCE CRITERION.  THE FACTOR IS A POWER OF THE BASE.
// PT - MODULUS OF COEFFICIENTS OF P
// ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC.
//------------------------------------------------------------------------------
double
Cpoly::scale( const int nn, const double pt[], const double eta,
              const double infin, const double smalno, const double base)
{
   int    i, l ;
   double hi, lo, max, min, x, sc ;
   double fn_val ;


   // Find largest and smallest moduli of coefficients
   hi = sqrt( infin) ;
   lo = smalno / eta ;
   max = 0 ;
   min = infin ;

   for ( i = 0 ; i <= nn ; i++) {
      x = pt[i] ;
      if ( x > max)
         max = x ;
      if ( x != 0  &&  x < min)
         min = x ;
   }

   // Scale only if there are very large or very small components
   fn_val = 1 ;
   if ( min >= lo  &&  max <= hi)
      return fn_val ;
   x = lo / min ;
   if ( x <= 1)
      sc = 1 / ( sqrt( max)* sqrt( min)) ;
   else {
      sc = x;
      if ( infin / sc > max)
         sc = 1 ;
   }
   l = (int)( log( sc) / log( base) + 0.5) ;
   fn_val = pow( base, l) ;
   return fn_val ;
}

//------------------------------------------------------------------------------
// COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW.
//------------------------------------------------------------------------------
void
Cpoly::cdivid( const double ar, const double ai, const double br, const double bi, double* cr, double* ci)
{
   double r, dinv, t, infin ;


   if ( br == 0  &&  bi == 0) {
     // Division by zero, c = infinity
      mcon( &t, &infin, &t, &t) ;
      *cr = infin ;
      *ci = infin ;
      return ;
   }

   if ( abs( br) < abs( bi)) {
      r = br / bi ;
      dinv = 1.0 / ( bi + r * br) ;
      *cr = ( ar * r + ai) * dinv ;
      *ci = ( ai * r - ar) * dinv ;
      return ;
   }

   r = bi / br ;
   dinv = 1.0 / ( br + r * bi) ;
   *cr = ( ar + ai * r) * dinv ;
   *ci = ( ai - ar * r) * dinv ;
}

//------------------------------------------------------------------------------
// MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW.
//------------------------------------------------------------------------------
double
Cpoly::cmod( const double r, const double i)
{
   double ar = abs( r) ;
   double ai = abs( i) ;
   if ( ar < ai)
      return ( ai * sqrt( 1.0 + ( ar * ar) / ( ai * ai))) ;

   if ( ar > ai)
      return ( ar * sqrt( 1.0 + ( ai * ai) / ( ar * ar))) ;

   return ( ar * sqrt( 2.0)) ;
}

//------------------------------------------------------------------------------
// MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM.
// THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO
// COMPUTE THEM. THE MEANING OF THE FOUR CONSTANTS ARE -
// ETA       THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED
//           AS THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
//           1.0_dp + ETA &gt; 1.0.
// INFINY    THE LARGEST FLOATING-POINT NUMBER
// SMALNO    THE SMALLEST POSITIVE FLOATING-POINT NUMBER
// BASE      THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED
//------------------------------------------------------------------------------
void
Cpoly::mcon( double* eta, double* infiny, double* smalno, double* base)
{
   *base = RADIX_DOUBLE ;
   *eta = EPS_DOUBLE ;
   *infiny = BIG_DOUBLE ;
   *smalno = SMALL_DOUBLE ;
}