<commit_before>/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkPolynomialSolversUnivariate.cxx

  Copyright (c) <NAME>, <NAME>, <NAME>
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================
  Copyright 2007 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
  license for use of this work by or on behalf of the
  U.S. Government. Redistribution and use in source and binary forms, with
  or without modification, are permitted provided that this Notice and any
  statement of authorship are reproduced on all copies.

  Contact: <EMAIL>,<EMAIL>

=========================================================================*/
#include "vtkPolynomialSolversUnivariate.h"
#include "vtkObjectFactory.h"
#include "vtkDataArray.h"
#include "vtkMath.h"

#if defined(_MSC_VER) || defined(__BORLANDC__) || defined (__SUNPRO_C) || defined(__SUNPRO_CC)
# define fmax(a,b) ( (a) >= (b) ? (a) : (b) )
#endif

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

vtkStandardNewMacro(vtkPolynomialSolversUnivariate);

static const double sqrt3 = sqrt( static_cast<double>( 3. ) );
static const 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: " << this->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 ) ? true : false;
}

//----------------------------------------------------------------------------
// 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 ? false : true;
    }
  else
    {
    return delta > rTol * absy ? false : true;
    }
}

//----------------------------------------------------------------------------
// 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 = 0 )
{
  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 = 0 )
{
  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 ofsets 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 k = 0;

  int deg = d;
  int degree = d - 1;
  int jp1 = 1;
  int ip1 = 0;
  while ( degree > 0  && j < d - 1 )
    {
    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;
}

extern "C" {

static int vtkPolynomialSolversUnivariateCompareRoots( const void* a, const void* b )
{
  return ( *( (const double*) a ) ) < ( *( (const double*) b ) ) ? -1 : 1;
}

} // extern "C"


// ------------------------------------------------------------
// 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() );
    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;
        }<commit_msg>Suppress warning about dead store<commit_after>/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkPolynomialSolversUnivariate.cxx

  Copyright (c) <NAME>, <NAME>, <NAME>
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================
  Copyright 2007 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
  license for use of this work by or on behalf of the
  U.S. Government. Redistribution and use in source and binary forms, with
  or without modification, are permitted provided that this Notice and any
  statement of authorship are reproduced on all copies.

  Contact: <EMAIL>,<EMAIL>

=========================================================================*/
#include "vtkPolynomialSolversUnivariate.h"
#include "vtkObjectFactory.h"
#include "vtkDataArray.h"
#include "vtkMath.h"

#if defined(_MSC_VER) || defined(__BORLANDC__) || defined (__SUNPRO_C) || defined(__SUNPRO_CC)
# define fmax(a,b) ( (a) >= (b) ? (a) : (b) )
#endif

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

vtkStandardNewMacro(vtkPolynomialSolversUnivariate);

static const double sqrt3 = sqrt( static_cast<double>( 3. ) );
static const 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: " << this->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 ) ? true : false;
}

//----------------------------------------------------------------------------
// 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 ? false : true;
    }
  else
    {
    return delta > rTol * absy ? false : true;
    }
}

//----------------------------------------------------------------------------
// 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 = 0 )
{
  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 = 0 )
{
  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 ofsets 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 k = 0;

  int deg = d;
  int degree = d - 1;
  int jp1 = 1;
  int ip1 = 0;
  while ( degree > 0  && j < d - 1 )
    {
    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;
}

extern "C" {

static int vtkPolynomialSolversUnivariateCompareRoots( const void* a, const void* b )
{
  return ( *( (const double*) a ) ) < ( *( (const double*) b ) ) ? -1 : 1;
}

} // extern "C"


// ------------------------------------------------------------
// 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;
        }
{'fraction_non_alphanumeric': 0.19334430614134349, 'fraction_numerical': 0.01846517847958029, 'mean_word_length': 1.6531795851528384, 'pattern_counts': {'":': 8, '<': 263, '<?xml version=': 0, '>': 85, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3425369', 'n_tokens_mistral': 14754, 'n_tokens_neox': 13826, 'n_words': 4520}