# Program Listing for File PolynomialRoots-Jenkins-Traub.cc

↰ Return to documentation for file (`PolynomialRoots-Jenkins-Traub.cc`)

    //rpoly_ak1.cpp - Program for calculating the roots of a polynomial of real coefficients.
    //Written in Visual C++ 2005 Express Edition
    //14 July 2007
    //
    //The sub - routines listed below are translations of the FORTRAN routines included in RPOLY.FOR,
    //posted off the NETLIB site as TOMS / 493:
    //
    //http: //www.netlib.org / toms / 493
    //
    //TOMS / 493 is based on the Jenkins - Traub algorithm.
    //
    //To distinguish the routines posted below from others, an _ak1 suffix has been appended to them.
    //
    //Following is a list of the major changes made in the course of translating the TOMS / 493 routines
    //to the C++ versions posted below:
    //1) All global variables have been eliminated.
    //2) The "FAIL" parameter passed into RPOLY.FOR has been eliminated.
    //3) RPOLY.FOR solves polynomials of degree up to 100, but      does    not explicitly state this limit.
    //rpoly_ak1 explicitly states this limit; uses the macro name MAXDEGREE to specify this limit;

    //and does a check to ensure that the user input variable Degree is not greater than MAXDEGREE
    //(if it is, an error message is output and rpoly_ak1 terminates).If a user wishes to compute
    //roots of polynomials of degree greater than MAXDEGREE, using a macro name like MAXDEGREE provides
    //the simplest way of offering this capability.
    //4) All "GO TO" statements have been eliminated.
    //
    //A small main program is included also, to provide an example of how to use rpoly_ak1.In this
    //example, data is input from a file to eliminate the need for a user to type data in via
    //the console.
    //

    #ifdef __GNUC__
    #pragma GCC diagnostic ignored "-Wunused-function"
    #endif
    #ifdef __clang__
    #pragma clang diagnostic ignored "-Wglobal-constructors"
    #pragma clang diagnostic ignored "-Wvla-extension"
    #pragma clang diagnostic ignored "-Wvla"
    #pragma clang diagnostic ignored "-Wunused-function"
    #pragma clang diagnostic ignored "-Wdocumentation-unknown-command"
    #endif

    #include "PolynomialRoots.hh"

    #include <iostream>
    #include <fstream>
    #include <cctype>
    #include <cmath>
    #include <cfloat>
    #include <limits>

    #ifndef DOXYGEN_SHOULD_SKIP_THIS
    using namespace std;
    #endif

    namespace PolynomialRoots {

      #ifndef DOXYGEN_SHOULD_SKIP_THIS

      using std::abs;
      using std::pow;
      using std::frexp;

      #ifndef M_PI
      #define M_PI 3.14159265358979323846264338328
      #endif

      //static valueType const maxValue   = std::numeric_limits<valueType>::max();
      //static valueType const minValue   = std::numeric_limits<valueType>::min();
      static valueType const epsilon    = std::numeric_limits<valueType>::epsilon();
      static valueType const epsilon10  = 10*epsilon;
      static valueType const epsilon100 = 100*epsilon;

      static valueType const RADFAC = M_PI / 180; // Degrees - to - radians conversion factor = pi / 180
      //static valueType const lb2    = log(2.0); // Dummy variable to avoid re - calculating this value in loop below
      static valueType const cosr   = cos(94.0 * RADFAC); //= -0.069756474
      static valueType const sinr   = sin(94.0 * RADFAC); //= 0.99756405

      //============================================================================
      // Divides p by the quadratic x^2+u*x+v
      // placing the quotient in q and the remainder in a, b
      static
      void
      QuadraticSyntheticDivision(
        indexType       NN,
        valueType       u,
        valueType       v,
        valueType const p[],
        valueType       q[],
        valueType     & a,
        valueType     & b
      ) {
        q[0] = b = p[0];
        q[1] = a = p[1] - (b*u);
        for ( indexType i = 2; i < NN; ++i )
          { q[i] = p[i]-(a*u+b*v); b = a; a = q[i]; }
      }

      //============================================================================
      // This routine calculates scalar quantities used to compute the next
      // K polynomial and new estimates of the quadratic coefficients.
      // calcSC - integer variable set here indicating how the calculations
      // are normalized to avoid overflow.
      static
      indexType
      calcSC(
        indexType       N,
        valueType       a,
        valueType       b,
        valueType     & a1,
        valueType     & a3,
        valueType     & a7,
        valueType     & c,
        valueType     & d,
        valueType     & e,
        valueType     & f,
        valueType     & g,
        valueType     & h,
        valueType const K[],
        valueType       u,
        valueType       v,
        valueType       qk[]
      ) {

        indexType dumFlag = 3; // TYPE = 3 indicates the quadratic is almost a factor of K

        // Synthetic division of K by the quadratic 1, u, v
        QuadraticSyntheticDivision( N, u, v, K, qk, c, d );

        if ( std::abs(c) <= epsilon100 * std::abs(K[N-1]) &&
             std::abs(d) <= epsilon100 * std::abs(K[N-2]) ) return dumFlag;

        h = v*b;
        if ( std::abs(d) >= std::abs(c) ) {
          dumFlag = 2; // TYPE = 2 indicates that all formulas are divided by d
          e  = a / d;
          f  = c / d;
          g  = u*b;
          a3 = e*(g+a) + h*(b/d);
          a1 = f*b - a;
          a7 = h + (f+u) * a;
        } else {
          dumFlag = 1; // TYPE = 1 indicates that all formulas are divided by c;
          e  = a/c;
          f  = d/c;
          g  = e*u;
          a3 = e*a + (g+h/c)*b;
          a1 = b - (a*(d/c));
          a7 = g*d + h*f + a;
        }
        return dumFlag;
      }

      //============================================================================
      // Computes the next K polynomials using the scalars computed in calcSC_ak1
      static
      void
      nextK(
        indexType       N,
        indexType       tFlag,
        valueType       a,
        valueType       b,
        valueType       a1,
        valueType     & a3,
        valueType     & a7,
        valueType       K[],
        valueType const qk[],
        valueType const qp[]
      ) {
        if ( tFlag == 3 ) {
          // Use unscaled form of the recurrence
          K[1] = K[0] = 0;
          for ( indexType i = 2; i < N; ++i ) K[i] = qk[i-2];
        } else {
          valueType temp = tFlag == 1 ? b : a;
          if ( std::abs(a1) > epsilon10 * std::abs(temp) ) {
            // Use scaled form of the recurrence
            a7 /= a1;
            a3 /= a1;
            K[0] = qp[0];
            K[1] = qp[1]-(a7*qp[0]);
            for ( indexType i = 2; i < N; ++i )
              K[i] = qp[i] - (a7*qp[i-1]) + a3*qk[i-2];
          } else {
            // If a1 is nearly zero, then use a special form of the recurrence
            K[0] = 0;
            K[1] = -a7 * qp[0];
            for ( indexType i = 2; i < N; ++i )
              K[i] = a3*qk[i-2] - a7*qp[i-1];
          }
        }
      }

      //============================================================================
      // Compute new estimates of the quadratic coefficients
      // using the scalars computed in calcSC
      static
      void
      newest(
        indexType       tFlag,
        valueType     & uu,
        valueType     & vv,
        valueType       a,
        valueType       a1,
        valueType       a3,
        valueType       a7,
        valueType       b,
        valueType       c,
        valueType       d,
        valueType       f,
        valueType       g,
        valueType       h,
        valueType       u,
        valueType       v,
        valueType const K[],
        indexType       N,
        valueType const p[]
      ) {

        vv = uu = 0; // The quadratic is zeroed
        if ( tFlag != 3 ) {
          valueType  a4, a5;
          if (tFlag != 2) {
            a4 = a + u * b + h * f;
            a5 = c + (u + v * f) * d;
          } else {
            a4 = (a + g) * f + h;
            a5 = (f + u) * c + v * d;
          }
          // Evaluate new quadratic coefficients
          valueType b1 = -K[N - 1] / p[N];
          valueType b2 = -(K[N - 2] + b1 * p[N - 1]) / p[N];
          valueType c1 = v * b2 * a1;
          valueType c2 = b1 * a7;
          valueType c3 = b1 * b1 * a3;
          valueType c4 = -(c2 + c3) + c1;
          valueType temp = -c4 + a5 + b1 * a4;
          if ( temp != 0.0) {
            uu = -((u * (c3 + c2) + v * (b1 * a1 + b2 * a7)) / temp) + u;
            vv = v * (1.0 + c4 / temp);
          }
        }
      }

      //============================================================================
      // 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
      static
      void
      RealIT(
        indexType     & iFlag,
        indexType     & NZ,
        valueType     & sss,
        indexType       N,
        valueType const p[],
        indexType       NN,
        valueType       qp[],
        valueType     & szr,
        valueType     & szi,
        valueType       K[],
        valueType       qk[]
      ) {
        iFlag = NZ = 0;
        valueType t=0, omp=0, s = sss;
        for (indexType j=0;;) {
          valueType pv = p[0]; // Evaluate p at s
          qp[0] = pv;
          for ( indexType i = 1; i < NN; ++i )
            qp[i] = pv = pv * s + p[i];
          valueType mp = std::abs(pv);
          // Compute a rigorous bound on the error in evaluating p
          valueType ms = std::abs(s);
          valueType ee = 0.5 * std::abs(qp[0]);
          for ( indexType i = 1; i < NN; ++i ) ee = ee * ms + std::abs(qp[i]);
          // Iteration has converged sufficiently if the polynomial
          // value is less than 20 times this bound
          if ( mp <= 20.0 * epsilon * (2*ee - mp) )
            { NZ = 1; szr = s; szi = 0; break; }

          // Stop iteration after 10 steps
          if ( ++j > 10 ) break;
          if ( j >= 2 ) {
            if ( (abs(t) <= 0.001 * std::abs(s-t)) && (mp > omp) ) {
              // A cluster of zeros near the real axis has been encountered;
              // Return with iFlag set to initiate a quadratic iteration
              iFlag = 1;
              sss   = s;
              break;
            }
          }
          // Return if the polynomial value has increased significantly
          omp = mp;
          // Compute t, the next polynomial and the new iterate
          valueType kv = qk[0] = K[0];
          for ( indexType i = 1; i < N; ++i )
            qk[i] = kv = kv * s + K[i];

          if ( std::abs(kv) > std::abs(K[N-1]) * epsilon10 ) {
            // Use the scaled form of the recurrence if the value of K at s is non - zero
            valueType tt = -(pv / kv);
            K[0] = qp[0];
            for ( indexType i = 1; i < N; ++i ) K[i] = tt * qk[i-1] + qp[i];
          } else { // Use unscaled form
            K[0] = 0.0;
            for ( indexType i = 1; i < N; ++i ) K[i] = qk[i-1];
          }
          kv = K[0];
          for ( indexType i = 1; i < N; ++i ) kv = kv * s + K[i];
          t = std::abs(kv) > std::abs(K[N-1])*epsilon10 ? -(pv/kv) : 0;
          s += t;
        }
      }

      //============================================================================
      // Variable - shift K - polynomial iteration for a quadratic
      // factor converges only if the zeros are equimodular or nearly so.
      static
      void
      QuadIT(
        indexType       N,
        indexType     & NZ,
        valueType       uu,
        valueType       vv,
        valueType     & szr,
        valueType     & szi,
        valueType     & lzr,
        valueType     & lzi,
        valueType       qp[],
        indexType       NN,
        valueType     & a,
        valueType     & b,
        valueType const p[],
        valueType       qk[],
        valueType     & a1,
        valueType     & a3,
        valueType     & a7,
        valueType     & c,
        valueType     & d,
        valueType     & e,
        valueType     & f,
        valueType     & g,
        valueType     & h,
        valueType       K[]
      ) {

        NZ = 0; // Number of zeros found
        indexType   j = 0, tFlag;
        valueType u = uu; // uu and vv are coefficients of the starting quadratic
        valueType v = vv;
        valueType relstp = 0; // initialize remove warning
        valueType omp    = 0; // initialize remove warning
        valueType ui, vi;
        bool triedFlag = false;
        do {
          Quadratic solve( 1.0, u, v );
          solve.getRoot0(szr,szi);
          solve.getRoot1(lzr,lzi);

          // Return if roots of the quadratic are real and not close
          // to multiple or nearly equal and of opposite sign.
          if ( std::abs(abs(szr) - std::abs(lzr)) > 0.01 * std::abs(lzr) ) break;
          // Evaluate polynomial by quadratic synthetic division
          QuadraticSyntheticDivision( NN, u, v, p, qp, a, b );
          valueType mp = std::abs(a-(szr*b)) + std::abs(szi*b);
          // Compute a rigorous bound on the rounding error in evaluating p
          valueType zm = std::sqrt(abs(v));
          valueType ee = 2 * std::abs(qp[0]);
          valueType t  = -szr*b;
          for ( indexType i = 1; i < N; ++i ) ee = ee * zm + std::abs(qp[i]);
          ee = ee * zm + std::abs(a+t);
          ee = (9*ee+2*abs(t)-7*(abs(a+t)+zm*abs(b)))*epsilon;
          // Iteration has converged sufficiently if the polynomial
          // value is less than 20 times this bound
          if ( mp <= 20*ee ) { NZ = 2; break; }
          // Stop iteration after 20 steps
          if ( ++j > 20 ) break;
          if ( j >= 2 ) {
            if ( (relstp <= 0.01) && (mp >= omp) && !triedFlag ) {
              // A cluster appears to be stalling the convergence.
              // Five fixed shift steps are taken with a u, v close to the cluster.
              relstp = (relstp < epsilon) ? std::sqrt(epsilon) : std::sqrt(relstp);
              u -= u * relstp;
              v += v * relstp;
              QuadraticSyntheticDivision(NN, u, v, p, qp, a, b);
              for ( indexType i = 0; i < 5; ++i ) {
                tFlag = calcSC(N, a, b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
                nextK(N, tFlag, a, b, a1, a3, a7, K, qk, qp);
              }
              triedFlag = true;
              j = 0;
            }
          }
          omp = mp;
          // Calculate next K polynomial and new u and v
          tFlag = calcSC(N, a, b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
          nextK(N, tFlag, a, b, a1, a3, a7, K, qk, qp);
          tFlag = calcSC(N, a, b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
          newest(tFlag, ui, vi, a, a1, a3, a7, b, c, d, f, g, h, u, v, K, N, p);

          // If vi is zero, the iteration is not converging
          if ( !isZero(vi) ) {
            relstp = std::abs((vi-v)/vi);
            u = ui;
            v = vi;
          }
        } while( !isZero(vi) );
      }

      //============================================================================
      /*
      // 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.
      //
      // L2 limit of fixed shift steps
      // NZ number of zeros found
      */
      static
      indexType
      FixedShift(
        indexType   L2,
        valueType   sr,
        valueType   v,
        valueType   K[],
        indexType   N,
        valueType   p[],
        indexType   NN,
        valueType   qp[],
        valueType   u,
        valueType & lzi,
        valueType & lzr,
        valueType & szi,
        valueType & szr
      ) {

        #ifdef _MSC_VER
        valueType * qk  = (valueType*)alloca( 2*(N+1)*sizeof(valueType) );
          valueType * svk = qk+N+1;
        #else
        valueType   qk[N+1], svk[N+1];
          #endif

        indexType iFlag = 1;
        indexType NZ    = 0;
        valueType betav = 0.25;
        valueType betas = 0.25;
        valueType oss   = sr;
        valueType ovv   = v;

        // Evaluate polynomial by synthetic division
        valueType a, b;
        QuadraticSyntheticDivision(NN, u, v, p, qp, a, b);
        valueType a1, a3, a7, c, d, e, f, g, h;
        indexType tFlag =   calcSC(N, a, b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);

        valueType otv = 0; // initialize remove warning
        valueType ots = 0; // initialize remove warning
        for ( indexType j = 0; j < L2; ++j ) {
          indexType fflag = 1;
          // Calculate next K polynomial and estimate v
          nextK( N, tFlag, a, b, a1, a3, a7, K, qk, qp );
          tFlag = calcSC(N, a, b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
          valueType ui, vi;
          newest(tFlag, ui, vi, a, a1, a3, a7, b, c, d, f, g, h, u, v, K, N, p);
          valueType vv = vi;
          // Estimate s
          valueType ss = (K[N-1] != 0.0) ? -(p[N]/K[N-1]) : 0.0;
          valueType ts = 1.0;
          valueType tv = 1.0;
          if ( (j != 0) && (tFlag != 3) ) {
            // Compute relative measures of convergence of s and v sequences
            tv = (vv != 0.0) ? std::abs((vv - ovv) / vv) : tv;
            ts = (ss != 0.0) ? std::abs((ss - oss) / ss) : ts;
            // If decreasing, multiply the two most recent convergence measures
            valueType tvv = (tv < otv) ? tv * otv : 1;
            valueType tss = (ts < ots) ? ts * ots : 1;
            // Compare with convergence criteria
            bool vpass = tvv < betav;
            bool spass = tss < betas;
            if ( spass || vpass ) {
              // At least one sequence has passed the convergence test.
              // Store variables before iterating
              for ( indexType i = 0; i < N; ++i ) svk[i] = K[i];
              valueType s = ss;
              // Choose iteration according to the fastest converging sequence
              indexType stry = 0;
              indexType vtry = 0;
              for (;;) {
                if ( (fflag && ((fflag = 0) == 0)) &&
                     ((spass) && (!vpass || (tss < tvv))) ) {
                          ;
                            //Do nothing.Provides a quick "short circuit".
                } else {
                  QuadIT(N, NZ, ui, vi, szr, szi, lzr, lzi, qp, NN, a, b, p, qk, a1, a3, a7, c, d, e, f, g, h, K);
                  if ( NZ > 0) return NZ;
                  // Quadratic iteration has failed.Flag that it has been
                  // tried and decrease the convergence criterion
                  iFlag = vtry = 1;
                  betav *= 0.25;
                  // Try linear iteration if it has not been tried and the s sequence is converging
                  if ( stry || !spass ) iFlag = 0;
                  else                  std::copy( svk, svk + N, K );
                }
                if ( iFlag != 0 ) {
                  RealIT(iFlag, NZ, s, N, p, NN, qp, szr, szi, K, qk);
                  if ( NZ > 0 ) return NZ;
                  // Linear iteration has failed.Flag that it has been
                  // tried and decrease the convergence criterion
                  stry = 1;
                  betas *= 0.25;
                  if (iFlag != 0) {
                    // If linear iteration signals an almost double real zero,                                attempt   quadratic iteration
                    ui = -(s + s);
                    vi = s * s;
                    continue;
                  }
                }
                // Restore variables
                std::copy( svk, svk+N, K );

                // Try quadratic iteration if it has not been tried
                // and the v sequence is converging
                if (!vpass || vtry) break; // Break out of infinite for loop
              }
              // Re - compute qp and scalar values to continue the second stage
              QuadraticSyntheticDivision(NN, u, v, p, qp, a, b);
              tFlag = calcSC(N, a, b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
            }
          }
          ovv = vv;
          oss = ss;
          otv = tv;
          ots = ts;
        }
        return NZ;
      }

      //============================================================================

      static
      inline
      valueType
      evalPoly( valueType x, valueType const p[], indexType Degree ) {
        valueType ff = p[0];
        for ( indexType i = 1; i <= Degree; ++i ) ff = ff * x + p[i];
        return ff;
      }

      //============================================================================

      static
      inline
      void
      evalPoly(
        valueType       x,
        valueType const p[],
        indexType       Degree,
        valueType     & f,
        valueType     & df
      ) {
        df = f = p[0];
        for ( indexType i = 1; i < Degree; ++i ) {
          f  = x * f + p[i];
          df = x * df + f;
        }
        f = x * f + p[Degree];
      }

      //============================================================================

      /*
      // Scale if there are large or very small coefficients
      // Computes 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.
      */
      static
      inline
      void
      scalePoly( valueType p[], indexType N ) {
        /*
        .. double ldexp(double x, int n)
        .. The ldexp() functions multiply x by 2 to the power n.
        ..
        .. double frexp(double value, int *exp);
        .. The frexp() functions break the floating-point number value into
        .. a normalized fraction and an integral power of 2.
        .. They store the integer in the int object pointed to by exp.
        .. The functions return a number x such that x has a magnitude in
        .. the interval [1/2, 1) or 0, and value = x*(2**exp).
        */
        int max_exponent = std::numeric_limits<int>::min();
        for ( int i = 0; i <= N; ++i ) {
          if ( !isZero(p[i]) ) {
            int exponent;
            frexp( p[i], &exponent );
            if ( exponent > max_exponent ) max_exponent = exponent;
          }
        }
        int l = -max_exponent;
        for ( indexType i = 0; i <= N; ++i ) p[i] = ldexp(p[i],l);
      }

      //============================================================================

      // Compute lower bound on moduli of zeros
      static
      inline
      valueType
      lowerBoundZeroPoly( valueType p[], indexType N ) {

        #ifdef _MSC_VER
        valueType * pt  = (valueType*)alloca( (N+1)*sizeof(valueType) );
        #else
        valueType pt[N+1];
          #endif

        for ( indexType i = 0; i < N; ++i ) pt[i] = std::abs(p[i]);
        pt[N] = -abs(p[N]);

        // Compute upper estimate of bound
        valueType x = exp((log(-pt[N]) - log(pt[0]))/N);
        if ( !isZero(pt[N-1]) ) { // If Newton step at the origin is better, use it
          valueType xm = -pt[N]/pt[N-1];
          if ( xm < x ) x = xm;
        }
        // Chop the interval(0, x) until f <= 0
        valueType xm = x;
        while ( evalPoly( xm, pt, N ) > 0 ) { x = xm; xm = 0.1 * x; }

        // Do Newton iteration until x converges to two decimal places
        valueType dx;
        do {
          valueType f, df;
          evalPoly( x, pt, N, f, df );
          dx = f / df;
          x -= dx;
        } while ( std::abs(dx) > std::abs(x)*0.005 );
        return x;
      }

      //============================================================================

      static
      inline
      void
      roots3(
        valueType const p[4],
        indexType       Degree,
        valueType       zeror[],
        valueType       zeroi[]
      ) {
        if ( Degree == 1 ) {
          zeror[0] = -(p[1]/p[0]);
          zeroi[0] = 0;
        } else if ( Degree == 2 ) {
          Quadratic solve( p[0], p[1], p[2] );
          switch ( solve.numRoots() ) {
            case 2: solve.getRoot1( zeror[1], zeroi[1] );
            case 1: solve.getRoot0( zeror[0], zeroi[0] );
          }
        } else if ( Degree == 3 ) {
          Cubic solve( p[0], p[1], p[2], p[3] );
          switch ( solve.numRoots() ) {
            case 3: solve.getRoot2( zeror[2], zeroi[2] );
            case 2: solve.getRoot1( zeror[1], zeroi[1] );
            case 1: solve.getRoot0( zeror[0], zeroi[0] );
          }
        }
      }

      #endif

      //============================================================================
      int
      roots(
        valueType const * op,
        indexType         Degree,
        valueType       * zeror,
        valueType       * zeroi
      ) {

        if ( Degree < 1 ) return -1;

        // Do a quick check to see if leading coefficient is 0
        // The leading coefficient is zero. No further action taken. Program terminated
        if ( isZero(op[0]) ) return -2;

        #ifdef _MSC_VER
        valueType * ptr = (valueType*)alloca( 4*(Degree+1)*sizeof(valueType) );
        valueType * K    = ptr; ptr += Degree+1;
        valueType * p    = ptr; ptr += Degree+1;
        valueType * qp   = ptr; ptr += Degree+1;
        valueType * temp = ptr; ptr += Degree+1;
        #else
        valueType K[Degree+1];
        valueType p[Degree+1];
        valueType qp[Degree+1];
        valueType temp[Degree+1];
        #endif

        int N = Degree;
        valueType xx = std::sqrt(0.5); //= 0.70710678
        valueType yy = -xx;

        // Remove zeros at the origin, if any
        for ( indexType j = 0; isZero(op[N]); ++j, --N ) zeror[j] = zeroi[j] = 0.0;
        std::copy( op, op+N+1, p ); // Make a copy of the coefficients

        while ( N > 0 ) {
          // Main loop
          // Start the algorithm for one zero
          if ( N < 4 ) { roots3( p, N, zeror+Degree-N, zeroi+Degree-N ); break; }

          // Find the largest and smallest moduli of the coefficients
          scalePoly( p, N );

          // Compute lower bound on moduli of zeros
          valueType bnd = lowerBoundZeroPoly( p, N );

          // Compute the derivative as the initial K polynomial and
          // do 5 steps with no shift
          for ( indexType i = 1; i < N; ++i ) K[i] = ((N-i) * p[i]) / N;
          K[0] = p[0];
          indexType NM1 = N-1;
          valueType aa = p[N];
          valueType bb = p[NM1];
          bool zerok = isZero(K[NM1]);
          for ( indexType iter = 0; iter < 5; ++iter ) {
            if ( zerok ) { // Use unscaled form of recurrence
              for ( indexType i = 0; i < NM1; ++i ) K[NM1-i] = K[NM1-i-1];
              K[0] = 0;
              zerok = isZero(K[NM1]);
            } else { // Used scaled form of recurrence if value of K at 0 is nonzero
              valueType t = -aa / K[NM1];
              for ( indexType i = 0; i < NM1; ++i ) {
                indexType j = NM1-i;
                K[j] = t * K[j-1] + p[j];
              }
              K[0] = p[0];
              zerok = std::abs(K[NM1]) <= std::abs(bb) * epsilon10;
            }
          }

          // Save K for restarts with new shifts
          std::copy( K, K+N, temp );

          // Loop to select the quadratic corresponding to each new shift
          bool ok = false;
          for ( indexType iter = 0; iter < 20 && !ok; ++iter ) {
            // 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.
            valueType tmp = -(sinr * yy) + cosr * xx;
            yy = sinr * xx + cosr * yy;
            xx = tmp;
            valueType sr = bnd * xx;
            valueType u = -2*sr;
            // Second stage calculation, fixed quadratic
            valueType lzi, lzr, szi, szr;
            indexType NZ = FixedShift( 20*(iter+1), sr, bnd, K, N, p, N+1, qp, u, lzi, lzr, szi, szr);
            ok = NZ != 0;
            if ( ok ) {
              //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.
              indexType j = Degree - N;
              zeror[j] = szr;
              zeroi[j] = szi;
              N -= NZ;
              for ( indexType i = 0; i <= N; i++) p[i] = qp[i];
              if ( NZ != 1 ) {
                zeror[j+1] = lzr;
                zeroi[j+1] = lzi;
              }
            } else {
              // If the iteration is unsuccessful, another quadratic is chosen after restoring K
              std::copy( temp, temp+N, K );
            }
          }
          // Return with failure if no convergence with 20 shifts
          if ( !ok ) return -2;
        }
        return 0;
      }
    }