// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-FileCopyrightText: Copyright 2007 Sandia Corporation
// SPDX-License-Identifier: LicenseRef-BSD-3-Clause-Sandia-USGov
#include "vtkPolynomialSolversUnivariate.h"
#include "vtkDataArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"

#include <cmath>

#define VTK_SIGN(x) (((x) < 0) ? (-1) : (1))

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkPolynomialSolversUnivariate);

static const double sqrt3 = sqrt(3.);
static constexpr double inv3 = 1 / 3.;
static const double absolute0 = 10. * VTK_DBL_MIN;

double vtkPolynomialSolversUnivariate::DivisionTolerance = 1e-8; // sqrt( VTK_DBL_EPSILON );

//------------------------------------------------------------------------------
void vtkPolynomialSolversUnivariate::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent
     << "(s) DivisionTolerance: " << vtkPolynomialSolversUnivariate::GetDivisionTolerance() << "\n";
}

//------------------------------------------------------------------------------
// Print the polynomial for debuggery
ostream& vtkPolynomialSolversUnivariate::PrintPolynomial(ostream& os, double* P, int degP)
{
  os << "\n";
  os << "The polynomial has degree " << degP << "\n";

  if (degP < 0)
  {
    os << "0\n";
    return os;
  }

  if (degP == 0)
  {
    os << P[0] << "\n";
    return os;
  }

  unsigned int degPm1 = degP - 1;
  for (unsigned int i = 0; i < degPm1; ++i)
  {
    if (P[i] > 0)
    {
      if (i)
      {
        os << "+";
      }
      if (P[i] != 1.)
      {
        os << P[i] << "*";
      }
      os << "x**" << degP - i;
    }
    else if (P[i] < 0)
    {
      os << P[i] << "*x**" << degP - i;
    }
  }

  if (degP > 0)
  {
    if (P[degPm1] > 0)
    {
      os << "+" << P[degPm1] << "*x";
    }
    else if (P[degPm1] < 0)
    {
      os << P[degPm1] << "*x";
    }
  }

  if (P[degP] > 0)
  {
    os << "+" << P[degP];
  }
  else if (P[degP] < 0)
  {
    os << P[degP];
  }

  os << "\n";
  return os;
}

//------------------------------------------------------------------------------
// Double precision comparison with 0
inline bool IsZero(double x)
{
  return fabs(x) < absolute0;
}

//------------------------------------------------------------------------------
// Double precision comparison
inline bool AreEqual(double x, double y, double rTol)
{
  double delta = fabs(x - y);

  // First, handle "absolute" zeros. This is to eliminate the
  // case (x+t - -x ) / x = 2 + t/x even when both x and t are small.
  if (delta < absolute0)
  {
    return true;
  }

  // Second, handle "relative" equalities.
  double absx = fabs(x);
  double absy = fabs(y);
  if (absx > absy)
  {
    return delta <= rTol * absx;
  }
  else
  {
    return delta <= rTol * absy;
  }
}

//------------------------------------------------------------------------------
// Polynomial Euclidean division of A (deg m) by B (deg n).
static int polynomialEucliDiv(double* A, int m, double* B, int n, double* Q, double* R, double rtol)
{
  // Note: for execution speed, no sanity checks are performed on A and B.
  // You must know what you are doing.

  int mMn = m - n;
  int i;

  if (mMn < 0)
  {
    Q[0] = 0.;
    for (i = 0; i <= m; ++i)
    {
      R[i] = A[i];
    }

    return m;
  }

  double iB0 = 1. / B[0];
  if (!n)
  {
    for (i = 0; i <= m; ++i)
    {
      Q[i] = A[i] * iB0;
    }

    return -1;
  }

  int nj;
  for (i = 0; i <= mMn; ++i)
  {
    nj = i > n ? n : i;
    Q[i] = A[i];
    for (int j = 1; j <= nj; ++j)
    {
      Q[i] -= B[j] * Q[i - j];
    }
    Q[i] *= iB0;
  }

  bool nullCoeff = false;
  int r = 0;
  for (i = 1; i <= n; ++i)
  {
    double sum = 0;
    nj = mMn + 1 > i ? i : mMn + 1;
    for (int j = 0; j < nj; ++j)
    {
      sum += B[n - i + 1 + j] * Q[mMn - j];
    }

    if (!AreEqual(A[m - i + 1], sum, rtol))
    {
      R[n - i] = A[m - i + 1] - sum;
      r = i - 1;
    }
    else
    {
      R[n - i] = 0.;
      if (n == i)
      {
        nullCoeff = true;
      }
    }
  }

  if (!r && nullCoeff)
  {
    return -1;
  }

  return r;
}

//------------------------------------------------------------------------------
// Polynomial Euclidean division of A (deg m) by B (deg n).
// Does not store Q and stores -R instead of R
static int polynomialEucliDivOppositeR(double* A, int m, double* B, int n, double* mR, double rtol)
{
  // Note: for execution speed, no sanity checks are performed on A and B.
  // You must know what you are doing.

  int mMn = m - n;
  int i;

  if (mMn < 0)
  {
    for (i = 0; i <= m; ++i)
    {
      mR[i] = A[i];
    }
    return m;
  }

  if (!n)
  {
    return -1;
  }

  int nj;
  double iB0 = 1. / B[0];
  double* Q = new double[mMn + 1];
  for (i = 0; i <= mMn; ++i)
  {
    nj = i > n ? n : i;
    Q[i] = A[i];
    for (int j = 1; j <= nj; ++j)
    {
      Q[i] -= B[j] * Q[i - j];
    }
    Q[i] *= iB0;
  }

  bool nullCoeff = false;
  int r = 0;
  for (i = 1; i <= n; ++i)
  {
    double sum = 0;
    nj = mMn + 1 > i ? i : mMn + 1;
    for (int j = 0; j < nj; ++j)
    {
      sum += B[n - i + 1 + j] * Q[mMn - j];
    }

    if (!AreEqual(A[m - i + 1], sum, rtol))
    {
      mR[n - i] = sum - A[m - i + 1];
      r = i - 1;
    }
    else
    {
      mR[n - i] = 0.;
      if (n == i)
      {
        nullCoeff = true;
      }
    }
  }
  delete[] Q;

  if (!r && nullCoeff)
  {
    r = -1;
  }

  return r;
}

//------------------------------------------------------------------------------
inline double vtkNormalizePolyCoeff(double d, double* div = nullptr)
{
  static const double high = 18446744073709551616.; // 2^64
  static const double reallyBig = 1.e300;
  static const double reallyBigInv = 1 / reallyBig;
  static const double notThatBig = 1.e30;
  static const double notThatBigInv = 1.e-30;

  if (fabs(d) < reallyBig)
  {
    while (fabs(d) > notThatBig)
    {
      d /= high;
      if (div)
      {
        *div /= high;
      }
    }
  }
  if (fabs(d) > reallyBigInv)
  {
    while (fabs(d) < notThatBigInv)
    {
      d *= high;
      if (div)
      {
        *div *= high;
      }
    }
  }
  return d;
}

//------------------------------------------------------------------------------
// Polynomial Euclidean division of A (deg m) by B (deg n).
// Does not store Q and stores -R instead of R. This premultiplies Ai by mul and
// then divides mR by div. mR MUST have at least the size of m+1, because it is
// used as temporary storage besides as the return value.
// This shouldn't be too much slower than the version without mul and div. It
// has to copy mul*Ai into mR at the beginning. So it has m+1 extra
// multiplications and assignments. When it assigns final values into mR it
// divides by div. So it has deg(mR) extra divisions.
static int polynomialEucliDivOppositeR(
  double mul, double* Ai, int m, double* B, int n, double div, double* mR, double rtol)
{
  // Note: for execution speed, no sanity checks are performed on A and B.
  // You must know what you are doing.

  int mMn = m - n;
  int i;

  // To save space we will use mR *instead* of Ai. Just have to be sure that
  // when we write to mR we are writing to a spot that will not be used again.
  for (i = 0; i <= m; ++i)
  {
    mR[i] = mul * Ai[i];
  }

  // Remember that mMn >= 0 implies m-n >= 0 implies m  >= n and so
  // m + 1 > n.
  if (mMn < 0)
  {
    return m;
  }

  if (!n)
  {
    return -1;
  }

  div = 1 / div;
  int nj;
  double iB0 = 1. / B[0];
  double* Q = new double[mMn + 1];

  for (i = 0; i <= mMn; ++i)
  {
    nj = i > n ? n : i;
    Q[i] = mR[i];
    for (int j = 1; j <= nj; ++j)
    {
      Q[i] -= B[j] * Q[i - j];
    }
    Q[i] *= iB0;
  }

  bool nullCoeff = false;
  int r = 0;
  for (i = n; i >= 1; --i)
  {
    double sum = 0;
    nj = mMn + 1 > i ? i : mMn + 1;
    for (int j = 0; j < nj; ++j)
    {
      sum += B[n - i + 1 + j] * Q[mMn - j];
    }

    if (!AreEqual(mR[m - i + 1], sum, rtol))
    {
      // Now we have m + 1 > n implies n - i < m - i + 1. Thus since i is
      // decreasing (hence m - i + 1 increasing) we will never use n - i again,
      // so we can safely write over it.
      mR[n - i] = (sum - mR[m - i + 1]) * div;
      // We want the non-zero term with the largest index. So we only write to r
      // once.
      if (r == 0)
      {
        mR[n - i] = vtkNormalizePolyCoeff(mR[n - i], &div);
        r = i - 1;
      }
    }
    else
    {
      // See previous comment.
      mR[n - i] = 0.;
      if (n == i)
      {
        nullCoeff = true;
      }
    }
  }
  delete[] Q;

  if (!r && nullCoeff)
  {
    r = -1;
  }

  return r;
}

//------------------------------------------------------------------------------
// Evaluate the value of the degree d univariate polynomial P at x
// using Horner's algorithm.
inline double evaluateHorner(double* P, int d, double x)
{
  if (d == -1)
  {
    return 0.;
  }
  double val = P[0];
  for (int i = 1; i <= d; ++i)
  {
    val = val * x + P[i];
  }

  return val;
}

static int vtkGetSignChanges(
  double* P, int* degP, int* offsets, int count, double val, int* fsign = nullptr)
{
  int oldVal = 0;
  double v;
  int changes = 0;

  for (int i = 0; i < count; ++i)
  {
    v = evaluateHorner(P + offsets[i], degP[i], val);

    if (fsign && i == 0)
    {
      if (IsZero(v))
      {
        *fsign = 0;
      }
      else if (v > 0)
      {
        *fsign = 1;
      }
      else
      {
        *fsign = -1;
      }
    }

    if (v == 0)
    {
      continue;
    }

    if (v * oldVal < 0)
    {
      ++changes;
      oldVal = -oldVal;
    }

    if (oldVal == 0)
    {
      oldVal = (v < 0 ? -1 : 1);
    }
  }

  return changes;
}

// ----------------------------------------------------------
// Gets the Habicht sequence. SSS and degrees and offsets are expected to be
// large enough and the number of non-zero items is returned. P is expected to
// have degree at least 1.
//
// This algorithm is modified from
// BPR, Algorithms in Real Algebraic Geometry, page 318
static int vtkGetHabichtSequence(
  double* P, int d, double* SSS, int* degrees, int* offsets, double rtol)
{
  // const static double high = pow(2.,256);
  degrees[0] = d;
  offsets[0] = 0;

  int dp1 = d + 1;
  double* t = new double[dp1];
  double* s = new double[dp1];

  degrees[1] = d - 1;
  offsets[1] = dp1;
  int offset = dp1;

  // Set the first two elements SSS = {P, P'}.
  for (int m = 0; m < d; ++m)
  {
    SSS[m] = P[m];
    SSS[m + offset] = static_cast<double>(d - m) * SSS[m];
  }
  SSS[d] = P[d];
  t[0] = s[0] = (P[0] > 0. ? 1. : -1.);
  t[1] = s[1] = SSS[offset];

  int j = 0;

  int deg = d;
  int degree = d - 1;
  int jp1 = 1;
  int ip1 = 0;
  while (degree > 0 && j < d - 1)
  {
    int k = deg - degree;
    if (k == jp1)
    {
      s[jp1] = t[jp1];

      degrees[k + 1] = polynomialEucliDivOppositeR(s[jp1] * s[jp1], SSS + offsets[ip1],
        degrees[ip1], SSS + offset, degree, s[j] * t[ip1], SSS + offsets[k] + degree + 1, rtol);
      offsets[k + 1] = offset + 2 * degree - degrees[k + 1];
    }
    else
    {
      s[jp1] = 0;
      for (int delta = 1; delta < k - j; ++delta)
      {
        t[jp1 + delta] = (t[jp1] * t[j + delta]) / s[j];
        t[jp1 + delta] = vtkNormalizePolyCoeff(t[jp1 + delta]);
        if (delta % 2)
        {
          t[jp1 + delta] *= -1;
        }
      }
      s[k] = t[k];

      offsets[k] = offsets[jp1] + degrees[jp1] + 1;
      degrees[k] = degrees[jp1];
      for (int dg = 0; dg <= degree; ++dg)
      {
        SSS[offsets[k] + dg] = (s[k] * SSS[offset + dg]) / t[jp1];
      }

      for (int l = j + 2; l < k; ++l)
      {
        degrees[l] = -1;
        offsets[l] = offsets[k];
        s[l] = 0;
      }

      degrees[k + 1] = polynomialEucliDivOppositeR(t[jp1] * s[k], SSS + offsets[ip1], degrees[ip1],
        SSS + offset, degree, s[j] * t[ip1], SSS + offsets[k] + degrees[k] + 1, rtol);
      offsets[k + 1] = offsets[k] + 2 * degrees[k] - degrees[k + 1];
    }
    t[k + 1] = SSS[offsets[k + 1]];
    ip1 = jp1;
    j = k;
    jp1 = j + 1;
    degree = degrees[jp1];
    offset = offsets[jp1];
  }

  delete[] s;
  delete[] t;

  if (degree == 0)
  {
    return jp1 + 1;
  }
  else
  {
    while (degrees[jp1] < 0)
    {
      --jp1;
    }
    return jp1 + 1;
  }
}

// ----------------------------------------------------------
// Gets the sturm sequence. SSS and degrees and offsets are expected to
// be large enough and the number of non-zero items
// is returned. P is expected to have degree at least 1.
static int vtkGetSturmSequence(
  double* P, int d, double* SSS, int* degrees, int* offsets, double rtol)
{
  degrees[0] = d;
  offsets[0] = 0;

  int dp1 = d + 1;
  int dm1 = d - 1;
  degrees[1] = dm1;
  offsets[1] = dp1;
  int offset = dp1;

  // nSSS will keep track of the number of the index of
  // the last item in our list.
  int nSSS = 1;

  // Set the first two elements SSS = {P, P'}.
  for (int k = 0; k < d; ++k)
  {
    SSS[k] = P[k];
    SSS[k + offset] = static_cast<double>(d - k) * P[k];
  }
  SSS[d] = P[d];

  int degree = dm1;
  while (degrees[nSSS] > 0)
  {
    nSSS++;
    degrees[nSSS] = polynomialEucliDivOppositeR(SSS + offsets[nSSS - 2], degrees[nSSS - 2],
      SSS + offset, degree, SSS + offset + degree + 1, rtol);

    offsets[nSSS] = offset + 2 * degree - degrees[nSSS];

    offset = offsets[nSSS];
    degree = degrees[nSSS];
  }

  // If the last element is zero then we ignore it.
  if (degrees[nSSS] < 0)
  {
    return nSSS;
  }

  // Otherwise we include it in our count, because it
  // is a constant. (We are returning the number of
  // nonzero items, so 1 plus the index of the last).
  return nSSS + 1;
}

int vtkPolynomialSolversUnivariateCompareRoots(const void* a, const void* b)
{
  double aa = *static_cast<const double*>(a);
  double bb = *static_cast<const double*>(b);
  if (aa < bb)
  {
    return -1;
  }
  if (aa > bb)
  {
    return 1;
  }
  return 0;
}

// ------------------------------------------------------------
// upperBnds is expected to be large enough.
// intervalType finds roots as follows
// 0 = ]a,b[
// 1 = [a,b[
// 2 = ]a,a]
// 3 = [a,b]
//
// divideGCD uses an integer in case in the future someone wants to add logic to
// whether the gcd is divided. For example something like
// divideGCD == 0 -> never divide
// divideGCD == 1 -> use logic
// divideGCD == 2 -> divide as long as non-constant (ignore the logic).
//
// It probably would have been better to originally have had the tolerance be a
// relative tolerance rather than an absolute tolerance.
static int vtkHabichtOrSturmBisectionSolve(double* P, int d, double* a, double* upperBnds,
  double tol, int intervalType, int divideGCD, int method)
{
  // Pretend to be one solver or the other (for error reporting)
  static const char title1[] = "vtkPolynomialSolversUnivariate::SturmBisectionSolve";
  static const char title2[] = "vtkPolynomialSolversUnivariate::HabichtBisectionSolve";
  const char* title = (method == 0 ? title1 : title2);

  // 0. Stupidity checks
  if (tol <= 0)
  {
    vtkGenericWarningMacro(<< title << ": Tolerance must be positive");
    return -1;
  }

  if (IsZero(P[0]))
  {
    vtkGenericWarningMacro(<< title << ": Zero leading coefficient");
    return -1;
  }

  if (d < 1)
  {
    vtkGenericWarningMacro(<< title << ": Degree (" << d << ") < 1");
    return -1;
  }

  if (a[1] < a[0] + tol)
  {
    vtkGenericWarningMacro(<< title << ": Erroneous interval endpoints and/or tolerance");
    return -1;
  }

  // Check for 0 as a root and reduce the degree if so.
  bool zeroroot = false;
  if (IsZero(P[d]))
  {
    zeroroot = true;
    while (IsZero(P[d]))
    {
      --d;
    }
  }

  // Take care of constant polynomials and polynomials of the form a*x^d.
  if (d == 0)
  {
    if (zeroroot)
    {
      upperBnds[0] = 0.;
      return 1;
    }
    else
    {
      return 0;
    }
  }

  // Create one large array to hold all the
  // polynomials.
  //
  // We need two extra spaces because the habicht division uses R to temporarily
  // hold mul*A. So there must be at least size(A) space when dividing, which can be
  // 2 floats larger than R.
  double* SSS = new double[((d + 1) * (d + 2)) / 2 + 2];
  int* degrees = new int[d + 1];
  int* offsets = new int[d + 1];

  double bounds[] = { a[0], a[1] };
  int nSSS;

  if (method == 0)
  {
    nSSS = vtkGetSturmSequence(
      P, d, SSS, degrees, offsets, vtkPolynomialSolversUnivariate::GetDivisionTolerance());
  }
  else
  {
    nSSS = vtkGetHabichtSequence(
      P, d, SSS, degrees, offsets, vtkPolynomialSolversUnivariate::GetDivisionTolerance());
  }

  // If degrees[count-1] > 0 then we have degenerate roots.
  // We could possibly then find the degenerate roots.
  // Maybe we could more or less remove the degenerate roots.
  // Lets do some testwork to see.
  if (degrees[nSSS - 1] > 0 && divideGCD == 1)
  {
    double* R = new double[d + 1];
    double* Q = new double[d + 1];

    // Get the quotient and call this function again using the
    // quotient.
    int deg = polynomialEucliDiv(SSS, d, SSS + offsets[nSSS - 1], degrees[nSSS - 1], Q, R,
      vtkPolynomialSolversUnivariate::GetDivisionTolerance());
    (void)deg; // above result not really needed, but useful in debugging
    deg = d - degrees[nSSS - 1];

    delete[] R;
    // The Habicht sequence will occasionally get infinite coeffs and cause
    // unpleasant events to happen with the sequence. In that case Q[0] == 0, thus
    // the division is not used.
    if (!IsZero(Q[0]))
    {
      delete[] SSS;
      delete[] degrees;
      delete[] offsets;

      int rval =
        vtkHabichtOrSturmBisectionSolve(Q, deg, a, upperBnds, tol, intervalType, 0, method);
      delete[] Q;
      if (zeroroot)
      {
        upperBnds[rval] = 0;
        return rval + 1;
      }
      return rval;
    }
    else
    {
      delete[] Q;
    }
  }

  // Move away from zeros on the edges. We can also slightly speed up
  // computation by keeping the fact that these are roots and
  // continuing on.

  // If perturbation is too small compared to the bounds then we won't move when
  // we perturb.
  double perturbation =
    fmax(fmax(fabs(bounds[0]) * 1e-12, fabs(bounds[1]) * 1e-12), .5 * tol / static_cast<double>(d));

  int varSgn[] = { 0, 0 };
  varSgn[0] = vtkGetSignChanges(SSS, degrees, offsets, nSSS, bounds[0]);
  varSgn[1] = vtkGetSignChanges(SSS, degrees, offsets, nSSS, bounds[1]);

  for (int k = 0; k <= 1; ++k)
  {
    if (IsZero(evaluateHorner(SSS, d, bounds[k])))
    {
      int leftVarSgn = varSgn[k];
      int rightVarSgn = varSgn[k];
      double leftx = bounds[k];
      double rightx = bounds[k];
      // Make sure we move far enough away that everything still works. That is
      // we needs to be non-zero and have the sequence realize that we've
      // got a zero in the interval. It is possible that this causes our code
      // to slow down. Should probably play around with how large the
      // perturbation should be. Also it is possible (but unlikely) that our inexact
      // sturm sequence doesn't realize there is a zero here.
      //
      // JUST AS WITH THE BISECTING, NEED TO MAKE SURE WE DON'T HAVE AN
      // INFINITE LOOP.
      while (IsZero(evaluateHorner(SSS, d, leftx)) || IsZero(evaluateHorner(SSS, d, rightx)) ||
        leftVarSgn <= rightVarSgn ||
        ((leftVarSgn == varSgn[k] || rightVarSgn == varSgn[k]) && leftVarSgn - rightVarSgn != 1))
      {
        leftx -= perturbation;
        rightx += perturbation;
        leftVarSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, leftx);
        rightVarSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, rightx);
      }

      // Move properly according to what kind of sequence we are searching.
      if ((!(intervalType & 2) && k == 1) || ((intervalType & 1) && k == 0))
      {
        bounds[k] = leftx;
        varSgn[k] = leftVarSgn;
      }
      else
      {
        bounds[k] = rightx;
        varSgn[k] = rightVarSgn;
      }
    }
  }

  // If we don't have roots then leave here.
  int nRoots = varSgn[0] - varSgn[1];
  if (nRoots < 1)
  {
    upperBnds[0] = 0;
    delete[] SSS;
    delete[] degrees;
    delete[] offsets;
    if (zeroroot)
      return 1;
    return 0;
  }

  // 2. Root bracketing

  // Initialize the bounds for the root intervals. The interval
  // ]lowerBnds[i], upperBnds[i][ contains the i+1 root. We will
  // see if we can completely separate the roots. Of course
  // the interval ]bounds[0], bounds[1][ contains all the roots.
  // Afterwards if some intervals are the same all but
  // one will be removed.
  int i;
  double* lowerBnds = new double[nRoots];

  for (i = 0; i < nRoots; ++i)
  {
    upperBnds[i] = bounds[1];
    lowerBnds[i] = bounds[0];
  }

  int leftVarSgn, rightVarSgn, tempSgn;
  double leftx, rightx;
  bool hitroot;
  int nloc = nRoots - 1;
  while (nloc >= 1)
  {
    // Only one root according to Sturm or the interval is
    // small enough to consider the same root.
    if (upperBnds[nloc] - lowerBnds[nloc] <= tol ||
      ((nloc < 1 || (upperBnds[nloc - 1] < lowerBnds[nloc] - tol)) &&
        ((nloc >= nRoots - 1) || (upperBnds[nloc]) < lowerBnds[nloc + 1] - tol)))
    {
      --nloc;
      continue;
    }

    // We begin with leftx and rightx being equal and change them only if
    // leftx (rightx) is a root. Then we can bracket the root. We do this
    // because roots can cause problems (our sequence is inexact so
    // even single roots can cause problems. Furthermore if we hit a

    // root we may as well bracket it within tol so that we don't have to
    // worry about it later.
    leftx = (upperBnds[nloc] + lowerBnds[nloc]) / 2;
    // If we are going nowhere then quit.
    if (leftx >= upperBnds[nloc] || leftx <= lowerBnds[nloc])
    {
      nloc--;
      continue;
    }
    rightx = leftx;
    hitroot = false;

    leftVarSgn = rightVarSgn = tempSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, rightx);

    // Move away if we have a root, just like we did with the initial endpoints.
    // Maybe could write a function to do this, but it may be slower.
    // After moving away we want the following things to be handled:
    // 1 Not zero at the endpoints (leftx or rightx)
    // 2 Sign changes at left != sign changes at right.
    // 3 No "crazy" values:
    //   a. sign[left] > sign[0]
    //   b. sign[right] < sign[1]
    //   c. sign[right] > sign[left].
    // 4 Does not take too long.
    //
    // It would be convenient if we could get sign[left]-sign[right] = 1, but
    // in reality we may be too close for our values to even be worthwhile.
    //
    // The question remains, "What happens if we move away too far?" There must
    // be some distance after which we quit perturbing and deal with the fact
    // that our sequence is not perfect. This will take care of 4, but we still
    // need to handle what has just happened.
    //
    // Option 2: In this step we are only worried about bracketing the
    // roots. If the midpoint is a root, we could try this. Suppose a=0, b=1.
    // Then if we get 1/2 is a root, we check whether 1/4 and 3/4 work. If
    // one of them does then call that our "mid". Then in the next step we
    // will be moved away from 1/2. If neither of them work then we could
    // do what we are doing here, or try 1/8, 3/8, 5/8, and 7/8. In this latter
    // case we should make sure we don't come back to that messy point later.
    // This option removes the use of one of leftx or rightx.
    // So we need to have the following not happen.
    //
    // 1. Not zero at mid.
    // 2. No "crazy" values:
    //    a. sign[mid] > sign[0]
    //    b. sign[mid] < sign[1].
    // 3. Does no take too long.
    //
    // This method is better because the "perturbation" is always relative to
    // the endpoints. Unlike the other method, because the perturbation
    // variable is relative to the larger endpoing.
    if (IsZero(leftx) || IsZero(evaluateHorner(SSS, d, leftx)) || (tempSgn > varSgn[0]) ||
      (tempSgn < nloc))
    {
      int step = 2;
      int pos = 1;
      double p2 = 4.;
      double mid = upperBnds[nloc] / p2 + (p2 - pos) * lowerBnds[nloc] / p2;
      bool found = false;
      leftVarSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, lowerBnds[nloc]);
      rightVarSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, upperBnds[nloc]);
      tempSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, mid);
      while (
        // 3.
        step < 10 &&
        // 2a.
        ((tempSgn > leftVarSgn) ||
          // 2b.
          (tempSgn < rightVarSgn) ||
          // 1
          IsZero(evaluateHorner(SSS, d, mid)) || IsZero(mid)))
      {
        pos += 2;
        if (pos > p2)
        {
          pos = 1;
          ++step;
          p2 *= 2.;
        }
        mid = pos * upperBnds[nloc] / p2 + (p2 - pos) * lowerBnds[nloc] / p2;
        tempSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, mid);
      }

      if (step < 10)
      {
        found = true;
        leftx = rightx = mid;
        leftVarSgn = rightVarSgn = tempSgn;
        // Set the new information for the current.
        if (varSgn[0] - leftVarSgn <= nloc)
          lowerBnds[nloc] = leftx;
        if (varSgn[0] - rightVarSgn > nloc)
          upperBnds[nloc] = rightx;
      }

      hitroot = !found;
      // Make sure all our measures change around the root. This is another place where
      // the perturbation may be too small and cause slowness. This could
      // theoretically cause an infinite loop. (It has!).
      while (!found &&
        (
          // 1
          IsZero(evaluateHorner(SSS, d, leftx)) ||
          IsZero(evaluateHorner(SSS, d, rightx))
          // 2 and 3c.
          || leftVarSgn <= rightVarSgn
          // 3a and 3b
          || leftVarSgn > varSgn[0] || rightVarSgn < varSgn[1]))
      {
        leftx -= perturbation;
        rightx += perturbation;
        // Take care of 4.
        if (rightx - leftx > 2 * tol)
          break;
        leftVarSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, leftx);
        rightVarSgn = vtkGetSignChanges(SSS, degrees, offsets, nSSS, rightx);
      }
      // Now we must take care of our possible blunders.
      if (rightx - leftx > 2 * tol)
      {
        // If leftx <= upperBnds[nloc-1] then we have to gracefully
        // clean up the mess.
        // For now assume we moved well enough...
        leftVarSgn = std::min(leftVarSgn, varSgn[0]);

        rightVarSgn = std::max(rightVarSgn, varSgn[1]);

        if (rightVarSgn > varSgn[0])
          rightVarSgn = varSgn[0] - nloc + 1;

        if (leftVarSgn < varSgn[1])
          leftVarSgn = varSgn[0] - nloc;

        rightx += tol;
        leftx -= tol;
      }

      if (hitroot)
      {
        lowerBnds[nloc] = mid;
        upperBnds[nloc] = mid;
      }
    }
    else
    {
      // Set the new information for the current.
      if (varSgn[0] - leftVarSgn <= nloc)
        lowerBnds[nloc] = leftx;
      if (varSgn[0] - rightVarSgn > nloc)
        upperBnds[nloc] = rightx;
    }

    // We have isolated the rightVarSgn to the leftVarSgn roots.
    if (rightx != leftx)
    {
      for (i = varSgn[0] - leftVarSgn; i <= varSgn[0] - rightVarSgn - 1; ++i)
      {
        if (i > 0 && lowerBnds[i - 1] < leftx)
        {
          lowerBnds[i] = leftx;
        }
        upperBnds[i] = std::min(upperBnds[i], rightx);
      }
    }

    // Set the new lower bounds for the intervals to the right of the new and the left of the
    // current. We have to do the check because we can get crazy information with some points.
    for (i = varSgn[0] - rightVarSgn; i >= 0 && i < nRoots; ++i)
      if (lowerBnds[i] < rightx && upperBnds[i] > rightx)
        lowerBnds[i] = rightx;

    // We have to do the check because we can get crazy information with some points.
    // Set the new upper bounds for the intervals to the left of the new interval.
    for (i = 0; i < varSgn[0] - leftVarSgn && varSgn[0] - leftVarSgn <= nloc; ++i)
      if (upperBnds[i] > leftx && lowerBnds[i] < leftx)
        upperBnds[i] = leftx;

    if (leftVarSgn - rightVarSgn == 1 || hitroot)
      --nloc;
  }

  int nIntervals = nRoots;
  bool* bisection = new bool[nRoots];

  // 3. Root polishing (if needed)
  for (nloc = 0; nloc < nRoots; ++nloc)
  {
    if (upperBnds[nloc] - lowerBnds[nloc] < tol)
    {
      continue;
    }

    double zv = evaluateHorner(P, d, upperBnds[nloc]);
    double lv = evaluateHorner(P, d, lowerBnds[nloc]);
    // These are for the possibility that we get sign change differences when
    // using the sequence.
    int ls, zs;
    double z;

    // Neither of the next two conditions should ever be met, but when using
    // floats we still want to be sure.
    if (IsZero(zv))
    {
      lowerBnds[nloc] = upperBnds[nloc];
      continue;
    }

    if (IsZero(lv))
    {
      upperBnds[nloc] = lowerBnds[nloc];
      continue;
    }

    int us = (zv > 0) ? 1 : -1;
    ls = (lv > 0) ? 1 : -1;

    bool bisect = false;
    // Check to see if we should use the sequence, which is MUCH slower.
    if (us * ls > 0)
    {
      int zc;
      while (upperBnds[nloc] - lowerBnds[nloc] > tol)
      {
        z = (upperBnds[nloc] + lowerBnds[nloc]) / 2;
        // Sometimes the tolerance can be poorly chosen causing an infinite
        // loop, this should fix that. That is, u - l > tol, but (u+l)/2 == l|u.
        if (z >= upperBnds[nloc] || z <= lowerBnds[nloc])
        {
          break;
        }
        zc = vtkGetSignChanges(SSS, degrees, offsets, nSSS, z, &zs);

        // If that point is a zero go ahead and quit.
        if (zs == 0)
        {
          upperBnds[nloc] = lowerBnds[nloc] = z;
          break;
        }

        // Decide which side it goes on.
        if ((varSgn[0] - zc) == nloc + 1)
        {
          us = zs;
          upperBnds[nloc] = z;
        }
        else
        {
          ls = zs;
          lowerBnds[nloc] = z;
        }

        // Maybe we are *now* close enough to use bisection. In which case DO THAT!
        // One extra multiplication is worth it.
        if (us * ls < 0)
        {
          bisect = true;
          // Quit this loop so that we can use bisection.
          break;
        }

      } // While ub[nloc] - lb[nloc] > tol

      bisection[nloc] = false;
      if (!bisect)
      {
        // Move on to next interval.
        continue;
      }
    } // If ub and lb have the same sign.
    else
    {
      bisect = true;
    }

    // If we can, use bisection.
    if (bisect)
    {
      double tempu = zv;
      while (upperBnds[nloc] - lowerBnds[nloc] > tol)
      {
        z = (upperBnds[nloc] + lowerBnds[nloc]) / 2;
        // Sometimes the tolerance can be poorly chosen causing an infinite
        // loop, this should fix that. That is, u - l > tol, but (u+l)/2 == l|u.
        if (z >= upperBnds[nloc] || z <= lowerBnds[nloc])
        {
          break;
        }
        zv = evaluateHorner(P, d, z);
        if (IsZero(zv))
        {
          upperBnds[nloc] = lowerBnds[nloc] = z;
          break;
        }

        if (zv * tempu > 0)
        {
          tempu = zv;
          upperBnds[nloc] = z;
        }
        else
          lowerBnds[nloc] = z;
      }
      bisection[nloc] = true;
    } // While ub[nloc] - lb[nloc] > tol
  }   // If we can bisect

  // Though theoretically this shouldn't happen, sometimes the roots are out of
  // order. Lets sort them just in case, because it does happen sometimes.
  qsort(upperBnds, nIntervals, sizeof(double), vtkPolynomialSolversUnivariateCompareRoots);
  qsort(lowerBnds, nIntervals, sizeof(double), vtkPolynomialSolversUnivariateCompareRoots);

  // Remove duplicate roots.
  for (int j = 1; j < nIntervals; ++j)
  {
    if (upperBnds[j] < (upperBnds[j - 1] + 2 * tol) ||
      lowerBnds[j] < (lowerBnds[j - 1] + 2 * tol) || (zeroroot && fabs(upperBnds[j]) < 2 * tol))
    {
      for (int k = j + 1; k < nIntervals; ++k)
      {
        upperBnds[k - 1] = upperBnds[k];
        lowerBnds[k - 1] = lowerBnds[k];
      }
      --j;
      --nIntervals;
      continue;
    }
  }

  delete[] bisection;
  delete[] lowerBnds;

  delete[] degrees;
  delete[] SSS;
  delete[] offsets;

  // Make sure the first root isn't zero.
  if (zeroroot && (fabs(upperBnds[0]) < 2 * tol))
  {
    for (int k = 1; k < nIntervals; ++k)
    {
      upperBnds[i - 1] = upperBnds[i];
    }
  }

  if (zeroroot)
  {
    upperBnds[nIntervals] = 0;
    ++nIntervals;
  }

  return nIntervals;
}

//------------------------------------------------------------------------------
// Find all real roots in ] a[0] ; a[1] [ of a real
// d-th degree polynomial using the Habicht sequence.
// intervalType specifies as follows (in binary)
// 0 = 00 = ]a,b[
// 1 = 10 = [a,b[
// 2 = 01 = ]a,b]
// 3 = 11 = [a,b]
int vtkPolynomialSolversUnivariate::HabichtBisectionSolve(
  double* P, int d, double* a, double* upperBnds, double tol)
{
  return vtkHabichtOrSturmBisectionSolve(P, d, a, upperBnds, tol, 0, 0, 1);
}

int vtkPolynomialSolversUnivariate::HabichtBisectionSolve(
  double* P, int d, double* a, double* upperBnds, double tol, int intervalType)
{
  return vtkHabichtOrSturmBisectionSolve(P, d, a, upperBnds, tol, intervalType, 0, 1);
}

int vtkPolynomialSolversUnivariate::HabichtBisectionSolve(
  double* P, int d, double* a, double* upperBnds, double tol, int intervalType, bool divideGCD)
{
  return vtkHabichtOrSturmBisectionSolve(
    P, d, a, upperBnds, tol, intervalType, divideGCD ? 1 : 0, 1);
}

//------------------------------------------------------------------------------
// Find all real roots in ] a[0] ; a[1] [ of a real
// d-th degree polynomial using the Sturm sequence.
// intervalType specifies as follows (in binary)
// 0 = 00 = ]a,b[
// 1 = 10 = [a,b[
// 2 = 01 = ]a,b]
// 3 = 11 = [a,b]
int vtkPolynomialSolversUnivariate::SturmBisectionSolve(
  double* P, int d, double* a, double* upperBnds, double tol)
{
  return vtkHabichtOrSturmBisectionSolve(P, d, a, upperBnds, tol, 0, 0, 0);
}

int vtkPolynomialSolversUnivariate::SturmBisectionSolve(
  double* P, int d, double* a, double* upperBnds, double tol, int intervalType)
{
  return vtkHabichtOrSturmBisectionSolve(P, d, a, upperBnds, tol, intervalType, 0, 0);
}

int vtkPolynomialSolversUnivariate::SturmBisectionSolve(
  double* P, int d, double* a, double* upperBnds, double tol, int intervalType, bool divideGCD)
{
  return vtkHabichtOrSturmBisectionSolve(
    P, d, a, upperBnds, tol, intervalType, divideGCD ? 1 : 0, 0);
}

// Assume that dP = {f} and p is the degree of f.
// Furthermore assume that dP is large enough.
// Stores
// {f,f',f^(2)/2!,f^(3)/3!,...,f^(p)/p!}.
static void vtkGetDerivativeSequence(double* dP, int p)
{
  int offsetA = 0;
  int offsetB = p + 1;

  for (int i = 1; i <= p; ++i)
  {
    for (int j = 0; j <= p - i; ++j)
      dP[offsetB + j] = static_cast<double>(p - i - j + 1) * dP[offsetA + j] / i;

    offsetA = offsetB;
    offsetB += p - i + 1;
  }
}

static int vtkGetSignChangesForDerivativeSequence(double* dP, int count, double val)
{
  int oldVal = 0;
  double v;
  int changes = 0;
  int offset = 0;

  for (int i = 0; i <= count; ++i)
  {
    v = evaluateHorner(dP + offset, count - i, val);

    if (v * oldVal < 0)
    {
      ++changes;
      oldVal = -oldVal;
    }
    if (oldVal == 0)
    {
      if (v < 0)
      {
        oldVal = -1;
      }
      else
      {
        oldVal = 1;
      }
    }
    offset += count - i + 1;
  }

  return changes;
}

int vtkPolynomialSolversUnivariate::FilterRoots(
  double* P, int d, double* upperBnds, int rootcount, double diameter)
{
  // Sort the roots.
  qsort(upperBnds, rootcount, sizeof(double), vtkPolynomialSolversUnivariateCompareRoots);

  // Remove duplicates.
  for (int j = 1; j < rootcount; ++j)
  {
    if (upperBnds[j] < upperBnds[j - 1] + diameter)
    {
      for (int k = j + 1; k < rootcount; ++k)
      {
        upperBnds[k - 1] = upperBnds[k];
      }
      --j;
      --rootcount;
      continue;
    }
  }

  if (rootcount == 0)
  {
    return 0;
  }

  // Ignore 0 as a root.
  for (int i = d; i >= 0; --i)
  {
    if (IsZero(P[i]))
    {
      --d;
    }
    else
    {
      break;
    }
  }

  double* dp = new double[((d + 2) * (d + 1)) / 2];
  for (int i = 0; i <= d; ++i)
  {
    dp[i] = P[i];
  }

  vtkGetDerivativeSequence(dp, d);

  // I could tag now and remove later. It'd be faster when there are a lot of
  // roots. But I doubt any of these root finders would work on a system with
  // enough roots for that method to be any better.
  for (int i = 0; i < rootcount; ++i)
  {
    if (fabs(upperBnds[i]) < diameter)
    {
      continue;
    }

    if (vtkGetSignChangesForDerivativeSequence(dp, d, upperBnds[i]) ==
      vtkGetSignChangesForDerivativeSequence(dp, d, upperBnds[i] - diameter))
    {
      // Remove the root.
      for (int j = i + 1; j < rootcount; ++j)
      {
        upperBnds[j - 1] = upperBnds[j];
      }
      --i;
      --rootcount;
    }
  }
  delete[] dp;
  return rootcount;
}

//------------------------------------------------------------------------------
// Solves a d-th degree polynomial equation using Lin-Bairstow's method.
//
int vtkPolynomialSolversUnivariate::LinBairstowSolve(double* c, int d, double* r, double& tolerance)
{
  if (IsZero(c[0]))
  {
    vtkGenericWarningMacro(
      << "vtkPolynomialSolversUnivariate::LinBairstowSolve: Zero leading coefficient");
    return 0;
  }

  int i;
  int dp1 = d + 1;
  for (i = 1; i < dp1; ++i)
  {
    c[i] /= c[0];
  }

  double* div1 = new double[dp1];
  double* div2 = new double[dp1];
  div1[0] = div2[0] = 1;
  for (i = d; i > 2; i -= 2)
  {
    double det, detR, detS;
    double R = 0.;
    double S = 0.;
    double dR = 1.;
    double dS = 0.;
    int nIterations = 1;

    while ((fabs(dR) + fabs(dS)) > tolerance)
    {
      // relax tolerance after 100 iterations did not suffice to converge
      // within the current tolerance
      if (!(nIterations % 100))
      {
        R = vtkMath::Random(0., 2.);
        if (!(nIterations % 200))
        {
          tolerance *= 4.;
        }
      }

      div1[1] = c[1] - R;
      div2[1] = div1[1] - R;

      for (int j = 2; j <= i; ++j)
      {
        div1[j] = c[j] - R * div1[j - 1] - S * div1[j - 2];
        div2[j] = div1[j] - R * div2[j - 1] - S * div2[j - 2];
      }

      double u = div2[i - 1] * div2[i - 3];
      double v = div2[i - 2] * div2[i - 2];
      if (AreEqual(u, v, 1.e-6))
      {
        det = detR = detS = 1.;
      }
      else
      {
        det = u - v;
        detR = div1[i] * div2[i - 3] - div1[i - 1] * div2[i - 2];
        detS = div1[i - 1] * div2[i - 1] - div1[i] * div2[i - 2];
      }

      dR = detR / det;
      dS = detS / det;

      // prevent Jacobian from exploding faster than tolerance can be relaxed
      // by the means of a crude limiter
      if (fabs(dR) + fabs(dS) > 10.)
      {
        dR = vtkMath::Random(-1., 1.);
        dS = vtkMath::Random(-1., 1.);
      }

      R += dR;
      S += dS;
      ++nIterations;
    }

    for (int j = 0; j < i - 1; ++j)
    {
      c[j] = div1[j];
    }
    c[i] = S;
    c[i - 1] = R;
  }

  int nr = 0;
  for (i = d; i >= 2; i -= 2)
  {
    double delta = c[i - 1] * c[i - 1] - 4. * c[i];
    if (delta >= 0)
    {
      // check whether there are 2 simple roots or 1 double root
      if (delta != 0.0)
      {
        delta = sqrt(delta);
        // we have 2 simple real roots
        r[nr++] = (-c[i - 1] - delta) / 2.;
        // insert 2nd simple real root
        r[nr++] = (-c[i - 1] + delta) / 2.;
      }
      else
      {
        // we have a double real root
        r[nr++] = -c[1];
        r[nr++] = -c[1];
      }
    }
  }
  if ((d % 2) == 1)
  {
    // what's left when degree is odd
    r[nr++] = -c[1];
  }

  delete[] div1;
  delete[] div2;
  return nr;
}

//------------------------------------------------------------------------------
// Algebraically extracts REAL roots of the quartic polynomial with
// REAL coefficients X^4 + c[0] X^3 + c[1] X^2 + c[2] X + c[3]
// and stores them (when they exist) and their respective multiplicities
// in the r and m arrays.
int vtkPolynomialSolversUnivariate::FerrariSolve(double* c, double* r, int* m, double tol)
{
  // step 0: eliminate trivial cases up to numerical noise
  if (fabs(c[3]) <= tol)
  {
    if (fabs(c[2]) <= tol)
    {
      if (fabs(c[1]) <= tol)
      {
        if (fabs(c[0]) <= tol)
        {
          r[0] = 0.;
          m[0] = 4;
          return 1;
        }
        else
        {
          r[0] = -c[1];
          m[0] = 1;
          r[1] = 0.;
          m[1] = 3;
          return 2;
        }
      }
      else
      {
        double cc[3];
        cc[0] = 1.;
        cc[1] = c[0];
        cc[2] = c[1];
        int nr = vtkPolynomialSolversUnivariate::SolveQuadratic(cc, r, m);
        r[nr] = 0.;
        m[nr] = 2;
        return nr + 1;
      }
    }
    else
    {
      int nr = vtkPolynomialSolversUnivariate::TartagliaCardanSolve(c, r, m, tol);
      r[nr] = 0.;
      m[nr] = 1;
      return nr + 1;
    }
  }
  if ((fabs(c[0]) <= tol) && (fabs(c[2]) <= tol))
  {
    if (fabs(c[1]) <= tol)
    {
      if (c[3] < 0.)
      {
        return 0;
      }
      r[0] = sqrt(sqrt(c[3]));
      m[0] = 4;
      return 1;
    }
    double cc[3], cr[2];
    int cm[2];
    cc[0] = 1.;
    cc[1] = c[1];
    cc[2] = c[3];
    int nr1 = vtkPolynomialSolversUnivariate::SolveQuadratic(cc, cr, cm);
    int nr = 0;
    int i;
    for (i = 0; i < nr1; ++i)
    {
      if (fabs(cr[i]) <= tol)
      {
        r[nr] = 0.;
        m[nr++] = 2 * cm[i];
      }
      else
      {
        if (cr[i] > tol)
        {
          r[nr] = sqrt(cr[i]);
          m[nr++] = cm[i];
          r[nr] = -sqrt(cr[i]);
          m[nr++] = cm[i];
        }
      }
    }
    return nr;
  }

  // step 1: reduce to X^4 + aX^2 + bX + d
  double p2d8 = c[0] * c[0] * .125;
  double qd2 = c[1] * .5;
  double a = c[1] - 3 * p2d8;
  double b = c[0] * (p2d8 - qd2) + c[2];
  double d = p2d8 * (qd2 - .75 * p2d8) - c[0] * c[2] * .25 + c[3];
  // expedite the case when the reduced equation is biquadratic
  if (fabs(b) <= tol)
  {
    double cc[3], cr[2];
    int cm[2];
    cc[0] = 1.;
    cc[1] = a;
    cc[2] = d;
    int nr1 = vtkPolynomialSolversUnivariate::SolveQuadratic(cc, cr, cm);
    int nr = 0;
    double shift = -c[0] * .25;
    int i;
    for (i = 0; i < nr1; ++i)
    {
      if (fabs(cr[i]) <= tol)
      {
        r[nr] = shift;
        m[nr++] = 2 * cm[i];
      }
      else
      {
        if (cr[i] > tol)
        {
          r[nr] = sqrt(cr[i]) + shift;
          m[nr++] = cm[i];
          r[nr] = -sqrt(cr[i]) + shift;
          m[nr++] = cm[i];
        }
      }
    }
    return nr;
  }

  // step 2: solve the companion cubic
  double cc[3], cr[3];
  double unsorted[8];
  int cm[3];
  cc[0] = 2. * a;
  cc[1] = a * a - 4. * d;
  cc[2] = -b * b;
  int nr = vtkPolynomialSolversUnivariate::TartagliaCardanSolve(cc, cr, cm, tol);

  // step 3: figure alpha^2
  double alpha2 = cr[--nr];
  while (alpha2 < 0. && nr)
    alpha2 = cr[--nr];

  // step 4: solve the quadratics
  cc[0] = 1.;
  cc[1] = sqrt(alpha2);
  double rho = -b / cc[1];
  cc[2] = (a + alpha2 + rho) * .5;
  int nr1 = vtkPolynomialSolversUnivariate::SolveQuadratic(cc, r, m);
  cc[1] = -cc[1];
  cc[2] -= rho;
  nr = nr1 + vtkPolynomialSolversUnivariate::SolveQuadratic(cc, r + nr1, m + nr1);
  if (!nr)
  {
    return 0;
  }

  // step 5: sort, filter and shift roots (if any)
  int i;
  for (i = 0; i < nr; ++i)
  {
    unsorted[2 * i] = r[i];
    unsorted[2 * i + 1] = m[i];
  }
  qsort(unsorted, nr, 2 * sizeof(double), vtkPolynomialSolversUnivariateCompareRoots);
  r[0] = unsorted[0];
  m[0] = static_cast<int>(unsorted[1]);
  nr1 = 1;
  for (i = 1; i < nr; ++i)
  {
    if (unsorted[2 * i] == unsorted[2 * i - 2])
    {
      m[i - 1] += static_cast<int>(unsorted[2 * i + 1]);
      continue;
    }
    r[nr1] = unsorted[2 * i];
    m[nr1++] = static_cast<int>(unsorted[2 * i + 1]);
  }
  double shift = -c[0] * .25;
  for (i = 0; i < nr1; ++i)
  {
    r[i] += shift;
  }

  return nr1;
}

//------------------------------------------------------------------------------
// Algebraically extracts REAL roots of the cubic polynomial with
// REAL coefficients X^3 + c[0] X^2 + c[1] X + c[2]
// and stores them (when they exist) and their respective multiplicities.
// The main differences with SolveCubic are that (1) the polynomial must have
// unit leading coefficient, (2) no information is returned regarding complex
// roots, and (3) non-simple roots are stored only once -- this is a
// specialized solver.
// Returns the number of roots.
//
int vtkPolynomialSolversUnivariate::TartagliaCardanSolve(double* c, double* r, int* m, double tol)
{
  // step 0: eliminate trivial cases up to numerical noise
  if (fabs(c[2]) <= tol)
  {
    r[0] = 0.;
    if (fabs(c[1]) <= tol)
    {
      if (fabs(c[0]) <= tol)
      {
        m[0] = 3;
        return 1;
      }
      else
      {
        m[0] = 2;
        r[1] = -c[0];
        m[1] = 1;
        return 2;
      }
    }
    else
    {
      m[0] = 1;
      double a2 = c[0] * c[0];
      double fourc1 = 4. * c[1];
      double delta = a2 - fourc1;
      double threshold = tol * (a2 > fabs(fourc1) ? a2 : fabs(fourc1));
      if (delta > threshold)
      {
        delta = sqrt(delta);
        r[1] = (-delta - c[0]) * 0.5;
        m[1] = 1;
        r[2] = (delta - c[0]) * 0.5;
        m[2] = 1;
        return 3;
      }
      else
      {
        if (delta < -threshold)
        {
          return 1;
        }
        r[1] = -c[0] * 0.5;
        m[1] = 2;
        return 2;
      }
    }
  }

  // step 1: reduce to X^3 + pX + q
  double shift = -c[0] / 3.;
  double a2 = c[0] * c[0];
  double p = c[1] - a2 / 3.;
  double q = c[0] * (2. * a2 / 9. - c[1]) / 3. + c[2];

  // step 2: compute the trivial real roots if p or q are 0
  // case 2.1: p = 0: 1 triple real root
  if (fabs(p) <= tol)
  {
    if (fabs(q) <= tol)
    {
      r[0] = +shift;
      m[0] = 3;
      return 1;
    }
    double x;
    x = q < 0 ? pow(-q, inv3) : -pow(q, inv3);
    r[0] = x + shift;
    m[0] = 3;
    return 1;
  }

  // case 2.2: q = 0: 1 ( p > 0 ) or 3 ( p < 0 ) simple real root(s)
  if (fabs(q) <= tol)
  {
    r[0] = +shift;
    m[0] = 1;
    if (p < 0)
    {
      double x = sqrt(-p);
      r[1] = x + shift;
      r[2] = -x + shift;
      m[1] = m[2] = 1;
      return 3;
    }
    return 1;
  }

  // step 3: compute discriminant
  double p_3 = p * inv3;
  double q_2 = q * 0.5;
  double D = p_3 * p_3 * p_3 + q_2 * q_2;

  // step 4: compute roots depending on the discriminant
  double u;
  // 4.1: case D = 0: 1 simple and 1 double real roots
  if (fabs(D) <= tol)
  {
    u = q > 0 ? -pow(q_2, inv3) : pow(-q_2, inv3);
    r[0] = 2. * u + shift;
    m[0] = 1;
    r[1] = -u + shift;
    m[1] = 2;
    return 2;
  }
  // 4.2: case D > 0: 1 simple real root
  if (D > 0)
  {
    u = sqrt(D) - q_2;
    u = u < 0 ? -pow(-u, inv3) : pow(u, inv3);
    r[0] = u - p_3 / u + shift;
    m[0] = 1;
    return 1;
  }
  // 5.3: case D < 0: 3 simple real roots
  double smp_3 = sqrt(-p_3);
  double argu = acos(q_2 / (p_3 * smp_3)) * inv3;
  double x1 = cos(argu);
  double x2 = sqrt3 * sqrt(1. - x1 * x1);
  x1 *= smp_3;
  x2 *= smp_3;
  r[0] = 2. * x1 + shift;
  r[1] = x2 - x1 + shift;
  r[2] = r[1] - 2. * x2;
  m[0] = m[1] = m[2] = 1;
  return 3;
}

//------------------------------------------------------------------------------
// Solves a cubic equation c0*t^3  + c1*t^2  + c2*t + c3 = 0 when
// c0, c1, c2, and c3 are REAL.
// Solution is motivated by Numerical Recipes In C 2nd Ed.
// Return array contains number of (real) roots (counting multiple roots as one)
// followed by roots themselves. The value in roots[4] is a integer giving
// further information about the roots (see return codes for int SolveCubic()).
double* vtkPolynomialSolversUnivariate::SolveCubic(double c0, double c1, double c2, double c3)
{
  static double roots[5];
  roots[1] = 0.0;
  roots[2] = 0.0;
  roots[3] = 0.0;
  int num_roots;

  roots[4] = vtkPolynomialSolversUnivariate::SolveCubic(
    c0, c1, c2, c3, &roots[1], &roots[2], &roots[3], &num_roots);
  roots[0] = num_roots;
  return roots;
}

//------------------------------------------------------------------------------
// Solves a cubic equation when c0, c1, c2, And c3 Are REAL.  Solution
// is motivated by Numerical Recipes In C 2nd Ed.  Roots and number of
// real roots are stored in user provided variables r1, r2, r3, and
// num_roots. Note that the function can return the following integer
// values describing the roots: (0)-no solution; (-1)-infinite number
// of solutions; (1)-one distinct real root of multiplicity 3 (stored
// in r1); (2)-two distinct real roots, one of multiplicity 2 (stored
// in r1 & r2); (3)-three distinct real roots; (-2)-quadratic equation
// with complex conjugate solution (real part of root returned in r1,
// imaginary in r2); (-3)-one real root and a complex conjugate pair
// (real root in r1 and real part of pair in r2 and imaginary in r3).
int vtkPolynomialSolversUnivariate::SolveCubic(
  double c0, double c1, double c2, double c3, double* r1, double* r2, double* r3, int* num_roots)
{
  double Q, R;
  double R_squared; /* R*R */
  double Q_cubed;   /* Q*Q*Q */
  double theta;
  double A, B;

  // Cubic equation: c0*t^3  + c1*t^2  + c2*t + c3 = 0
  //
  //   r1, r2, r3 are roots and num_roots is the number
  //   of real roots

  // Make Sure This Is A Bonafide Cubic Equation
  if (c0 != 0.0)
  {
    // Put Coefficients In Right Form
    c1 = c1 / c0;
    c2 = c2 / c0;
    c3 = c3 / c0;

    Q = ((c1 * c1) - 3 * c2) / 9.0;

    R = (2.0 * (c1 * c1 * c1) - 9.0 * (c1 * c2) + 27.0 * c3) / 54.0;

    R_squared = R * R;
    Q_cubed = Q * Q * Q;

    if (R_squared <= Q_cubed)
    {
      if (Q_cubed == 0.0)
      {
        *r1 = -c1 / 3.0;
        *r2 = *r1;
        *r3 = *r1;
        *num_roots = 1;
        return 1;
      }
      else
      {
        theta = acos(R / (sqrt(Q_cubed)));

        *r1 = -2.0 * sqrt(Q) * cos(theta / 3.0) - c1 / 3.0;
        *r2 = -2.0 * sqrt(Q) * cos((theta + 2.0 * vtkMath::Pi()) / 3.0) - c1 / 3.0;
        *r3 = -2.0 * sqrt(Q) * cos((theta - 2.0 * vtkMath::Pi()) / 3.0) - c1 / 3.0;

        *num_roots = 3;

        // Reduce Number Of Roots To Two
        if (*r1 == *r2)
        {
          *num_roots = 2;
          *r2 = *r3;
        }
        else if (*r1 == *r3)
        {
          *num_roots = 2;
        }

        if ((*r2 == *r3) && (*num_roots == 3))
        {
          *num_roots = 2;
        }

        // Reduce Number Of Roots To One
        if (*r1 == *r2)
        {
          *num_roots = 1;
        }
      }
      return *num_roots;
    }
    else // single real and complex conjugate pair
    {
      A = -VTK_SIGN(R) * pow(fabs(R) + sqrt(R_squared - Q_cubed), 1.0 / 3);

      if (A == 0.0)
      {
        B = 0.0;
      }
      else
      {
        B = Q / A;
      }

      *r1 = (A + B) - c1 / 3.0;
      *r2 = -0.5 * (A + B) - c1 / 3.0;
      *r3 = sqrt(3.0) / 2.0 * (A - B);

      *num_roots = 1;
      return (-3);
    }
  } // if cubic equation

  else // Quadratic Equation: c1*t^2  + c2*t + c3 = 0
  {
    // Okay this was not a cubic - lets try quadratic
    return vtkPolynomialSolversUnivariate::SolveQuadratic(c1, c2, c3, r1, r2, num_roots);
  }
}

//------------------------------------------------------------------------------
// Solves a quadratic equation c0*t^2 + c1*t + c2 = 0 when c0, c1, and
// c2 are REAL.  Solution is motivated by Numerical Recipes In C 2nd
// Ed.  Return array contains number of (real) roots (counting
// multiple roots as one) followed by roots themselves. Note that
// roots[3] contains a return code further describing solution - see
// documentation for SolveCubic() for meaining of return codes.
double* vtkPolynomialSolversUnivariate::SolveQuadratic(double c0, double c1, double c2)
{
  static double roots[4];
  roots[0] = 0.0;
  roots[1] = 0.0;
  roots[2] = 0.0;
  int num_roots;

  roots[3] =
    vtkPolynomialSolversUnivariate::SolveQuadratic(c0, c1, c2, &roots[1], &roots[2], &num_roots);
  roots[0] = num_roots;
  return roots;
}

//------------------------------------------------------------------------------
// Solves A Quadratic Equation c0*t^2  + c1*t  + c2 = 0 when
// c0, c1, and c2 are REAL.
// Solution is motivated by Numerical Recipes In C 2nd Ed.
// Roots and number of roots are stored in user provided variables
// r1, r2, num_roots
int vtkPolynomialSolversUnivariate::SolveQuadratic(
  double c0, double c1, double c2, double* r1, double* r2, int* num_roots)
{
  double Q;
  double determinant;

  // Quadratic equation: c0*t^2 + c1*t + c2 = 0

  // Make sure this is a quadratic equation
  if (c0 != 0.0)
  {
    determinant = c1 * c1 - 4 * c0 * c2;

    if (determinant >= 0.0)
    {
      Q = -0.5 * (c1 + VTK_SIGN(c1) * sqrt(determinant));

      *r1 = Q / c0;

      if (Q == 0.0)
      {
        *r2 = 0.0;
      }
      else
      {
        *r2 = c2 / Q;
      }

      *num_roots = 2;

      // Reduce Number Of Roots To One
      if (*r1 == *r2)
      {
        *num_roots = 1;
      }
      return *num_roots;
    }
    else // Equation Does Not Have Real Roots
    {
      *num_roots = 0;
      return (-2);
    }
  }

  else // Linear Equation: c1*t + c2 = 0
  {
    // Okay this was not quadratic - lets try linear
    return vtkPolynomialSolversUnivariate::SolveLinear(c1, c2, r1, num_roots);
  }
}

//------------------------------------------------------------------------------
// Algebraically extracts REAL roots of the quadratic polynomial with
// REAL coefficients c[0] X^2 + c[1] X + c[2]
// and stores them (when they exist) and their respective multiplicities.
// Returns either the number of roots, or -1 if infinite number of roots.
int vtkPolynomialSolversUnivariate::SolveQuadratic(double* c, double* r, int* m)
{
  if (c[0] == 0.0)
  {
    if (c[1] != 0.0)
    {
      r[0] = -c[2] / c[1];
      m[0] = 1;
      return 1;
    }
    else
    {
      if (c[2] != 0.0)
        return 0;
      else
        return -1;
    }
  }

  double delta = c[1] * c[1] - 4. * c[0] * c[2];

  if (delta >= 0.0)
  {
    double fac = 1. / (2. * c[0]);
    // check whether there are 2 simple or 1 double root(s)
    if (delta != 0.0)
    {
      delta = sqrt(delta);
      // insert 1st simple real root
      r[0] = (-delta - c[1]) * fac;
      m[0] = 1;
      // insert 2nd simple real root
      r[1] = (delta - c[1]) * fac;
      m[1] = 1;
      return 2;
    }
    else
    {
      // insert single double real root
      r[0] = -c[1] * fac;
      m[0] = 2;
      return 1;
    }
  }
  else
  {
    return 0;
  }
}

//------------------------------------------------------------------------------
// Solves a linear equation c0*t  + c1 = 0 when c0 and c1 are REAL.
// Solution is motivated by Numerical Recipes In C 2nd Ed.
// Return array contains number of roots followed by roots themselves.
double* vtkPolynomialSolversUnivariate::SolveLinear(double c0, double c1)
{
  static double roots[3];
  int num_roots;
  roots[1] = 0.0;
  roots[2] = vtkPolynomialSolversUnivariate::SolveLinear(c0, c1, &roots[1], &num_roots);
  roots[0] = num_roots;
  return roots;
}

//------------------------------------------------------------------------------
// Solves a linear equation c0*t + c1 = 0 when c0 and c1 are REAL.
// Solution is motivated by Numerical Recipes In C 2nd Ed.
// Root and number of (real) roots are stored in user provided variables
// r2 and num_roots.
int vtkPolynomialSolversUnivariate::SolveLinear(double c0, double c1, double* r1, int* num_roots)
{
  // Linear equation: c0*t + c1 = 0
  // Now this had better be linear
  if (c0 != 0.0)
  {
    *r1 = -c1 / c0;
    *num_roots = 1;
    return *num_roots;
  }
  else
  {
    *num_roots = 0;
    if (c1 == 0.0)
    {
      return (-1);
    }
  }

  return *num_roots;
}

//------------------------------------------------------------------------------
void vtkPolynomialSolversUnivariate::SetDivisionTolerance(double tol)
{
  vtkPolynomialSolversUnivariate::DivisionTolerance = tol;
}

//------------------------------------------------------------------------------
double vtkPolynomialSolversUnivariate::GetDivisionTolerance()
{
  return vtkPolynomialSolversUnivariate::DivisionTolerance;
}
VTK_ABI_NAMESPACE_END