/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2017                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Universita` degli Studi di Trento                                   |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/

#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 <cmath>
#include <iostream>
#include <algorithm>
#include <limits>

#ifndef DOXYGEN_SHOULD_SKIP_THIS

#define MAX_ITER_SAFETY 50

namespace PolynomialRoots {

  using std::abs;
  static valueType const machepsi = std::numeric_limits<valueType>::epsilon();

  indexType
  Quartic::getRealRoots( valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx0() ) r[nr++] = r0;
    if ( !cplx1() ) r[nr++] = r1;
    if ( !cplx2() ) r[nr++] = r2;
    if ( !cplx3() ) r[nr++] = r3;
    return nr;
  }

  indexType
  Quartic::getPositiveRoots( valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx0() && r0 > 0 ) r[nr++] = r0;
    if ( !cplx1() && r1 > 0 ) r[nr++] = r1;
    if ( !cplx2() && r2 > 0 ) r[nr++] = r2;
    if ( !cplx3() && r3 > 0 ) r[nr++] = r3;
    return nr;
  }

  indexType
  Quartic::getNegativeRoots( valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx0() && r0 < 0 ) r[nr++] = r0;
    if ( !cplx1() && r1 < 0 ) r[nr++] = r1;
    if ( !cplx2() && r2 < 0 ) r[nr++] = r2;
    if ( !cplx3() && r3 < 0 ) r[nr++] = r3;
    return nr;
  }

  indexType
  Quartic::getRootsInRange( valueType a, valueType b, valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx0() && r0 >= a && r0 <= b ) r[nr++] = r0;
    if ( !cplx1() && r1 >= a && r1 <= b ) r[nr++] = r1;
    if ( !cplx2() && r2 >= a && r2 <= b ) r[nr++] = r2;
    if ( !cplx3() && r3 >= a && r3 <= b ) r[nr++] = r3;
    return nr;
  }

  indexType
  Quartic::getRootsInOpenRange( valueType a, valueType b, valueType r[] ) const {
    indexType nr = 0;
    if ( !cplx0() && r0 > a && r0 < b ) r[nr++] = r0;
    if ( !cplx1() && r1 > a && r1 < b ) r[nr++] = r1;
    if ( !cplx2() && r2 > a && r2 < b ) r[nr++] = r2;
    if ( !cplx3() && r3 > a && r3 < b ) r[nr++] = r3;
    return nr;
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  //
  // x^4 + A x^3 + B x^2 + C x + D
  static
  inline
  void
  scaleQuarticMonicPolynomial(
    valueType   A,
    valueType   B,
    valueType   C,
    valueType   D,
    valueType & AS,
    valueType & BS,
    valueType & CS,
    valueType & DS,
    indexType & i_case,
    valueType & scale
  ) {

    valueType a = std::abs(A);
    valueType b = std::sqrt(abs(B));
    valueType c = std::cbrt(abs(C));
    valueType d = std::sqrt(sqrt(abs(D)));

    if ( a < b ) {
      if ( b < c ) {
        if ( c < d ) i_case = 0; // a < b < c < d        --> d MAX
        else         i_case = 1; // a < b < c and d <= c --> c MAX
      } else {
        if ( b < d ) i_case = 0; // a < b and c <= b < d --> d MAX
        else         i_case = 2; // a < b and c <= b and d <= b --> b MAX
      }
    } else {
      if ( a < c ) {
        if ( c < d ) i_case = 0; // b <= a < c < d        --> d MAX
        else         i_case = 1; // b <= a < c and d <= c --> c MAX
      } else {
        if ( a < d ) i_case = 0; // b <= a and c <= a and a < d --> d MAX
        else         i_case = 3; // b <= a and c <= a and d <= a --> a MAX
      }
    }

    switch ( i_case ) {
      case 0:
        scale = d;
        AS    = A/d;
        BS    = (B/d)/d;
        CS    = ((C/d)/d)/d;
        DS    = D > 0 ? 1 : -1;
      break;
      case 1:
        scale = c;
        AS    = A/c;
        BS    = (B/c)/c;
        CS    = C > 0 ? 1 : -1;
        DS    = (((D/c)/c)/c)/c;
      break;
      case 2:
        scale = b;
        AS    = A/b;
        BS    = B > 0 ? 1 : -1;
        CS    = ((C/b)/b)/b;
        DS    = (((D/b)/b)/b)/b;
      break;
      case 3:
        scale = a;
        AS    = A > 0 ? 1 : -1;
        BS    = (B/a)/a;
        CS    = ((C/a)/a)/a;
        DS    = (((D/a)/a)/a)/a;
      break;
    }
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  static
  inline
  valueType
  evalHexic(
    valueType x,
    valueType q3,
    valueType q2,
    valueType q1,
    valueType q0
  ) {
    valueType t1 = x + q3;
    valueType t2 = x + t1;
    valueType t3 = x + t2;
    t1 = t1 * x + q2;
    t2 = t2 * x + t1;
    t1 = t1 * x + q1;
    valueType Q    = t1 * x + q0; // Q(x)
    valueType dQ   = t2 * x + t1; // Q'(x)
    valueType ddQ  = t3 * x + t2; // Q''(x) / 2
    valueType dddQ = x + t3;      // Q'''(x) / 6
    return  Q * dddQ * dddQ - dQ * ddQ * dddQ + dQ * dQ; // H(x), usually < 0
  }

  static
  inline
  void
  evalHexic(
    valueType   x,
    valueType   q3,
    valueType   q2,
    valueType   q1,
    valueType   q0,
    valueType & p,
    valueType & dp
  ) {
    valueType t1 = x + q3;
    valueType t2 = x + t1;
    valueType t3 = x + t2;
    t1 = t1 * x + q2;
    t2 = t2 * x + t1;
    t1 = t1 * x + q1;
    valueType Q    = t1 * x + q0; // Q(x)
    valueType dQ   = t2 * x + t1; // Q'(x)
    valueType ddQ  = t3 * x + t2; // Q''(x) / 2
    valueType dddQ = x + t3;      // Q'''(x) / 6
    p  = Q * dddQ * dddQ - dQ * ddQ * dddQ + dQ * dQ; // H(x), usually < 0
    dp = 2 * dddQ * (4 * Q - dQ * dddQ - ddQ * ddQ); // H'(x)
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  //
  // a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = (x-r)*(a4*x^3+b2*x^2+b1*x+b0)
  static
  void
  deflateQuarticPolynomial(
    valueType   a4,
    valueType   a3,
    valueType   a2,
    valueType   a1,
    valueType   a0,
    valueType   r,
    valueType & b2,
    valueType & b1,
    valueType & b0
  ) {
    indexType i_cross  = 0;
    valueType r2       = r*r;
    valueType v_cross  = std::abs(a0);
    valueType v_cross1 = std::abs(a1*r);
    if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 1; }
    v_cross1 = std::abs(a2*r2);
    if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 2; }
    v_cross1 = std::abs(a3*r*r2);
    if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 3; }
    v_cross1 = std::abs(a4*r2*r2);
    if ( v_cross1 > v_cross ) i_cross = 4;
    switch ( i_cross ) {
      case 0: b2 = a3+a4*r; b1 = a2+r*b2; b0 = a1+r*b1;   break;
      case 1: b2 = a3+a4*r; b1 = a2+r*b2; b0 = -a0/r;     break;
      case 2: b2 = a3+a4*r; b0 = -a0/r;   b1 = (b0-a1)/r; break;
      case 3:
      case 4: b0 = -a0/r; b1 = (b0-a1)/r; b2 = (b1-a2)/r; break;
    }
  }

  /*
  ||   _   _               _              ____  _               _   _
  ||  | \ | | _____      _| |_ ___  _ __ | __ )(_)___  ___  ___| |_(_) ___  _ __
  ||  |  \| |/ _ \ \ /\ / / __/ _ \| '_ \|  _ \| / __|/ _ \/ __| __| |/ _ \| '_ \
  ||  | |\  |  __/\ V  V /| || (_) | | | | |_) | \__ \  __/ (__| |_| | (_) | | | |
  ||  |_| \_|\___| \_/\_/  \__\___/|_| |_|____/|_|___/\___|\___|\__|_|\___/|_| |_|
  */

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Translate to C from Polynomial234RootSolvers
  static
  indexType
  zeroQuarticByNewtonBisection(
    valueType  a,
    valueType  b,
    valueType  c,
    valueType  d,
    valueType & x
  ) {

    valueType p, dp;
    evalMonicQuartic( x, a, b, c, d, p, dp );
    valueType t = p; // save p(x) for sign comparison
    x -= p/dp; // 1st improved root

    indexType iter       = 1;
    indexType oscillate  = 0;
    indexType nconverged = 0;
    bool      bisection  = false;
    bool      converged  = false;
    valueType s(0), u(0); // to mute warning
    while ( ! ( nconverged > 1 || bisection ) && iter < MAX_ITER_SAFETY ) {
      ++iter;
      valueType ddp;
      evalMonicQuartic( x, a, b, c, d, p, dp, ddp );
      if ( p*t < 0 ) { // does Newton start oscillating ?
        if ( p < 0 ) {
          ++oscillate; // increment oscillation counter
          s = x;       // save lower bisection bound
        } else {
          u = x; // save upper bisection bound
        }
        t = p; // save current p(x)
      }
      //dp = (p*dp)/(dp*dp-p*ddp); // double zero correction
      dp = (p*dp)/(dp*dp-0.5*p*ddp); // Halley correction
      //dp = p/dp; // Newton correction
      x -= dp; // new Newton root
      bisection = oscillate > 2; // activate bisection
      converged = std::abs(dp) <= (1+std::abs(x)) * machepsi; // Newton convergence indicator
      if ( converged ) ++nconverged; else nconverged = 0;
    }
    if ( bisection || !converged ) {
      t = u - s; // initial bisection interval
      while ( std::abs(t) > (1+std::abs(x)) * machepsi && iter < MAX_ITER_SAFETY ) { // bisection iterates
        ++iter;
        p = evalMonicQuartic( x, a, b, c, d );
        if ( p < 0 ) s = x;
        else         u = x; // keep bracket on root
        t = (u-s)/2; // new bisection interval
        x = s + t;   // new bisection root
      }
    }
    return iter;
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Translate to C from Polynomial234RootSolvers
  static
  indexType
  zeroHexicByNewtonBisection(
    valueType   q3,
    valueType   q2,
    valueType   q1,
    valueType   q0,
    valueType & x
  ) {
    valueType p, dp;
    evalHexic( x, q3, q2, q1, q0, p, dp );
    valueType t = p; // save p(x) for sign comparison
    x -= p/dp; // 1st improved root

    indexType iter      = 1;
    indexType oscillate = 0;
    bool      bisection = false;
    bool      converged = false;
    valueType s(0), u(0); // to mute warning
    while ( ! (converged||bisection) ) {
      ++iter;
      evalHexic( x, q3, q2, q1, q0, p, dp );
      if ( p*t < 0 ) { // does Newton start oscillating ?
        if ( p < 0 ) {
          ++oscillate; // increment oscillation counter
          s = x;       // save lower bisection bound
        } else {
          u = x; // save upper bisection bound
        }
        t = p; // save current p(x)
      }
      dp = p/dp; // Newton correction
      x -= dp; // new Newton root
      bisection = oscillate > 2; // activate bisection
      converged = std::abs(dp) <= std::abs(x) * machepsi; // Newton convergence indicator
    }
    if ( bisection ) {
      t = u - s; // initial bisection interval
      while ( std::abs(t) > std::abs(x) * machepsi ) { // bisection iterates
        ++iter;
        p = evalHexic( x, q3, q2, q1, q0 );
        if ( p < 0 ) s = x;
        else         u = x; // keep bracket on root
        t = (u-s)/2; // new bisection interval
        x = s + t;   // new bisection root
      }
    }
    return iter;
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  /*
  ||       _       __ _       _
  ||    __| | ___ / _| | __ _| |_ ___
  ||   / _` |/ _ \ |_| |/ _` | __/ _ \
  ||  | (_| |  __/  _| | (_| | ||  __/
  ||   \__,_|\___|_| |_|\__,_|\__\___|
  */
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  //
  // a3*x^3 + a2*x^2 + a1*x + a0 = (x-r)*(a3*x^2+b1*x+b0)


  /*\
   *  Calculate the zeros of the quartic a*z^4 + b*z^3 + c*z^2 + d*x + e.
  \*/

  void
  Quartic::findRoots() {
    valueType const & A = ABCDE[0];
    valueType const & B = ABCDE[1];
    valueType const & C = ABCDE[2];
    valueType const & D = ABCDE[3];
    valueType const & E = ABCDE[4];

    iter = nreal = ncplx = 0;

    // special cases
    if ( isZero(A) ) {
      Cubic csolve( B, C, D, E );
      nreal = csolve.numRoots();
      switch ( nreal ) {
        case 3: r2 = csolve.real_root2();
        case 2: r1 = csolve.real_root1();
        case 1: r0 = csolve.real_root0();
      }
      if ( csolve.complexRoots() ) { ncplx = 2; nreal -= 2; }
      return;
    }
    if ( isZero(E) ) {
      Cubic csolve( A, B, C, D );
      nreal = csolve.numRoots();
      r0 = r1 = r2 = r3 = 0;
      switch ( nreal ) {
        case 3: r2 = csolve.real_root2();
        case 2: r1 = csolve.real_root1();
        case 1: r0 = csolve.real_root0();
      }
      ++nreal;
      if ( csolve.complexRoots() ) { ncplx = 2; nreal -= 2; }
      if ( nreal == 4 ) { // caso regolare, 4 radici reali maintain order
        if ( r3 < r2 ) std::swap(r2,r3);
        if ( r2 < r1 ) std::swap(r1,r2);
        if ( r1 < r0 ) std::swap(r0,r1);
      }
      return;
    }
    if ( isZero(B) && isZero(D) ) { // biquadratic case
      // A x^4 + C x^2 + E
      Quadratic qsolve( A, C, E ) ;
      valueType x = qsolve.real_root0() ;
      valueType y = qsolve.real_root1() ;
      if ( qsolve.complexRoots() ) {
        // complex conjugate pair biquadratic roots x +/- iy.
        ncplx = 4;
        x /= 2; // re
        y /= 2; // im
        valueType z = hypot(x,y);
        y = std::sqrt(z - x);
        x = std::sqrt(z + x);
        r0 = -x;
        r1 = y;
        r2 = x;
        r3 = y;
      } else {
        // real roots of quadratic are ordered x <= y
        if ( x >= 0 ) { // y >= 0
          nreal = 4;
          x = std::sqrt(x); y = std::sqrt(y);
          r0 = -y; r1 = -x; r2 =  x; r3 =  y;
        } else if ( y >= 0 ) { // x < 0 && y >= 0
          nreal = ncplx = 2;
          x = std::sqrt(-x); y = std::sqrt(y);
          r0 =  0; r1 = x; // (real,imaginary)
          r2 = -y; r3 = y;
        } else { // x < 0 && y < 0
          ncplx = 4;
          x = std::sqrt(-x); y = std::sqrt(-y);
          r0 = 0; r1 = x; r2 = 0; r3 = y; // 2 x (real,imaginary)
        }
      }
      return;
    }

    /*
    .. The general case. Rescale quartic polynomial, such that largest absolute
    .. coefficient is (exactly!) equal to 1. Honor the presence of a special
    .. quartic case that might have been obtained during the rescaling process
    .. (due to underflow in the coefficients).
    */

    valueType const A3 = B/A;
    valueType const A2 = C/A;
    valueType const A1 = D/A;
    valueType const A0 = E/A;
    valueType q0, q1, q2, q3, scale;
    indexType i_case;
    scaleQuarticMonicPolynomial( A3, A2, A1, A0, q3, q2, q1, q0, i_case, scale );

    // Q(x)    = q0 + q1 * x + q2 * x^2 + q3 * x^3 + x^4
    // Q'(x)/4 = q1/4 + (q2/2) * x + (3/4)*q3 * x^2 + x^3

    /*
    ..  The general quartic case. Search for stationary points. Set the first
    ..  derivative polynomial (cubic) equal to zero and find its roots.
    ..  Check, if any minimum point of Q(x) is below zero, in which case we
    ..  must have real roots for Q(x). Hunt down only the real root, which
    ..  will potentially converge fastest during Newton iterates. The remaining
    ..  roots will be determined by deflation Q(x) -> cubic.
    ..
    ..  The best roots for the Newton iterations are the two on the opposite
    ..  ends, i.e. those closest to the +2 and -2. Which of these two roots
    ..  to take, depends on the location of the Q(x) minima x = s and x = u,
    ..  with s > u. There are three cases:
    ..
    ..     1) both Q(s) and Q(u) < 0
    ..        ----------------------
    ..
    ..        The best root is the one that corresponds to the lowest of
    ..        these minima. If Q(s) is lowest -> start Newton from +2
    ..        downwards (or zero, if s < 0 and a0 > 0). If Q(u) is lowest
    ..        -> start Newton from -2 upwards (or zero, if u > 0 and a0 > 0).
    ..
    ..     2) only Q(s) < 0
    ..        -------------
    ..
    ..        With both sides +2 and -2 possible as a Newton starting point,
    ..        we have to avoid the area in the Q(x) graph, where inflection
    ..        points are present. Solving Q''(x) = 0, leads to solutions
    ..        x = -a3/4 +/- discriminant, i.e. they are centered around -a3/4.
    ..        Since both inflection points must be either on the r.h.s or l.h.s.
    ..        from x = s, a simple test where s is in relation to -a3/4 allows
    ..        us to avoid the inflection point area.
    ..
    ..     3) only Q(u) < 0
    ..        -------------
    ..
    ..        Same of what has been said under 2) but with x = u.
    */

    valueType c2 = 3 * (q3/4);
    valueType c1 = q2/2;
    valueType c0 = q1/4;

    Cubic qsolve( 1, c2, c1, c0 );
    valueType u = qsolve.real_root0(); // root according to paper
    valueType t = qsolve.real_root1();
    valueType s = qsolve.real_root2();

    bool must_refine_r3 = true;
    if ( !qsolve.complexRoots() ) {
      valueType Qs = evalMonicQuartic( s, q3, q2, q1, q0 );
      valueType Qu = evalMonicQuartic( u, q3, q2, q1, q0 );
      bool q0pos = q0 > 0;
      nreal = 1;
      if ( Qs < 0 && Qu < 0 ) {
        if ( Qs < Qu ) {
          r3 = 2;
          if ( q0pos && s < 0 ) r3 = 0;
        } else {
          r3 = -2;
          if ( q0pos && u > 0 ) r3 = 0;
        }
      } else if ( Qs < 0 ) {
        if ( 4*s < -q3 ) {
          r3 = -2;
          if ( q0pos && s > 0 ) r3 = 0;
        } else {
          r3 = 2;
          if ( q0pos && s < 0 ) r3 = 0;
        }
      } else if ( Qu < 0 ) {
        if ( 4*u < -q3 ) {
          r3 = -2;
          if ( q0pos && u > 0 ) r3 = 0;
        } else {
          r3 = 2;
          if ( q0pos && u < 0 ) r3 = 0;
        }
      } else {
        // check for astrological combination when s or u are root of the quartic
        /*
        DO NOT WORK
        valueType epsi = 10 * ( 1 + std::abs(q3) +
                                    std::abs(q2) +
                                    std::abs(q1) +
                                    std::abs(q0) ) * machepsi;
        if      ( Qs <= epsi ) r3 = s;
        else if ( Qu <= epsi ) r3 = u;
        else                   nreal = 0;
        */
        if ( isZero(Qs) ) {
          must_refine_r3 = false; r3 = s;
        } else if ( isZero(Qu) ) {
          must_refine_r3 = false; r3 = u;
        } else {
          nreal = 0;
        }
      }
    } else {
      // one single real root (only 1 minimum)
      valueType Qs = evalMonicQuartic( s, q3, q2, q1, q0 );
      if ( Qs <= 0 ) {
        valueType tmp = q0 >= 0 ? 0 : 2;
        r3 = u > 0 ? -tmp : -2;
        nreal = 1;
      }
    }

    /*
    ..  Do all necessary Newton iterations. In case we have more than 2 oscillations,
    ..  exit the Newton iterations and switch to bisection. Note, that from the
    ..  definition of the Newton starting point, we always have Q(x) > 0 and Q'(x)
    ..  starts (-ve/+ve) for the (-2/+2) starting points and (increase/decrease) smoothly
    ..  and staying (< 0 / > 0). In practice, for extremely shallow Q(x) curves near the
    ..  root, the Newton procedure can overshoot slightly due to rounding errors when
    ..  approaching the root. The result are tiny oscillations around the root. If such
    ..  a situation happens, the Newton iterations are abandoned after 3 oscillations
    ..  and further location of the root is done using bisection starting with the
    ..  oscillation brackets.
    */
    if ( nreal > 0 ) {
      if ( must_refine_r3 )
        iter += zeroQuarticByNewtonBisection( q3, q2, q1, q0, r3 );
      r3 *= scale;

      /*
      ..  Find remaining roots -> reduce to cubic. The reduction to a cubic polynomial
      ..  is done using composite deflation to minimize rounding errors. Also, while
      ..  the composite deflation analysis is done on the reduced quartic, the actual
      ..  deflation is being performed on the original quartic again to avoid enhanced
      ..  propagation of root errors.
      */
      deflateQuarticPolynomial( A, B, C, D, E, r3, q2, q1, q0 );

      Cubic csolve( A, q2, q1, q0 );
      r0 = csolve.real_root0();
      r1 = csolve.real_root1();
      r2 = csolve.real_root2();
      if ( csolve.complexRoots() ) {
        nreal = ncplx = 2;
        if ( r2 > r3 ) std::swap( r2, r3 );
      } else {
        nreal = 4;
        if ( r2 > r3 ) std::swap( r2, r3 );
        if ( r1 > r2 ) std::swap( r1, r2 );
        if ( r0 > r1 ) std::swap( r0, r1 );
      }
    } else {
      /*
      .. If no real roots have been found by now, only complex roots are
      .. possible. Find real parts of roots first, followed by imaginary
      .. components.
      */
      s = q3/2;
      t = s * s - q2;
      u = s * t + q1; // value of Q'(-q3/4) at stationary point -q3/4
      bool notZero = (abs(u) >= machepsi); // H(-a3/4) is considered > 0 at stationary point
      bool minimum;
      if ( !isZero(q3) ) {
        s = q1 / q3;
        minimum = q0 > s * s; // H''(-q3/4) > 0 -> minimum
      } else {
        minimum = 4 * q0 > q2 * q2; // H''(-q3/4) > 0 -> minimum
      }

      bool iterate = notZero || (!notZero && minimum);

      valueType a, b, c, d;
      if ( iterate ) {
        valueType x = q3 >= 0 ? 2 : -2; // initial root -> target = smaller mag root
        iter += zeroHexicByNewtonBisection( q3, q2, q1, q0, x );

        a = x*scale;   // 1st real component -> a
        b = -A3/2 - a; // 2nd real component -> b

        x = 4 * a + A3; // Q'''(a)
        valueType y = x + 2*A3;
        y = y * a + 2*A2; // Q'(a)
        y = y * a + A1;
        y /= x;
        if ( y < 0 ) y = 0; // ensure Q'(a) / Q'''(a) >= 0

        x = 4 * b + A3; // Q'''(b)
        valueType z = x + 2*A3;
        z = z * b + 2*A2; // Q'(b)
        z = z * b + A1;
        z /= x;
        if ( z < 0 ) z = 0; // ensure Q'(b) / Q'''(b) >= 0

        c = a * a;              // store a^2 for later
        d = b * b;              // store b^2 for later
        valueType ss = c + y;   // magnitude^2 of (a + iy) root
        valueType tt = d + z;   // magnitude^2 of (b + iz) root

        if ( ss > tt ) {         // minimize imaginary error
          c = std::sqrt(y);           // 1st imaginary component -> c
          d = std::sqrt(A0 / ss - d); // 2nd imaginary component -> d
        } else {
          c = std::sqrt(A0 / tt - c); // 1st imaginary component -> c
          d = std::sqrt(z);           // 2nd imaginary component -> d
        }

      } else { // no bisection -> real components equal

        a = -A3/4; // 1st real component -> a
        b = a;     // 2nd real component -> b = a

        valueType x = a + A3;
        x = x * a + A2; // Q(a)
        x = x * a + A1;
        x = x * a + A0;
        valueType y = A2/2 - 3*(a*a); // Q''(a) / 2
        valueType z = y * y - x;
        z = z > 0 ? std::sqrt(z) : 0;     // force discriminant to be >= 0
                                     // square root of discriminant
        y = y > 0 ? y + z : y - z;   // larger magnitude root
        x /= y;                      // smaller magnitude root
        c = y < 0 ? 0 : std::sqrt(y);     // ensure root of biquadratic > 0
        d = x < 0 ? 0 : std::sqrt(x);     // large magnitude imaginary component
      }

      ncplx = 4;
      if      (a > b) { r0 = a; r1 = c; r2 = b; r3 = d; }
      else if (a < b) { r0 = b; r1 = d; r2 = a; r3 = c; }
      else            { r0 = a; r1 = c; r2 = a; r3 = d; }
    }
  }

  void
  Quartic::info( std::ostream & s ) const {
    valueType const & A = ABCDE[0];
    valueType const & B = ABCDE[1];
    valueType const & C = ABCDE[2];
    valueType const & D = ABCDE[3];
    valueType const & E = ABCDE[4];

    s << "\npoly a=" << A << " b=" << B << " c=" << C << " d=" << D << " e=" << E
      << "\nn. complex = " << ncplx
      << "\nn. real    = " << nreal;
    if ( ncplx > 0 ) {
      s << "\nx0 = (" << r0 << "," <<  r1 << ')'
        << "\nx1 = (" << r0 << "," << -r1 << ')';
    } else {
      if ( nreal > 0 ) s << "\nx0 = " << r0;
      if ( nreal > 1 ) s << "\nx1 = " << r1;
    }
    if ( ncplx > 2 ) {
      s << "\nx2 = (" << r2 << "," <<  r3 << ')'
        << "\nx3 = (" << r2 << "," << -r3 << ')';
    } else {
      if ( nreal > 2 || (ncplx > 0 && nreal > 0) ) s << "\nx2 = " << r2;
      if ( nreal > 3 || (ncplx > 0 && nreal > 1) ) s << "\nx3 = " << r3;
    }
    s << '\n';
  }

  bool
  Quartic::check( std::ostream & s ) const {
    valueType const & A = ABCDE[0];
    valueType const & B = ABCDE[1];
    valueType const & C = ABCDE[2];
    valueType const & D = ABCDE[3];
    valueType const & E = ABCDE[4];
    bool ok = true;
    valueType epsi = 1000 * ( ( std::abs(A) +
                                std::abs(B) +
                                std::abs(C) +
                                std::abs(D) +
                                std::abs(E) )*machepsi) ;
    if ( ncplx > 0 ) {
      valueType z0 = std::abs(eval( root0() ));
      valueType z1 = std::abs(eval( root1() ));
      s << "|p(r0)| = " << z0 << "\n|p(r1)| = " << z1 << '\n';
      ok = ok && std::abs(z0) < epsi && std::abs(z1) < epsi;
    } else {
      if ( nreal > 0 ) {
        valueType z0 = eval( real_root0() );
        s << "p(r0) = " << z0 << '\n';
        ok = ok && std::abs(z0) < epsi;
      }
      if ( nreal > 1 ) {
        valueType z1 = eval( real_root1() );
        s << "p(r1) = " << z1 << '\n';
        ok = ok && std::abs(z1) < epsi;
      }
    }
    if ( ncplx > 2 ) {
      valueType z2 = std::abs(eval( root2() ));
      valueType z3 = std::abs(eval( root3() ));
      s << "|p(r2)| = " << z2 << "\n|p(r3)| = " << z3 << '\n';
      ok = ok && std::abs(z2) < epsi && std::abs(z3) < epsi;
    } else {
      if ( nreal > 2 || (ncplx > 0 && nreal > 0)  ) {
        valueType z2 = eval( real_root2() );
        s << "p(r2) = " << z2 << '\n';
        ok = ok && std::abs(z2) < epsi;
      }
      if ( nreal > 3 || (ncplx > 0 && nreal > 1)  ) {
        valueType z3 = eval( real_root3() );
        s << "p(r3) = " << z3 << '\n';
        ok = ok && std::abs(z3) < epsi;
      }
    }
    return ok;
  }

}

#endif

// EOF: PolynomialRoots-3-Quartic.cc
