text: /*
   * Copyright 2011 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.
   */
  #include "vtkMath.h"
  #include "vtkPolynomialSolversUnivariate.h"
  #include "vtkTimerLog.h"
  
  static int vtkRunPolynomial(
    double* poly, int degree, double* rootInt, double* upperBnds, 
    double tolSturm, double divtol, double* expected, int expectedLength, 
    double expectedTol, const char* name, bool divideGCD, bool useHabichtSolver, int intType)
  {
    int rootcount;
    const char* sname = useHabichtSolver ? "Habicht" : "Sturm";
    vtkTimerLog* timer = vtkTimerLog::New();
  
    vtkPolynomialSolversUnivariate::SetDivisionTolerance( divtol );
    if ( useHabichtSolver )
      {
      timer->StartTimer();
      rootcount = vtkPolynomialSolversUnivariate::HabichtBisectionSolve( 
        poly, degree, rootInt, upperBnds, tolSturm, intType, divideGCD );
      timer->StopTimer();
      }
    else
      {
      timer->StartTimer();
      rootcount = vtkPolynomialSolversUnivariate::SturmBisectionSolve( 
        poly, degree, rootInt, upperBnds, tolSturm, intType, divideGCD );
      timer->StopTimer();
      }
  
    int trueexp = expectedLength;
    cout << "divtol is: " << divtol << ", " << timer->GetElapsedTime() << "s" << endl;
    timer->Delete();
  
    int j = 0;
    for(int i = 0; i < rootcount; i++, j++)
      {
      cout << upperBnds[i] << endl;
      while(j < expectedLength && 
          (
            ( ((intType & 1) == 0) && expected[j] == rootInt[0] ) ||
            ( ((intType & 2) == 0) && expected[j] == rootInt[1] ) ||
            expected[j] < rootInt[0] ||
            expected[j] > rootInt[1]
          ))
        {
        j++;
        trueexp--;
        }
      if ( j >= expectedLength )
        {
        vtkGenericWarningMacro( <<  sname << "BisectionSolve( " << name << ", " << ((intType&1)?"[":"]") << rootInt[0] << ", " << rootInt[1] << " " << ((intType&2)?"]":"[") << " )," 
          << upperBnds[i] << " found but not expected." << "." );
        return 1;
        }
      else if ( fabs( upperBnds[i] - expected[j] ) > expectedTol)
        {
        vtkGenericWarningMacro( << sname << "BisectionSolve( " << name << ", " << ((intType&1)?"[":"]") << rootInt[0] << ", " << rootInt[1] << " " << ((intType&2)?"]":"[") << " ),"
          << " found: "
          << upperBnds[i] << ", expected " << expected[j] << ".");
        return 1;
        }
      }
    while(j < expectedLength)
      {
      if ( expected[j] <= rootInt[0] || expected[j] >= rootInt[1] )
        {
        trueexp--;
        }
      j++;
      }
    if( rootcount != trueexp )
      {
      vtkGenericWarningMacro( << sname << "BisectionSolve( " << name << ", " << ((intType&1)?"[":"]") << rootInt[0] << ", " << rootInt[1] << " " << ((intType&2)?"]":"[") << " ), found: "
        << rootcount << " roots, expected " << trueexp << " roots.");
      return 1;
      }
  
    return 0;
  }
  
  //=============================================================================
  static int vtkTestPolynomials(
    double * poly, int degree, double* rootInt, double* upperBnds, 
    double tolSturm, double* divtols, int len, double* expectd, int expectedLength,
    double expectTol, const char* name, bool divideGCD, int methods = 3, int intType = 0 )
  {
    int rval = 0;
    if( methods & 1 )
      {
      cout << endl << name << " (Sturm)" << endl;
      for ( int i = 0; i < len; ++ i )
        {
        rval |= vtkRunPolynomial(poly, degree, rootInt, upperBnds, tolSturm, divtols[i], 
          expectd, expectedLength, expectTol, name, divideGCD, false, intType );
        }
      }
  
    if( methods & 2 )
      {
      cout << endl << name << " (Habicht)" << endl;
      for ( int i = 0; i < len; ++ i )
        {
        rval |= vtkRunPolynomial(poly, degree, rootInt, upperBnds, tolSturm, divtols[i], 
          expectd, expectedLength, expectTol, name, divideGCD, true, intType );
        }
      }
    return rval;
  }
  
  //=============================================================================
  int TestPolynomialSolversUnivariate( int len, char * c[] )
  {
    (void) len;
    (void) c;
    double tolSturm = 1.e-18;
    double tolRoots = 1.e-15;
    double roots[3000];
    double rootInt[] = { -14., 28.1 };
    double divtol[] = { 1e-5,1e-6, 1e-7, 1e-8, 1e-9, 1e-10, 1e-11, 1e-12, 1e-13 };
    int limit = 6;
    int stat = 0;
    vtkTimerLog* timer = vtkTimerLog::New();
    timer->AllocateLog();
  
    cout << "------Testing set 'Jenkins Traub'------" << endl;
    cout << 13 << " polynomials in the set." << endl;
    // Testing mul(x-i, i=1..19)
    double JT2[] = {
      1.000000000000000e+00, -1.900000000000000e+02, 1.681500000000000e+04, 
      -9.205500000000000e+05, 3.491694600000000e+07, -9.739419000000000e+08, 
      2.069293363000000e+10, -3.422525119000000e+11, 4.465226757381000e+12, 
      -4.628064775191000e+13, 3.819220555021950e+14, -2.503858755467550e+15, 
      1.295363698994390e+16, -5.226090336251272e+16, 1.614297365301190e+17, 
      -3.713847873452280e+17, 6.101160757404918e+17, -6.686097303411533e+17, 
      4.315651468176384e+17, -1.216451004088320e+17
    };
    double JT2RTS[] = {
      1.000000000000000e+00, 2.000000000000000e+00, 3.000000000000000e+00, 
      4.000000000000000e+00, 5.000000000000000e+00, 6.000000000000000e+00, 
      7.000000000000000e+00, 8.000000000000000e+00, 9.000000000000000e+00, 
      1.000000000000000e+01, 1.100000000000000e+01, 1.200000000000000e+01, 
      1.300000000000000e+01, 1.400000000000000e+01, 1.500000000000000e+01, 
      1.600000000000000e+01, 1.700000000000000e+01, 1.800000000000000e+01, 
      1.900000000000000e+01
    };
    stat |= vtkTestPolynomials(JT2, 19, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT2RTS, 19, 2.500000e-03, "mul(x-i, i=1..19)", false, 3, 0);
  
    double JT22[] = {
      1., -171., 13566., -662796., 22323822., -549789282., 10246937272., 
      -147560703732., 1661573386473., -14710753408923., 102417740732658., 
      -557921681547048., 2.35312504054998e+015, -7.55152759206302e+015, 
      1.79507122809215e+016, -3.03212540077194e+016, 3.40122495938227e+016, 
      -2.23769880585216e+016, 6.402373705728e+015
    };
    stat |= vtkTestPolynomials(JT22, 18, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT2RTS, 18, 0.500000e-03, "mul(x-i, i=1..19)", false, 3, 0);
  
    double JT23[] = {
     1., -153., 10812., -468180., 13896582., -299650806., 4853222764., 
     -60202693980., 577924894833., -4308105301929., 24871845297936., 
     -110228466184200., 369012649234384., -909299905844112., 1.58331397572749e+015, 
     -1.82160244462464e+015, 1.2234055905792e+015, -355687428096000. 
    };
    stat |= vtkTestPolynomials(JT23, 17, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT2RTS, 17, 1e-4, "mul(x-i, i=1..19)", false, 3, 0);
  
  
  
    // Testing mul(x-10^(-i),i=1..8)
    double JT3[] = {
      1.000000000000000e+00, -1.111111100000000e-01, 1.122334332211000e-03, 
      -1.123456666543211e-06, 1.123557787755321e-10, -1.123456666543211e-15, 
      1.122334332211000e-21, -1.111111100000000e-28, 1.000000000000000e-36
    };
    double JT3RTS[] = {
      1.000000000000000e-08, 1.000000000000000e-07, 1.000000000000000e-06, 
      1.000000000000000e-05, 1.000000000000000e-04, 1.000000000000000e-03, 
      1.000000000000000e-02, 1.000000000000000e-01
    };
    stat |= vtkTestPolynomials(JT3, 8, rootInt, roots,
                  1e-9, divtol+1, limit, 
                  JT3RTS, 8, 1.000000e-08, "mul(x-10^(-i),i=1..8)", false, 3, 0);
  
    // Testing (x-1/10)^3*(x-5/10)*(x-6/10)*(x-7/10)
    double JT4[] = {
      1.000000000000000e+00, -2.100000000000000e+00, 1.640000000000000e+00, 
      -5.860000000000000e-01, 9.690000000000000e-02, -7.370000000000000e-03, 
      2.100000000000000e-04
    };
    double JT4RTS[] = {
      1.000000000000000e-01, 5.000000000000000e-01, 6.000000000000000e-01, 
      7.000000000000000e-01
    };
    stat |= vtkTestPolynomials(JT4, 6, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT4RTS, 4, 1.000000e-07, "(x-1/10)^3*(x-5/10)*(x-6/10)*(x-7/10)", true, 3, 0);
  
    // Testing (x-1/10)^4*(x-2/10)^3*(x-3/10)^2*(x-4/10)
    double JT5[] = {
      1.000000000000000e+00, -2.000000000000000e+00, 1.750000000000000e+00, 
      -8.820000000000000e-01, 2.835000000000000e-01, -6.072000000000000e-02, 
      8.777000000000000e-03, -8.458000000000000e-04, 5.204000000000000e-05, 
      -1.848000000000000e-06, 2.880000000000000e-08
    };
    double JT5RTS[] = {
      1.000000000000000e-01, 2.000000000000000e-01, 3.000000000000000e-01, 
      4.000000000000000e-01
    };
    stat |= vtkTestPolynomials(JT5, 10, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT5RTS, 4, 1.000000e-07, "(x-1/10)^4*(x-2/10)^3*(x-3/10)^2*(x-4/10)", true, 3, 0);
  
  
    // Testing 25*10^9*(x-1/10)*(x-1001/1000)*(x-998/1000)*(x-100002/100000)*(x-99999/100000)
    double JT62[] = {
      2.500000000000000e+10, -1.024752500000000e+11, 1.599232247450000e+11, 
      -1.149182194640050e+11, 3.496776469345551e+10, -2.497519974450501e+09
    };
    double JT62RTS[] = {
      1.000000000000000e-01, 9.980000000000000e-01, 9.999900000000000e-01, 
      1.000020000000000e+00, 1.001000000000000e+00
    };
    stat |= vtkTestPolynomials(JT62, 5, rootInt, roots,
                  tolSturm, divtol+5, limit-2, 
                  JT62RTS, 5, 1.000000e-05, "25*10^9*(x-1/10)*(x-1001/1000)*(x-998/1000)*(x-100002/100000)*(x-99999/100000)", false, 3, 0);
  
  #if !defined(_MSC_VER) || ( defined(_MSC_VER) && (_MSC_VER != 1310) )
    // Testing (x-1/10)*(x-1001/1000)*(x-998/1000)*(x-100002/100000)*(x-99999/100000)
    // This only works with the Habicht Sequence, and not on MSVC 7.1
    double JT[] = {
      JT62[0]/JT62[0], JT62[1]/JT62[0], JT62[2]/JT62[0],
      JT62[3]/JT62[0], JT62[4]/JT62[0], JT62[5]/JT62[0]
    };
    stat |= vtkTestPolynomials(JT, 5, rootInt, roots,
                  tolSturm, divtol+5, limit-2, 
                  JT62RTS, 5, 1.000000e-05, "(x-1/10)*(x-1001/1000)*(x-998/1000)*(x-100002/100000)*(x-99999/100000)", false, 2, 0);
  #endif
  
  
    // Testing (x+1)^5
    double JT8[] = {
      1.000000000000000e+00, 5.000000000000000e+00, 1.000000000000000e+01, 
      1.000000000000000e+01, 5.000000000000000e+00, 1.000000000000000e+00
    };
    double JT8RTS[] = {
      -1.000000000000000e+00
    };
    stat |= vtkTestPolynomials(JT8, 5, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT8RTS, 1, 1.000000e-07, "(x+1)^5", true, 3, 0);
  
    rootInt[0] = -1e14;
    // Testing (x-10^(-13))*(x+10^(13))
    double JT9[] = {
      1.000000000000000e+00, 1.000000000000000e+13, -1.000000000000000e+00
    };
    double JT9RTS[] = {
      -1.000000000000000e+13, 1.000000000000000e-13
    };
    stat |= vtkTestPolynomials(JT9, 2, rootInt, roots,
                  tolSturm, divtol, limit, 
                  JT9RTS, 2, 1.000000e-07, "(x-10^(-13))*(x+10^(13))", false, 3, 0);
  
    rootInt[0] = -14;
    rootInt[1] = 1e4;
    // Testing (x-10^3)*(x-1)*(x-10^(-3))
    double P10a[] = {
      1.000000000000000e+00, -1.001001000000000e+03, 1.001001000000000e+03, 
      -1.000000000000000e+00
    };
    double P10aRTS[] = {
      1.000000000000000e-03, 1.000000000000000e+00, 1.000000000000000e+03
    };
    stat |= vtkTestPolynomials(P10a, 3, rootInt, roots,
                  tolSturm, divtol, limit, 
                  P10aRTS, 3, 1.000000e-07, "(x-10^3)*(x-1)*(x-10^(-3))", false, 3, 0);
  
    rootInt[1] = 1e7;
    // Testing (x-10^6)*(x-1)*(x-10^(-6))
    double P10b[] = {
      1.000000000000000e+00, -1.000001000001000e+06, 1.000001000001000e+06, 
      -1.000000000000000e+00
    };
    double P10bRTS[] = {
      1.000000000000000e-06, 1.000000000000000e+00, 1.000000000000000e+06
    };
    stat |= vtkTestPolynomials(P10b, 3, rootInt, roots,
                  tolSturm, divtol+1, limit-1, 
                  P10bRTS, 3, 1.000000e-07, "(x-10^6)*(x-1)*(x-10^(-6))", false, 3, 0);
  
    rootInt[1] = 1e10;
    // Testing (x-10^9)*(x-1)*(x-10^(-9))
    double P10c[] = {
      1.000000000000000e+00, -1.000000001000000e+09, 1.000000001000000e+09, 
      -1.000000000000000e+00
    };
    double P10cRTS[] = {
      1.000000000000000e-09, 1.000000000000000e+00, 1.000000000000000e+09
    };
    stat |= vtkTestPolynomials(P10c, 3, rootInt, roots,
                  tolSturm, divtol+4, limit-2, 
                  P10cRTS, 3, 1.000000e-07, "(x-10^9)*(x-1)*(x-10^(-9))", false, 3, 0);
  
    rootInt[1] = 28.1;
    cout << "------Testing set 'Igarashi Ympa'------" << endl;
    cout << 2 << " polynomials in the set." << endl;
    // Testing (x-2.35)*(x-2.37)*(x-2.39)
    double igyp00[] = {
      1.000000000000000e+00, -7.110000000000000e+00, 1.685030000000000e+01, 
      -1.331110500000000e+01
    };
    double igyp00RTS[] = {
      2.350000000000000e+00, 2.370000000000000e+00, 2.390000000000000e+00
    };
    stat |= vtkTestPolynomials(igyp00, 3, rootInt, roots,
                  tolSturm, divtol, limit, 
                  igyp00RTS, 3, 1.000000e-07, "(x-2.35)*(x-2.37)*(x-2.39)", false, 3, 0);
  
    // Testing (x-2.35)^3*(x-2.37)
    double igyp01[] = {
      1.000000000000000e+00, -9.420000000000000e+00, 3.327600000000000e+01, 
      -5.224285000000000e+01, 3.075756375000000e+01
    };
    double igyp01RTS[] = {
      2.350000000000000e+00, 2.370000000000000e+00
    };
    stat |= vtkTestPolynomials(igyp01, 4, rootInt, roots,
                  tolSturm, divtol+1, limit, 
                  igyp01RTS, 2, 1.000000e-07, "(x-2.35)^3*(x-2.37)", true, 3, 0);
  
    cout << "------Testing set 'Iliev'------" << endl;
    cout << 2 << " polynomials in the set." << endl;
    // Testing (x-1)*(x+2)^2*(x+3)^3
    double iliev01[] = {
      1.000000000000000e+00, 1.200000000000000e+01, 5.400000000000000e+01, 
      1.040000000000000e+02, 4.500000000000000e+01, -1.080000000000000e+02, 
      -1.080000000000000e+02
    };
    double iliev01RTS[] = {
      -3.000000000000000e+00, -2.000000000000000e+00, 1.000000000000000e+00
    };
    stat |= vtkTestPolynomials(iliev01, 6, rootInt, roots,
                  tolSturm, divtol, limit, 
                  iliev01RTS, 3, 1.000000e-07, "(x-1)*(x+2)^2*(x+3)^3", true, 3, 0);
  
    // Testing (x-1)^2*(x-2)^4*(x-3)^6
    double iliev02[] = {
      1.000000000000000e+00, -2.800000000000000e+01, 3.560000000000000e+02, 
      -2.716000000000000e+03, 1.383800000000000e+04, -4.956400000000000e+04, 
      1.278520000000000e+05, -2.390760000000000e+05, 3.212730000000000e+05, 
      -3.021840000000000e+05, 1.885680000000000e+05, -6.998400000000000e+04, 
      1.166400000000000e+04
    };
    double iliev02RTS[] = {
      1.000000000000000e+00, 2.000000000000000e+00, 3.000000000000000e+00
    };
    stat |= vtkTestPolynomials(iliev02, 12, rootInt, roots,
                  tolSturm, divtol, limit-1, 
                  iliev02RTS, 3, 1.000000e-07, "(x-1)^2*(x-2)^4*(x-3)^6", true, 3, 0);
  
    cout << "------Testing other sets------" << endl;
    // High multiplicities
    // (x-1)^6(x-2)^6(x-3)^6
    double polymult[] = {1., -36., 606., -6336., 46095., -247716., 1018816., -3278016., 8361951., -17033580., 27767046., 
      -36128736., 37235521., -29981196., 18442620., -8362656., 2632176., -513216., 46656.};
    double expectedmult[] = {1.,2.,3.};
    stat |= vtkTestPolynomials(polymult, 18, rootInt, roots, tolSturm, divtol, limit-1, expectedmult, 3, 1e-7, 
      "(x-1)^6(x-2)^6(x-3)^6", true, 1);
  
    stat |= vtkTestPolynomials(polymult, 18, rootInt, roots, tolSturm, divtol, limit, expectedmult, 3, 1e-7, 
      "(x-1)^6(x-2)^6(x-3)^6", true, 2);
  
    // Zeng Polynomials
    double twin01[] = {1., -3.96, 7.1366, -7.736316, 5.61584481, -2.874527064, 1.0633461368, 
      -0.286299144, 0.055658703376, -0.007616429184, 0.00069613420032, -3.8146387968e-005, 
      9.475854336e-007};
    double twinroots[] = {.2,.39,.4};
  
    stat |= vtkTestPolynomials(twin01, 12, rootInt, roots, tolSturm, divtol, 3, twinroots, 3, 1e-5,
      "(x-.39)^4*(x-.4)^4*(x-.2)^4", true, 1);
  
    stat |= vtkTestPolynomials(twin01, 12, rootInt, roots, tolSturm, divtol, 3, twinroots, 3, 1e-5,
      "(x-.39)^4*(x-.4)^4*(x-.2)^4", true, 2);
  
    double toh06a[] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
  
    stat |= vtkTestPolynomials(toh06a, 20, rootInt, roots, tolSturm, divtol, 2, 0, 0, 1e-7,
      "1+x+x^2+\\cdots+x^20", false);
  
    stat |= vtkTestPolynomials(toh06a, 20, rootInt, roots, tolSturm, divtol+3, 2, 0, 0, 1e-7,
      "1+x+x^2+\\cdots+x^20", false, 1);
  
    stat |= vtkTestPolynomials(toh06a, 20, rootInt, roots, tolSturm, divtol+4, 1, 0, 0, 1e-7,
      "1+x+x^2+\\cdots+x^20", false, 2);
  
    cout << "Test non-Sequence solvers" << endl;
    // 1. find the roots of a degree 4 polynomial with a 1 double root (1) and 2
    // simple roots (2 and 3) using:
    // 1.a FerrariSolve
    // 1.b SturmBissectionSolve
    double tolLinBairstow = 1.e-12;
    double tolDirectSolvers = VTK_DBL_EPSILON;
    int mult[4];
    double P4[] = { 1., -7., 17., -17., 6. } ;
    vtkPolynomialSolversUnivariate::PrintPolynomial(cout, P4, 4 );
  
    // 1.a FerrariSolve
    timer->StartTimer();
    int testIntValue = vtkPolynomialSolversUnivariate::FerrariSolve( P4 + 1, roots, mult, tolDirectSolvers );
    timer->StopTimer();
  
    if ( testIntValue != 3 )
      {
      vtkGenericWarningMacro("FerrariSolve(x^4 -7x^3 +17x^2 -17 x +6 ) = "<<testIntValue<<" != 3");
      timer->Delete();
      return 1;
      }
    cout << "Ferrari tol=" << tolDirectSolvers
                 << ", " << testIntValue << " " 
                 << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      cout << roots[i]; 
      if ( mult[i] > 1 )
        {
        cout << "(" << mult[i] << ")";
        }
      cout << "\n";
      }
    double actualRoots[] = { 1., 2., 3. };
    int actualMult[] = { 2, 1, 1 };
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      if ( fabs ( roots[i] - actualRoots[i] ) > tolRoots )
        {
        vtkGenericWarningMacro("FerrariSolve(x^4 -7x^3 +17x^2 -17 x +6, ]-4;4] ): root "<<roots[i]<<" != "<<actualRoots[i]);
        timer->Delete();
        return 1;  
        }
      if ( mult[i] != actualMult[i] )
        {
        vtkGenericWarningMacro("FerrariSolve(x^4 -7x^3 +17x^2 -17 x +6, ]-4;4] ): multiplicity "<<mult[i]<<" != "<<actualMult[i]);
        timer->Delete();
        return 1;  
        }
      }
  
    // 1.b SturmBisectionSolve
    double P4RTS[] = {1.,2.,3.};
  
    stat |= vtkTestPolynomials(P4, 4, rootInt, roots, tolSturm, divtol, limit-1, 
                  P4RTS, 3, 1e-7, "(x-1)^2*(x-2)*(x-3)", true, 3);
  
    // 2. find the roots of a degree 5 polynomial with LinBairstowSolve
    double P5[] = { 1., -10., 35., -50., 24., 0. } ;
    vtkPolynomialSolversUnivariate::PrintPolynomial(cout, P5, 5 );
  
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::LinBairstowSolve( P5, 5, roots, tolLinBairstow );
    timer->StopTimer();
  
    if ( testIntValue != 5 )
      {
      vtkGenericWarningMacro("LinBairstowSolve(x^5 -10x^4 +35x^3 -50x^2 +24x ) = "<<testIntValue<<" != 5");
      timer->Delete();
      return 1;
      }
    cout << "LinBairstow tol=" << tolLinBairstow << ", " 
         << testIntValue << " " 
         << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) cout << roots[i] << "\n";
  
    // 3. find the roots of a quadratic trinomial with the BisectionSolvers
    double P2[] = { 1., -2., 1. };
    double P2RTS[] = {1.};
  
    stat |= vtkTestPolynomials(P2, 2, rootInt, roots, tolSturm, divtol, limit-1, 
                  P2RTS, 1, 1e-7, "(x-1)^2", true, 3);
  
    // 4. find the roots of a biquadratic trinomial with SturmBisectionSolve,
    // whose 2 double roots (-4 and 4) are also the bounds of the interval, thus
    // being a limiting case of Sturm's theorem, using:
    // 4.a FerrariSolve
    // 4.b SturmBisectionSolve
  
    double P4_2[] = { 1., 0., -32., 0., 256. };
    vtkPolynomialSolversUnivariate::PrintPolynomial(cout, P4_2, 4 );
  
    // 4.a FerrariSolve
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::FerrariSolve( P4_2 + 1, roots, mult, tolDirectSolvers );
    timer->StopTimer();
  
    if ( testIntValue != 2 )
      {
      vtkGenericWarningMacro("FerrariSolve(x^4 -32x^2 +256 ) = "<<testIntValue<<" != 2");
      timer->Delete();
      return 1;
      }
    cout << "Ferrari tol=" << tolDirectSolvers
                 << ", " << testIntValue << " " 
                 << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      cout << roots[i]; 
      if ( mult[i] > 1 ) cout << "(" << mult[i] << ")";
      cout << "\n";
      }
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      if ( fabs ( roots[i] ) - 4. > tolRoots )
        {
        vtkGenericWarningMacro("FerrariSolve(1*x**4-32*x**2+256, ]-4;4] ): root "<<roots[i]<<" != +/-4");
        timer->Delete();
        return 1;  
        }
      if ( mult[i] != 2 )
        {
        vtkGenericWarningMacro("FerrariSolve(1*x**4-32*x**2+256, ]-4;4] ): multiplicity "<<mult[i]<<" != 2");
        timer->Delete();
        return 1;  
        }
      }
  
    // 4.b SturmBisectionSolve
    rootInt[0] = -4.;
    rootInt[1] = 4.;
    double P4_2RTS[] = {-4.,4.};
    stat |= vtkTestPolynomials(P4_2, 4, rootInt, roots, tolSturm, divtol, limit-1, 
                  P4_2RTS, 2, 1e-7, "(x+4)^2*(x-4)^2 on [-4,4]", true, 3, 3);
    stat |= vtkTestPolynomials(P4_2, 4, rootInt, roots, tolSturm, divtol, limit-1, 
                  P4_2RTS, 2, 1e-7, "(x+4)^2*(x-4)^2 on [-4,4)", true, 3, 1);
    stat |= vtkTestPolynomials(P4_2, 4, rootInt, roots, tolSturm, divtol, limit-1, 
                  P4_2RTS, 2, 1e-7, "(x+4)^2*(x-4)^2 on (-4,4]", true, 3, 2);
    stat |= vtkTestPolynomials(P4_2, 4, rootInt, roots, tolSturm, divtol, limit-1, 
                  P4_2RTS, 2, 1e-7, "(x+4)^2*(x-4)^2 on (-4,4)", false, 3, 0);
  
    // 5. find the quadruple roots of the degree 12 polynomial (x-1)^4 (x-2)^4 (x-3)^4
    // All roots are quadruple roots, making it challenging for solvers using floating
    // point arithmetic.
    rootInt[0] = 0.;
    rootInt[1] = 20.;
    double P12[] = {
      1,
      -24,
      260,
      -1680,
      7206,
      -21600,
      46364,
      -71760,
      79441,
      -61320,
      31320,
      -9504,
      1296
    };
    double P12RTS[] = {1.,2.,3.};
  
     stat |= vtkTestPolynomials(P12, 12, rootInt, roots, tolSturm, divtol, limit, 
                  P12RTS, 3, 1e-7, "(x-1)^4*(x-2)^4*(x-3)^4 on (0,20)", true, 3, 0);
  
    // 6. Find the roots of a degree 22 polynomial with SturmBisectionSolve
    rootInt[0] = -10.;
    rootInt[1] = 10.;
    double P22[] = {
      -0.0005, -0.001, 0.05, 0.1, -0.2,
      1., 0., -5.1, 0., 4., 
      -1., .2, 3., 2.2, 2.,
      -7., -.3, 3.8, 14., -16.,
      80., -97.9, 5. };
    double P22RTS[] = {-9.5799184021307155490, -4.1659457357018254697, -1.6909764377051033080, 0.053377023199573159218, 1.1281470811227336526, 1.2804679030668921769, 1.7510169549495913088, 9.9112562522641383112};
  
     stat |= vtkTestPolynomials(P22, 22, rootInt, roots, tolSturm, divtol, limit, 
                  P22RTS, 8, 1e-7, "-0.0005*x^22 -0.001*x^21 +0.05*x^20 +0.1*x^19 -0.2*x^18 +1*x^17 -5.1*x^15 +4*x^13 -1*x^12 +0.2*x^11 +3*x^10 +2.2*x^9 +2*x^8 -7*x^7 -0.3*x^6 +3.8*x^5 +14*x^4 -16*x^3 +80*x^2 -97.9*x +5", false, 3, 0);
  
  
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::LinBairstowSolve( P22, 22, roots, tolLinBairstow );
    timer->StopTimer();
  
    if ( testIntValue != 8 )
      {
      vtkGenericWarningMacro("LinBairstowSolve( -0.0005x^22 -0.001x^21 +0.05x^20 +0.1x^19 -0.2x^18 +1x^17 -5.1x^15 +4x^13 -1x^12 +0.2x^11 +3x^10 +2.2x^9 +2x^8 -7x^7 -0.3x^6 +3.8x^5 +14x^4 -16x^3 +80x^2 -97.9x +5, ]-10;10] ): "<<testIntValue<<" root(s) instead of 8");
      timer->Delete();
      return 1;
      }
    cout << "LinBairstow tol=" << tolLinBairstow << ", " 
         << testIntValue << " " 
         << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) cout << roots[i] << "\n";
  
    // 7. Solving x^4 + 3x^3 - 4x + 1e-18 = 0 illustrates how the Ferrari solver
    // filters some numerical noise by noticing there is a double root.
    // This also exercises a case not otherwise tested.
    double P4_3[] = { 1., 3., -4., 0., 1.e-18 };
  
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::FerrariSolve( P4_3 + 1, roots, mult, tolDirectSolvers );
    timer->StopTimer();
  
    if ( testIntValue != 3 )
      {
      vtkGenericWarningMacro("FerrariSolve(x^4 +3x^3 -3x^2 +1e-18 ) = "<<testIntValue<<" != 3");
      timer->Delete();
      return 1;
      }
    cout << "Ferrari tol=" << tolDirectSolvers
                 << ", " << testIntValue << " " 
                 << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      cout << roots[i]; 
      if ( mult[i] > 1 ) cout << "(" << mult[i] << ")";
      cout << "\n";
      }
  
    // 8. Solving x(x - 10^-4)^2 = 0 illustrates how the Tartaglia-Cardan solver
    // filters some numerical noise by noticing there is a double root (that
    // SolveCubic does not notice).
    double P3[] = { 1., -2.e-4, 1.e-8, 0.};
  
  #if 0
    double r1, r2, r3;
    int nr;
    testIntValue = vtkPolynomialSolversUnivariate::SolveCubic( P3[0], P3[1], P3[2], P3[3], &r1, &r2, &r3, &nr );
    if ( testIntValue != 2 )
      {
      vtkGenericWarningMacro("SolveCubic returned "<<testIntValue<<" != 3");
      timer->Delete();
      return 1;
      }
  #endif // 0  
  
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::TartagliaCardanSolve( P3 + 1, roots, mult, tolDirectSolvers );
    timer->StopTimer();
  
    if ( testIntValue != 2 )
      {
      vtkGenericWarningMacro("TartagliaCardanSolve returned "<<testIntValue<<" != 2");
      timer->Delete();
      return 1;
      }
    cout << "TartagliaCardan tol=" << tolDirectSolvers
                 << ", " << testIntValue << " " 
                 << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      cout << roots[i]; 
      if ( mult[i] > 1 ) cout << "(" << mult[i] << ")";
      cout << "\n";
      }
  
    // 9. Solving x^3+x^2+x+1 = 0 to exercise a case not otherwise tested.
    double P3_2[] = { 1., 1., 1., 1.};
  
  #if 0
    double r1, r2, r3;
    int nr;
    testIntValue = vtkPolynomialSolversUnivariate::SolveCubic( P3_2[0], P3_2[1], P3_2[2], P3_2[3], &r1, &r2, &r3, &nr );
    if ( testIntValue != 1 )
      {
      vtkGenericWarningMacro("SolveCubic returned "<<testIntValue<<" != 1");
      timer->Delete();
      return 1;
      }
  #endif // 0  
  
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::TartagliaCardanSolve( P3_2 + 1, roots, mult, tolDirectSolvers );
    timer->StopTimer();
  
    if ( testIntValue != 1 )
      {
      vtkGenericWarningMacro("TartagliaCardanSolve returned "<<testIntValue<<" != 1");
      timer->Delete();
      return 1;
      }
    cout << "TartagliaCardan tol=" << tolDirectSolvers
                 << ", " << testIntValue << " " 
                 << timer->GetElapsedTime() << "s\n";
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      cout << roots[i]; 
      if ( mult[i] > 1 ) cout << "(" << mult[i] << ")";
      cout << "\n";
      }
  
    // 10. Solving x^3 - 2e-6 x^2 + 0.999999999999999e-12 x = 0 to test a nearly degenerate case.
    double P3_3[] = { 1., -2.e-6 , .999999999999999e-12, 0.};
  
  #if 0
    double r1, r2, r3;
    int nr;
    testIntValue = vtkPolynomialSolversUnivariate::SolveCubic( P3_3[0], P3_3[1], P3_3[2], P3_3[3], &r1, &r2, &r3, &nr );
    if ( testIntValue != 3 )
      {
      vtkGenericWarningMacro("SolveCubic returned "<<testIntValue<<" != 3");
      timer->Delete();
      return 1;
      }
  #endif // 0  
  
    timer->StartTimer();
    testIntValue = vtkPolynomialSolversUnivariate::TartagliaCardanSolve( P3_3 + 1, roots, mult, tolDirectSolvers );
    timer->StopTimer();
  
    if ( testIntValue != 3 )
      {
      vtkGenericWarningMacro("TartagliaCardanSolve returned "<<testIntValue<<" != 3");
      timer->Delete();
      return 1;
      }
    cout << "TartagliaCardan tol=" << tolDirectSolvers
                 << ", " << testIntValue << " " 
                 << timer->GetElapsedTime() << "s\n";
    cout.precision( 9 );
    for ( int i = 0; i < testIntValue ; ++ i ) 
      {
      cout << roots[i]; 
      if ( mult[i] > 1 ) cout << "(" << mult[i] << ")";
      cout << "\n";
      }
  
    // 11. Find the roots of a sparse degree 10 polynomial with SturmBisectionSolve to exercise a particular
    // case of the Euclidean division routine, where the remainder does not have maximal degree. 
    rootInt[0] = -10.;
    rootInt[1] = 10.;
    double P10[] = { 76., 0., 0., 0., 0., 0., 0., 0., 95., 0., -14. };
    double P10RTS[] = {-0.38381365387613186613, 0.38381365386885590851};
  
     stat |= vtkTestPolynomials(P10, 10, rootInt, roots, tolSturm, divtol, limit, 
                  P10RTS, 2, 1e-7, "76*x^10 +95*x^2 -14", false, 3, 0);
  
  
    // 11. Find the roots of a sparse degree 10 polynomial with SturmBisectionSolve to exercise a particular
    // case of the Euclidean division routine, where the remainder does not have maximal degree. 
    rootInt[0] = -10.;
    rootInt[1] = 10.;
    double P84[] = 
      {
        55.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        -79.000000,
        90.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        37.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        49.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        0.000000,
        -12.000000
      };
    double P84RTS[] = {-1.0923995943740010262, 0.92872986290603876114};
  
    stat |= vtkTestPolynomials(P84, 84, rootInt, roots, tolSturm, divtol, limit, 
                  P84RTS, 2, 1e-7, "55*x^84-79*x^72+90*x^71+37*x^44+49*x^21-12", false, 3, 0);
  
  
    timer->Delete();
    return stat;
  }
  