#include "Racine.h"
#include "MathUtil.h"

#include <stdio.h>
#include <cmath>

#define NCASES 40

// constructeur permettant d'initialiser les variables
Quartic::Quartic():
      debug(6),iterate(false),
      doub0((double)0),
      doub1((double)1),
      doub2((double)2),
      doub3((double)3),
      doub4((double)4),
      doub5((double)5),
      doub6((double)6),
      doub8((double)8),
      doub9((double)9),
      doub12((double)12),
      doub16((double)16),
      doub24((double)24),
      doub27((double)27),
      doub64((double)64),
      inv2(doub1/doub2),
      inv3(doub1/doub3),
      inv4(doub1/doub4),
      inv8(doub1/doub8),
      inv16(doub1/doub16),
      inv32(doub1/(double)32),
      inv64(doub1/doub64),
      inv128(doub1/(double)128),
      d3o8(doub3/doub8),
      d3o256(doub3/(double)256),
      rt3(sqrt(doub3))
  {   doubtol = doub1;
      for (int j = 1; doub1+doubtol > doub1; ++ j)
      {
          doubtol *= inv2;
      }
      doubtol = sqrt(doubtol);

      doubmin = inv2;
      for (int j = 1; j <= 100; ++ j)
      {
          doubmin=doubmin*doubmin;
          if ((doubmin*doubmin) <= (doubmin*doubmin*inv2))
              break;
      }
      doubmax=doub1/sqrt(doub2*doubmin);
   };


/* ************************************************** */
//   find the errors
//
//   called by quartictest, docoeff, compare,
//            chris, descartes, ferrari, neumark, yacfraid.
/* ************************************************** */

double Quartic::errors(double a,double b,double c,double d,double rts[4],double rterr[4],int nrts)
{
   double deriv,test,worst;

   worst = doubmax;
   if (nrts > 0)
   {
      worst =  doub0;
      for (int   k = 0 ; k < nrts ; ++ k )
      {
         test = (((rts[k]+a)*rts[k]+b)*rts[k]+c)*rts[k]+d;
         if (test == doub0) rterr[k] = doub0;
         else
         {
            deriv =
               ((doub4*rts[k]+doub3*a)*rts[k]+doub2*b)*rts[k]+c;
            if (deriv != doub0)
               rterr[k] = Dabs(test/deriv);
            else
            {
               deriv = (doub12*rts[k]+doub6*a)*rts[k]+doub2*b;
               if (deriv != doub0)
                   rterr[k] = sqrt(Dabs(test/deriv));
               else
               {
                  deriv = doub24*rts[k]+doub6*a;
                  if (deriv != doub0)
                     rterr[k] = curoot(Dabs(test/deriv));
                  else
                     rterr[k] = sqrt(sqrt(Dabs(test)/doub24));
               }
            }
         }
         if (rts[k] != doub0) rterr[k] /= rts[k];
         if (rterr[k] < doub0) rterr[k] = -rterr[k];
         if (rterr[k] > worst) worst = rterr[k];
      }
   }
   return(worst);
}; /* errors */

/**********************************************/
//     find cos(acos(x)/3) 
//    
//     16 Jul 1981   Don Herbison-Evans 
//
//     called by cubic . 
/**********************************************/

double Quartic::acos3(double x)
{
   double value;
   double y = acos(x);
   
//   value = std::cos(acos(x)*inv3);  pb de nan
   // on traite dans le cas d'une erreur (qui ça arrive !!)
   // comme ça ne marche pas pour l'instant sur la lo-lg2m-045 on spécifie les machines !!!
  #ifdef SYSTEM_MAC_OS_X
   switch(fpclassify(y))
    { case FP_NAN: case FP_INFINITE:
        { // cas d'un nombre infini on fait quelques tests
          if (Dabs(x-1.) <= ConstMath::pasmalpetit)
          	{ // dans ce cas on a affaire à une erreur d'arrondi
          	  y=0.;
          	}
          else if (x > 1.)
          	{ // il y a un pb dans l'algo, on le signal mais on continu
          	  cout << "\n erreur dans la fonction Qartic::acos3, l'argument est > 1 , on limite ! "
          	       << " x-1= " << (x-1.);
          	  y=0.;
          	};
        };
        break;
    };
   #endif     
  #ifdef SYSTEM_MAC_OS_X_unix
   switch(fpclassify(y))
    { case FP_NAN: case FP_INFINITE:
        { // cas d'un nombre infini on fait quelques tests
          if (Dabs(x-1.) <= ConstMath::pasmalpetit)
          	{ // dans ce cas on a affaire à une erreur d'arrondi
          	  y=0.;
          	}
          else if (x > 1.)
          	{ // il y a un pb dans l'algo, on le signal mais on continu
          	  cout << "\n erreur dans la fonction Qartic::acos3, l'argument est > 1 , on limite ! "
          	       << " x-1= " << (x-1.);
          	  y=0.;
          	};
        };
        break;
    };
   #endif     
   // calcul de la valeur de retour
   value = std::cos(y*inv3);
   return(value);
}; /* acos3 */


/***************************************/
//     find cube root of x.
//
//     30 Jan 1989   Don Herbison-Evans 
//
//     called by cubic . 
/***************************************/

double Quartic::curoot(double x)
{
   double value;
   double absx;
   int neg;

   neg = 0;
   absx = x;
   if (x < doub0)
   {
      absx = -x;
      neg = 1;
   }
   if (absx != doub0) value = exp( log(absx)*inv3 );
      else value = doub0;
   if (neg == 1) value = -value;
   return(value);
}; /* curoot */

/****************************************************/
//     solve the quadratic equation - 
//
//         x**2 + b*x + c = 0 
//
//     14 Jan 2004   cut determinant in quadratic call
//     29 Nov 2003   improved
//     16 Jul 1981   Don Herbison-Evans
//
//     called by  cubic,quartic,chris,descartes,ferrari,neumark.
/****************************************************/

int Quartic::quadratic(double b,double c,double rts[4])
{
   int nquad;
   double dis,rtdis ;

   dis = b*b - doub4*c;
   rts[0] = doub0;
   rts[1] = doub0;
   if (b == doub0)
   {
      if (c == doub0)
      {
         nquad = 2;
         ++nqud[0];
      }
      else
      {
         if (c < doub0)
         {
            nquad = 2;
            rts[0] = sqrt(-c);
            rts[1] = -rts[0];
            ++nqud[1];
         }
         else
         {
            nquad = 0;
            ++nqud[2];
         }         
      }
   }
   else
   if (c == doub0)
   {
      nquad = 2;
      rts[0] = -b;
      ++nqud[3];
   }
   else
   if (dis >= doub0)
   {
      nquad = 2 ;
      rtdis = sqrt(dis);
      if (b > doub0) 
      {
         rts[0] = ( -b - rtdis)*inv2;
         ++nqud[4];
      }
      else
      {
         rts[0] = ( -b + rtdis)*inv2;
         ++nqud[5];
      }
      if (rts[0] == doub0)
      {
         rts[1] =  -b;
         ++nqud[6];
      }
      else
      {
         rts[1] = c/rts[0];
         ++nqud[7];
      }
   }
   else
   {
      nquad = 0;
      ++nqud[8];
   }
   if (debug < 1)
   {
      printf("quad  b %g   c %g  dis %g\n",
         b,c,dis);
      printf("      %d %g %g\n",
         nquad,rts[0],rts[1]);
   }
   return(nquad);
}; /* quadratic */

/**************************************************/
//   find the real roots of the cubic - 
//       x**3 + p*x**2 + q*x + r = 0 
//
//     12 Dec 2003 initialising n,m,po3
//     12 Dec 2003 allow return of 3 zero roots if p=q=r=0
//      2 Dec 2003 negating j if p>0
//      1 Dec 2003 changing v from (sinsqk > doub0) to (sinsqk >= doub0)
//      1 Dec 2003 test changing v from po3sq+po3sq to doub2*po3sq
//     16 Jul 1981 Don Herbison-Evans
//
//   input parameters - 
//     p,q,r - coeffs of cubic equation. 
//
//   output- 
//     the number of real roots
//     v3 - the roots. 
//
//   global constants -
//     rt3 - sqrt(3) 
//     inv3 - 1/3 
//     doubmax - square root of largest number held by machine 
//
//     method - 
//     see D.E. Littlewood, "A University Algebra" pp.173 - 6 
//
//     15 Nov 2003 output 3 real roots: Don Herbison-Evans
//        Apr 1981 initial version: Charles Prineas
//
//     called by  cubictest,quartic,chris,yacfraid,neumark,descartes,ferrari.
//     calls      quadratic,acos3,curoot,cubnewton. 
/**************************************************/

int Quartic::cubic(double p,double q,double r,double v3[4]) 
{
   int    j,n3;
   double po3,po3sq,qo3,po3q;
   double uo3,u2o3,uo3sq4,uo3cu4;
   double v,vsq,wsq;
   double m1,m2,mcube;
   double muo3,s,scube,t,cosk,rt3sink,sinsqk;

   m1=doub0;  m2=doub0;  po3=doub0;
   v=doub0;  uo3=doub0; cosk=doub0;
   if (r == doub0)
   {
      ++ncub[0];
      n3 = quadratic(p,q,v3);
      v3[n3++] = doub0;
      goto done;
   }
   if ((p == doub0) && (q == doub0))
   {
      ++ncub[1];
      v3[0] = curoot(-r);
      v3[1] = v3[0];
      v3[2] = v3[0];
      n3 = 3;
      goto done;
   }
   n3 = 1;
   if ((p > doubmax) || (p <  -doubmax))
   {
      v3[0] = -p;
      ++ncub[2];
      goto done;
   }
   if ((q > doubmax) || (q <  -doubmax))
   {
       if (q > doub0)
       {
          v3[0] =  -r/q;
          ++ncub[3];
          goto done;
       }
       else
       if (q < doub0)
       {
          v3[0] = -sqrt(-q);
          ++ncub[4];
          goto done; 
       }
       else
       {
          v3[0] = doub0;
          ++ncub[5];
          goto done;
       }
   }
   else
   if ((r > doubmax)|| (r < -doubmax))
   {
      v3[0] =  -curoot(r);
      ++ncub[6];
      goto done;
   }
   else
   {
      po3 = p*inv3;
      po3q = po3*q;
      po3sq = po3*po3;
      if (po3sq > doubmax)
      {
         v3[0] = -p;
         ++ncub[7];
         goto done;
      }
      else
      {
         v = r + po3*(po3sq+po3sq - q);
         if ((v > doubmax) || (v < -doubmax))
         {
            v3[0] = -p;
            ++ncub[8];
            goto done;
         }
         else
         {
            vsq = v*v;
            qo3 = q*inv3;
            uo3 = qo3 - po3sq;
            u2o3 = uo3 + uo3;
            if ((u2o3 > doubmax) || (u2o3 < -doubmax))
            {
               if (p == doub0)
               {
                  if (q > doub0)
                  {
                     v3[0] =  -r/q;
                     ++ncub[9];
                     goto done;
                  }
		      else
                  if (q < doub0)
                  {
                     v3[0] =  -sqrt(-q);
                     ++ncub[10];
                     goto done;
                  }
                  else
                  {
                     v3[0] = doub0;
                     ++ncub[11];
                     goto done;
                  }
               }
               else
               {
                  v3[0] = -q/p;
                  ++ncub[12];
                  goto done;
               }
            }
            uo3sq4 = u2o3*u2o3;
            if (uo3sq4 > doubmax)
            {
               if (p == doub0)
               {
                  if (q > doub0)
                  {
                     v3[0] = -r/q;
                     ++ncub[13];
                     goto done;
                  }
                  else
	  	      if (q < doub0)
                  {
                     v3[0] = -sqrt(-q);
                     ++ncub[14];
                     goto done;
                  }
		      else
                  {
                     v3[0] = doub0;
                     ++ncub[15];
                     goto done;
                  }
               }
               else
               {
                  v3[0] = -q/p;
                  ++ncub[16];
                  goto done;
               }
            }
            uo3cu4 = uo3sq4*uo3;
            wsq = uo3cu4 + vsq;
            if (wsq > doub0)
            {
/* 
     cubic has one real root -
*/
               if (v <= doub0)
               {
                  mcube = ( -v + sqrt(wsq))*inv2;
                  ++ncub[17];
               }
               else
               {
                  mcube = ( -v - sqrt(wsq))*inv2;
                  ++ncub[18];
               }
               m1 = curoot(mcube);
               if (m1 != doub0)
               {
                  m2 = -uo3/m1;
                  ++ncub[19];
               }
               else
               {
                  m2 = doub0;
                  ++ncub[20];
               }
               v3[0] = m1 + m2 - po3;
            }
            else
            {
/* 
     cubic has three real roots -
*/
               if (uo3 < doub0)
               {
                  muo3 = -uo3;
                  if (muo3 > doub0)
                  {
                     s = sqrt(muo3);
                     ++ncub[21];
                     if (p > doub0)
                     {
                        s = -s;
                        ++ncub[22];
                     }
                  }
		      else
                  {
                     s = doub0;
                     ++ncub[23];
                  }
                  scube = s*muo3;
		      if (scube == doub0)
		      {
                     v3[0] = m1 + m2 - po3;
                     n3 = 1;
                     ++ncub[24];
		      }
                  else
                  {
                     t =  -v/(scube+scube);
                     cosk = acos3(t);
                     v3[0] = (s+s)*cosk - po3;
                     n3 = 1 ;
                     sinsqk = doub1 - cosk*cosk;
		         if (sinsqk >= doub0)
		         {
                        rt3sink = rt3*sqrt(sinsqk);
                        v3[1] = s*(-cosk + rt3sink) - po3;
                        v3[2] = s*(-cosk - rt3sink) - po3;
                        n3 = 3;
                        ++ncub[25];
                     }
                     else ++ncub[26];
                  }
               }
               else
/* 
     cubic has multiple root -  
*/
               {
                  ++ncub[27];
                  v3[0] = curoot(v) - po3;
                  v3[1] = v3[0];
                  v3[2] = v3[0];
                  n3 = 3;
               }
            }
         }
      }
   }
done:
   if (debug < 1)
   {
      printf("cubic %d %g %g %g\n",n3,p,q,r);
      for (j = 0; j < n3; ++j)
         printf("   %d %13g %13g\n",j,v3[j],r+v3[j]*(q+v3[j]*(p+v3[j])));
      printf("v %g,  uo3 %g,  m1 %g,   m2 %g,  po3 %g, cosk %g\n",
         v,uo3,m1,m2,po3,cosk);
      for (j = 0; j < 28; ++j)
      {
         printf("  %d",ncub[j]);
         if ((j%10) == 9) printf("\n");
      }
      printf("\n");
   }
   if (iterate == true) cubnewton(p,q,r,n3,v3);
   return(n3) ;
}; /* cubic */



/****************************************************/
//    improve roots of cubic by Newton-Raphson iteration
//
//    5 Jan 2004  Don Herbison-Evans
//
//    called by cubic.
/****************************************************/

void Quartic::cubnewton(double p,double q,double r,int n3,double v3[4])
{
   int j,k;
   double corr,deriv,err,root;

   if (debug < 2) printf("cubnewtona %d %g\n",
        n3,v3[0]);
   for (j = 0; j < n3; ++j)
   {
      for (k = 0; k < 4; ++k)
      {
         root = v3[j];
         err = ((root + p)*root + q)*root + r;
         deriv = (doub3*root + doub2*p)*root + q;
         if (deriv != doub0) corr = err/deriv; else corr = doub0;
         v3[j] -= corr;
         if (debug < 1) printf("cubnewtonb %d %d %g %g %g %g %g\n",
              j,k,root,err,deriv,corr,v3[j]);
      }
   }
}; /* cubnewton */

/****************************************************/
//   Solve quartic equation using either
//   quadratic, Ferrari's or Neumark's algorithm.
//
//   called by 
//   calls  descartes, ferrari, neumark, yacfraid.
//
//     15 Dec 2003  added yacfraid
//     10 Dec 2003  added descartes with neg coeffs
//     21 Jan 1989  Don Herbison-Evans
/****************************************************/
int Quartic::quartic(double a,double b,double c,double d,double rts[4])
{
   int j,k,nq,nr;
   double roots[4];

   if (Dabs(a) > doubmax)
      nr = yacfraid(a,b,c,d,rts);
   else
   if ((a == doub0) && (c == doub0))
   {
      nq = quadratic(b,d,roots);
      nr = 0;
      for (j = 0; j < nq; ++j)
      {
         if (roots[0] >= doub0)
         {
            rts[0] = sqrt(roots[0]);
            rts[1] = -rts[0];
            nr = 2;
         }
         if (roots[1] >= doub0)
         {
            rts[nr] = sqrt(roots[1]);
            rts[nr+1] = -rts[nr];
            nr += 2;
         }
      }
   }
   else
   {
      k = 0;
      if (a < doub0) k += 2; 
      if (b < doub0) k += 1; 
      if (c < doub0) k += 8; 
      if (d < doub0) k += 4;
      switch (k)
      {
              case 0 : nr = neumark(a,b,c,d,rts); break;
              case 1 : nr = neumark(a,b,c,d,rts); break;
              case 2 : nr = neumark(a,b,c,d,rts); break;
              case 3 : nr = ferrari(a,b,c,d,rts); break;
              case 4 : nr = neumark(a,b,c,d,rts); break;
              case 5 : nr = descartes(a,b,c,d,rts); break;
              case 6 : nr = neumark(a,b,c,d,rts); break;
              case 7 : nr = neumark(a,b,c,d,rts); break;
              case 8 : nr = neumark(a,b,c,d,rts); break;
              case 9 : nr = ferrari(a,b,c,d,rts); break;
              case 10 : nr = neumark(a,b,c,d,rts); break;
              case 11 : nr = neumark(a,b,c,d,rts); break;
              case 12 : nr = neumark(a,b,c,d,rts); break;
              case 13 : nr = neumark(a,b,c,d,rts); break;
              case 14 : nr = neumark(a,b,c,d,rts); break;
              case 15 : nr = descartes(-a,b,-c,d,rts); break;
      }
      if (k == 15)
        for (j = 0; j < nr; ++j)
           rts[j] = -rts[j];
   }
   return(nr);
}; /* quartic */


//*****************************************/
//   Solve quartic equation using
//   Descartes-Euler-Cardano algorithm
//
//   called by quartic, compare, quartictest.
//
//   Strong, T. "Elemementary and Higher Algebra"
//      Pratt and Oakley, p. 469 (1859)
//
//   16 Jul 1981  Don Herbison-Evans
//*****************************************/
int Quartic::descartes(double a,double b,double c,double d,double rts[4])
{
   int j;
   double v1[4],v2[4],v3[4];
   double k,y;
   double p,q,r;
   double e0,e1,e2;
   double g,h;
   double asq;
   double ainv4;
   double e1invk;

   asq = a*a;
   e2 = b - asq*d3o8;
   e1 = c + a*(asq*inv8 - b*inv2);
   e0 = d + asq*(b*inv16 - asq*d3o256) - a*c*inv4;

   p = doub2*e2;
   q = e2*e2 - doub4*e0;
   r = -e1*e1;

   n3 = cubic(p,q,r,v3);
   for (j3 = 0; j3 < n3; ++j3)
   {
      y = v3[j3];
      if (y <= doub0)
      { 
         n4[j3] = 0;
         ++ndes[0];
      } /* y<0 */
      else
      {
         k = sqrt(y);
         ainv4 = a*inv4;
         e1invk = e1/k;
         g = (y + e2 + e1invk)*inv2;
         h = (y + e2 - e1invk)*inv2 ;
         n1 = quadratic(-k, g, v1);
         n2 = quadratic( k, h, v2);
         qrts[0][j3] = v1[0] - ainv4;
         qrts[1][j3] = v1[1] - ainv4;
         qrts[n1][j3] = v2[0] - ainv4;
         qrts[n1+1][j3] = v2[1] - ainv4;
         n4[j3]= n1 + n2;  
         ++ndes[1];
      } /* y>=0 */
donej3:
      for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
      worst3[j3] = errors(a,b,c,d,rts,rterd,n4[j3]);
   } /* j3 loop */
done:
   j3 = 0;
   if (n3 != 1)
   {
      if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
         ++ndes[NCASES-2];
      else 
         ++ndes[NCASES-1];
      if ((n4[1] > n4[j3]) || 
         ((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
      if ((n4[2] > n4[j3]) ||
         ((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
   }
   for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
   if (debug < 1)
       printf("descartes chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
   ++ndes[30+n4[j3]];
   ++ndes[35+j3];
   return(n4[j3]);
}; /* descartes */

//****************************************************/
//     solve the quartic equation - 
//
//     x**4 + a*x**3 + b*x**2 + c*x + d = 0 
//
//     calls     cubic, quadratic.
//
//     input - 
//     a,b,c,e - coeffs of equation. 
//
//     output - 
//     n4 - number of real roots. 
//     rts - array of root values. 
//
//     method :  Ferrari - Lagrange
//     Theory of Equations, H.W. Turnbull p. 140 (1947)
//
//     16 Jul 1981 Don Herbison-Evans
//
//     calls  cubic, quadratic 
//****************************************************/
int Quartic::ferrari(double a,double b,double c,double d,double rts[4])
{
   int j;
   double asqinv4;
   double ainv2;
   double d4;
   double yinv2;
   double v1[4],v2[4],v3[4];
   double p,q,r;
   double y;
   double e,f,esq,fsq,ef;
   double g,gg,h,hh;

   ainv2 = a*inv2;
   asqinv4 = ainv2*ainv2;
   d4 = d*doub4 ;

   p = b;
   q = a*c-d4;
   r = (asqinv4 - b)*d4 + c*c;
   n3 = cubic(p,q,r,v3);
   for (j3 = 0; j3 < n3; ++j3)
   {
      y = v3[j3];
      yinv2 = y*inv2;
      esq = asqinv4 - b - y;
      fsq = yinv2*yinv2 - d;
      if ((esq < doub0) && (fsq < doub0))
      {
         n4[j3] = 0;
         ++nfer[0];
      }
      else
      {
         ef = -(inv4*a*y + inv2*c);
         if ( ((a > doub0)&&(y > doub0)&&(c > doub0))
           || ((a > doub0)&&(y < doub0)&&(c < doub0))
           || ((a < doub0)&&(y > doub0)&&(c < doub0))
           || ((a < doub0)&&(y < doub0)&&(c > doub0))
           ||  (a == doub0)||(y == doub0)||(c == doub0))
/* use ef - */
         {
               if ((b < doub0)&&(y < doub0))
               {
                  e = sqrt(esq);
                  f = ef/e;
                  ++nfer[1];
               }
               else if (d < doub0)
               {
                  f = sqrt(fsq);
                  e = ef/f;
                  ++nfer[2];
               }
               else
               {
                  if (esq > doub0)
                  {
                     e = sqrt(esq);
                     ++nfer[3];
                  }
                  else
                  {
                     e = doub0;
                     ++nfer[4];
                  }
                  if (fsq > doub0)
                  {
                     f = sqrt(fsq);
                     ++nfer[5];
                  }
                  else
                  {
                     f = doub0;
                     ++nfer[6];
                  }
                  if (ef < doub0)
                  {
                     f = -f;
                     ++nfer[7];
                  }
               }
         }
         else
/* use esq and fsq - */
         {
               if (esq > doub0)
               {
                  e = sqrt(esq);
                  ++nfer[8];
               }
               else
               {
                  e = doub0;
                  ++nfer[9];
               }
               if (fsq > doub0)
               {
                  f = sqrt(fsq);
                  ++nfer[10];
               }
               else
               {
                  f = doub0;
                  ++nfer[11];
               }
               if (ef < doub0)
               {
                  f = -f;
                  ++nfer[12];
               }
            }
/* note that e >= doub0 */
            g = ainv2 - e;
            gg = ainv2 + e;
            if ( ((b > doub0)&&(y > doub0))
               || ((b < doub0)&&(y < doub0)) )
            {
               if ( ((a > doub0) && (e > doub0))
	                || ((a < doub0) && (e < doub0)) )
               {
                  g = (b + y)/gg;
                  ++nfer[13];
               }
               else
	         if (    ((a > doub0) && (e < doub0))
	              || ((a < doub0) && (e > doub0)) )
               {
                  gg = (b + y)/g;
                  ++nfer[14];
               }
               else
                  ++nfer[15];
            }
            hh = -yinv2 + f;
            h = -yinv2 - f;
            if ( ((f > doub0)&&(y < doub0))
              || ((f < doub0)&&(y > doub0)) )
            {
               h = d/hh;
               ++nfer[16];
            }
            else
            if ( ((f < doub0)&&(y < doub0))
               || ((f > doub0)&&(y > doub0)) )
            {
               hh = d/h;
               ++nfer[17];
            }
            else
               ++nfer[18];

            n1 = quadratic(gg,hh,v1);
            n2 = quadratic(g,h,v2);
            n4[j3] = n1+n2;
            qrts[0][j3] = v1[0];
            qrts[1][j3] = v1[1];
            qrts[n1+0][j3] = v2[0];
            qrts[n1+1][j3] = v2[1];
      } 
donej3:
      for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
      worst3[j3] = errors(a,b,c,d,rts,rterf,n4[j3]);
   } /* j3 loop */
done:
   j3 = 0;
   if (n3 != 1)
   {
      if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
         ++nfer[NCASES-2];
      else 
         ++nfer[NCASES-1];
      if ((n4[1] > n4[j3]) || 
         ((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
      if ((n4[2] > n4[j3]) ||
         ((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
   }
   for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
   if (debug < 1)
       printf("ferrari chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
   ++nfer[30+n4[j3]];
   ++nfer[35+j3];
   return(n4[j3]);
}; /* ferrari */


//*****************************************/
//     solve the quartic equation - 
//
//     x**4 + a*x**3 + b*x**2 + c*x + d = 0 
//
//    calls     cubic, quadratic.
//
//    input parameters - 
//     a,b,c,e - coeffs of equation. 
//
//     output parameters - 
//     n4 - number of real roots. 
//     rts - array of root values. 
//
//     method -  S. Neumark 
//        "Solution of Cubic and Quartic Equations" - Pergamon 1965 
//
//      1 Dec 1985   translated to C with help of Shawn Neely
//     16 Jul 1981   Don Herbison-Evans
//*****************************************/
int Quartic::neumark(double a,double b,double c,double d,double rts[4])
{
   int j;
   double y,g,gg,h,hh,gdis,gdisrt,hdis,hdisrt,g1,g2,h1,h2;
   double bmy,gerr,herr,y4,bmysq;
   double v1[4],v2[4],v3[4];
   double asq;
   double d4;
   double p,q,r;
   double hmax,gmax;

   if (d == doub0)
   {
      n3 = 0;
      n4[0] = cubic(a,b,c,rts);
      for (j = 0; j < n4[0]; ++j)
         qrts[j][0] = rts[j];
      qrts[n4[0]++][0] = doub0;
      goto done;
   }
   asq = a*a;
   d4 = d*doub4;
   p =  -b*doub2;
   q = b*b + a*c - d4;
   r = (c - a*b)*c + asq*d;
   if (debug < 3)
      printf("neumarka %g %g %g %g,  %g %g %g\n",
      a,b,c,d,p,q,r);
   n3 = cubic(p,q,r,v3);
   for (j3 = 0; j3 < n3; ++j3)
   {
      y = v3[j3];
      bmy = b - y;
      y4 = y*doub4;
      bmysq = bmy*bmy;
      gdis = asq - y4;
      hdis = bmysq - d4;
      if (debug < 3)
         printf("neumarkb %g %g\n",
            gdis,hdis);
      if ((gdis < doub0) || (hdis < doub0))
      {
         n4[j3] = 0;
         ++nneu[0];
      }
      else   
      {
         g1 = a*inv2;
         h1 = bmy*inv2;
         gerr = asq + y4;
         herr = hdis;
         if (d > doub0)
         {
            herr = bmysq + d4;
            ++nneu[1];
         }
         if ((y < doub0) || (herr*gdis > gerr*hdis))
         {
            gdisrt = sqrt(gdis);
            g2 = gdisrt*inv2;
            if (gdisrt != doub0)
            {
               h2 = (a*h1 - c)/gdisrt;
               ++nneu[2];
            }
            else
            {
               h2 = doub0;
               ++nneu[3];
            }   
         }
         else
         {
            hdisrt = sqrt(hdis);
            h2 = hdisrt*inv2;
            if (hdisrt != doub0)
            {
               g2 = (a*h1 - c)/hdisrt;
               ++nneu[4];
            }
            else
            {
               g2 = doub0;
               ++nneu[5];
            }
         }
/* 
     note that in the following, the tests ensure non-zero 
     denominators -  
*/
         h = h1 - h2;
         hh = h1 + h2;
         hmax = hh;
         if (hmax < doub0)
         {
            hmax =  -hmax;
            ++nneu[6];
         }
         if (hmax < h)
         {
            hmax = h;
            ++nneu[7];
         }   
         if (hmax <  -h)
         {
            hmax =  -h;
            ++nneu[8];
         }
         if ((h1 > doub0)&&(h2 > doub0))
         {
            h = d/hh;
            ++nneu[9];
         }
         if ((h1 < doub0)&&(h2 < doub0))
         {
            h = d/hh;
            ++nneu[10];
         }
         if ((h1 > doub0)&&(h2 < doub0))
         {
            hh = d/h;
            ++nneu[11];
         }
         if ((h1 < doub0)&&(h2 > doub0))
         {
            hh = d/h;
            ++nneu[12];
         }
         if (h > hmax)   
         {
            h = hmax;
            ++nneu[13];
         }
         if (h <  -hmax)
         {
            h =  -hmax;
            ++nneu[14];
         }
         if (hh > hmax)
         {
            hh = hmax;
            ++nneu[15];
         }
         if (hh < -hmax)
         {
            hh =  -hmax;
            ++nneu[16];
         }

         g = g1 - g2;
         gg = g1 + g2;
         gmax = gg;
         if (gmax < doub0)
         {
            gmax = -gmax;
            ++nneu[17];
         }
         if (gmax < g)
         {
            gmax = g;
            ++nneu[18];
         }
         if (gmax < -g)
         {
            gmax = -g;
            ++nneu[19];
         }
         if ((g1 > doub0)&&(g2 > doub0))
         {
            g = y/gg;
            ++nneu[20];
         }
         if ((g1 < doub0)&&(g2 < doub0))
         {
            g = y/gg;
            ++nneu[21];
         }
         if ((g1 > doub0)&&(g2 < doub0))
         {
            gg = y/g;
            ++nneu[22];
         }
         if ((g1 < doub0)&&(g2 > doub0))
         {
            gg = y/g;
            ++nneu[23];
         }
         if (g > gmax)
         {
            g = gmax;
            ++nneu[24];
         }
         if (g <  -gmax)
         {
            g = -gmax;
            ++nneu[25];
         }
         if (gg > gmax)
         {   
            gg = gmax;
            ++nneu[26];
         }
         if (gg <  -gmax)
         {
            gg = -gmax;
            ++nneu[27];
         }
 
         n1 = quadratic(gg,hh,v1);
         n2 = quadratic(g,h,v2);
         n4[j3] = n1 + n2;
         qrts[0][j3] = v1[0];
         qrts[1][j3] = v1[1];
         qrts[n1+0][j3] = v2[0];
         qrts[n1+1][j3] = v2[1];
      } 
donej3:
      for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
      worst3[j3] = errors(a,b,c,d,rts,rtern,n4[j3]);
   } /* j3 loop */
done:
   j3 = 0;
   if (n3 > 1)
   {
      if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
         ++nneu[NCASES-2];
      else 
         ++nneu[NCASES-1];
      if ((n4[1] > n4[j3]) || 
         ((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
      if ((n4[2] > n4[j3]) ||
         ((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
   }
   for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
   if (debug < 1)
       printf("neumark chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
   ++nneu[30+n4[j3]];
   ++nneu[35+j3];
   return(n4[j3]);
}; /* neumark */



//****************************************************/ 
//     solve the quartic equation - 
//     x**4 + a*x**3 + b*x**2 + c*x + d = 0 
//
//     calls     cubic, quadratic.
//
//     input parameters - 
//     a,b,c,e - coeffs of equation. 
//
//     output parameters - 
//     n4 - number of real roots. 
//     rts - array of root values. 
//
//     method - 
//        K.S. Brown 
//        Reducing Quartics to Cubics,
//        http://www.seanet.com/~ksbrown/kmath296.htm (1967)
// 
//        Michael Daoud Yacoub & Gustavo Fraidenraich
//        "A new simple solution of the general quartic equation"
//        Revised 16 Feb 2004
//
//     14 Nov 2003 Don Herbison-Evans
//*****************************************/
int Quartic::yacfraid(double a,double b,double c,double d,double rts[4])
{
   int j;
   double y;
   double v3[4];
   double asq,acu;
   double b4;
   double det0,det1,det2,det3;
   double det0rt,det1rt,det2rt,det3rt;
   double e,f,g,h,k;
   double fsq,gsq,hsq,invk;
   double P,Q,R,U;

   if (d == doub0)
   {
      n3 = 0;
      n4[0] = cubic(a,b,c,rts);
      for (j = 0; j < n4[0]; ++j)
         qrts[j][0] = rts[j];
      qrts[n4[0]++][0] = doub0;
      goto done;
   }
   asq = a*a;
   acu = a*asq;
   b4 = b*doub4;
   n3 = 0;

   P = asq*b - b4*b + doub2*a*c + doub16*d ;
   Q = asq*c - b4*c + doub8*a*d;
   R = asq*d - c*c ;
   U = acu - b4*a + doub8*c;
   n4[0] = 0; 
   if (U == doub0)
   {
      if (P == doub0)
      {
         det0 = doub3*asq - doub8*b;
         if (det0 < doub0)
         {
            ++nyac[0];
            goto done;
         }
         det0rt = sqrt(det0);
         qrts[0][0] = (-a + det0rt)*inv4;
         qrts[1][0] = qrts[0][0];
         qrts[2][0] = (-a - det0rt)*inv4;
         qrts[3][0] = qrts[2][0];
         ++nyac[1];
         n4[0] = 4;
         goto done;
      } /* P=0 */
      else
      {
         det1 = asq*asq - doub8*asq*b + doub16*b*b - doub64*d;
         if (det1 < doub0)
         {
            ++nyac[2];
            goto done;;
         }
         n4[0] = 0;
         det1rt =  sqrt(det1);
         det2 = doub3*asq - doub8*b + doub2*det1rt;
         if (det2 >= doub0)
         {
            det2rt = sqrt(det2);
            qrts[0][0] = (-a + det2rt)*inv4;
            qrts[1][0] = (-a - det2rt)*inv4;
            n4[0] = 2;
            ++nyac[3];
         }
         det3 = doub3*asq - doub8*b - doub2*det1rt;
         if (det3 >= doub0)
         {
            det3rt = sqrt(det3);
            qrts[n4[0]++][0] = (-a + det3rt)*inv4;
            qrts[n4[0]++][0] = (-a - det3rt)*inv4;
            ++nyac[5];
         }
         if (n4[0] == 0) ++nyac[6];
         goto done;
      } /* P<>0 */
   }

   n3 = cubic(P/U,Q/U,R/U,v3);
   for (j3 = 0; j3 < n3; ++j3)
   {
      y = v3[j3];
      j = 0;
      k = a + doub4*y;
      if (k == doub0)
      {
         ++nyac[9];
         goto donej3;
      }
      invk = doub1/k;
      e = (acu - doub4*c - doub2*a*b + (doub6*asq - doub16*b)*y)*invk;
      fsq = (acu + doub8*c - doub4*a*b)*invk;
      if (fsq < doub0)
      {
         ++nyac[10];
         goto donej3;
      }
      f = sqrt(fsq);
      gsq = doub2*(e + f*k);
      hsq = doub2*(e - f*k);
      if (gsq >= doub0)
      {
         ++nyac[11];
         g = sqrt(gsq);
         qrts[j++][j3] = (-a - f - g)*inv4;
         qrts[j++][j3] = (-a - f + g)*inv4;
      }
      if (hsq >= doub0)
      {
         ++nyac[12];
         h = sqrt(hsq);
         qrts[j++][j3] = (-a + f - h)*inv4;
         qrts[j++][j3] = (-a + f + h)*inv4;
      }
      if (debug < 1)
      {
         printf("j3 %d y %g k %g fsq %g gsq %g hsq %g\n",
            j3,y,k,fsq,gsq,hsq);
         printf("e %g f %g g %g h %g\n",e,f,g,h);
      }
donej3:
      n4[j3] = j;
      for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
      worst3[j3] = errors(a,b,c,d,rts,rtery,n4[j3]);
   } /* j3 loop */
done:
   j3 = 0;
   if (n3 > 1)
   {
      if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
         ++nyac[NCASES-2];
      else 
      {
         ++nyac[NCASES-1];
         if ((n4[0] != n4[1]) && (n4[0] != n4[2]) && (n4[1] != n4[2]))
            printf("yace %d %d %d %g %g %g %g\n",
               n4[0],n4[1],n4[2],a,b,c,d);
      }
      if ((n4[1] > n4[j3]) || 
         ((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
      if ((n4[2] > n4[j3]) ||
         ((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
   }
   for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
   ++nyac[30+n4[j3]];
   ++nyac[35+j3];
   if (debug < 1)
       printf("yacfraid chose cubic %d %g  %g\n\n",j3,v3[j3],worst3[j3]);
   return(n4[j3]);
}; /* yacfraid */


//*****************************************/
//   Solve quartic equation using
//   Christianson's algorithm
//   calls errors, quadratic, cubic.
//
//   Bruce Christianson, Solving Quartics Using Palindromes,
//      Mathematical Gazette, Vol. 75, pp. 327-328 (1991)
//
//   14 Jan 2004  Don Herbison-Evans
//*****************************************/

int Quartic::chris(double a,double b,double c,double d,double rts[4])
{
   int j;
   int n0,n1,n2;
   double v0[4],v1[4],v2[4],v3[4];
   double y,ysq,ycu;
   double p,q,r;
   double asq,acu,aqu,ao4;
   double e0,e1,e2;
   double k,ksq,kcu,kqu,kquinv,k1,k2;
   double g,h,g1,Z1,Z2;

   if (d == doub0)
   {
      n3 = 0;
      n4[0] = cubic(a,b,c,rts);
      for (j = 0; j < n4[0]; ++j)
         qrts[j][0] = rts[j];
      qrts[n4[0]++][0] = doub0;
      goto done;
   }
   asq = a*a;
   acu = asq*a;
   aqu = acu*a;
   ao4 = a*inv4;
   e2 = b - d3o8*asq;
   e1 = c - inv2*b*a + inv8*acu;
   e0 = d - inv4*c*a + inv16*b*asq - d3o256*acu*a;
   if (debug < 1)
      printf("chrisa e0 %g   e1 %g   e2 %g\n",e0,e1,e2);
   if (e1 == doub0)
   { 
      n3 = 0;
      n4[0] = 0;
      ++nchr[0];
      n3 = 3;
      v3[2] = doub0;
      v3[1] = doub0;
      v3[0] = -inv8*(doub16*e0 - doub4*e2*e2);
   }
   else
   {
      p = (inv2*b*asq - inv2*b*b - inv2*c*a + doub2*d - doub3*inv32*aqu)/e1;
      q = doub3*inv16*asq - inv2*b;
      r = inv16*b*a - inv8*c - inv64*acu;
      if (debug < 1) printf("chrisb %g %g %g\n",p,q,r);
      ++nchr[1];
      n3 = cubic(p,q,r,v3);
   }
   for (j3 = 0; j3 < n3; ++j3)
   {
      y = v3[j3];
      n4[j3] = 0;
      ysq = y*y;
      ycu = y*ysq;
      if ((y == doub0) ||
         (( y < doub0) && (e1 <= doub0) && (e0 >= doub0) && (e2 >= doub0)))
      {
         kqu = y*ycu + e2*ysq + e1*y + e0;
         if (kqu <= doub0)
         {
             ksq = doub0;
             ++nchr[2];
             goto donej3;
         }
         ksq = sqrt(kqu);
      }
      else
      {
         ksq = ysq + inv2*e2 + inv4*e1/y;
         if (ksq <= doub0)
         {
             kqu = doub0;
             ++nchr[3];
             goto donej3;
         }
         kqu = ksq*ksq;
      }
      k = sqrt(ksq);
      kcu = k*ksq;
      kquinv = doub1/kqu;
      g = doub4*y*kcu;
      h = (doub6*ysq + e2)*ksq;
      if (debug < 1)
      {
         k2 = sqrt(sqrt(Dabs(y*ycu + e2*ysq + e1*y + e0)));
         if (y == doub0)
             k1 = doub0;
         else
             k1 = sqrt(Dabs(ysq + inv2*e2 + inv4*e1/y));
         g1 = (doub4*ycu + doub2* y*e2 + e1)*k;
         printf("chrisc %g %g %g %g %g %g\n",k,k1,y,g,g1,h);
      }
      n0 = quadratic(g*kquinv,h*kquinv-doub2,v0);
      if (n0 < 1)
      {
         ++nchr[4]; 
         goto donej3;
      }
      Z1 = v0[0];
      Z2 = v0[1];
      n1 = quadratic(-Z1,doub1,v1);
      if (n1 > 0)
      {
          ++nchr[5];
          n4[j3] = n1;
          qrts[0][j3] = y + k*v1[0] - ao4;
          qrts[1][j3] = y + k*v1[1] - ao4;
      }
      n2 = quadratic(-Z2,doub1,v2);
      if (n2 > 0)
      {
          ++nchr[6];
          n4[j3] += n2;
          qrts[n1][j3] = y + k*v2[0] - ao4;
          qrts[n1+1][j3] = y + k*v2[1] - ao4;
      }
donej3:
      if (debug < 1)
         printf("chrisd %g %g\n",ksq,kqu);
      for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
      worst3[j3] = errors(a,b,c,d,rts,rterc,n4[j3]);
   } /* j3 loop */
done:
   j3 = 0;
   if (n3 > 1)
   {
      if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
         ++nchr[NCASES-2];
      else 
      {
         ++nchr[NCASES-1];
         if (debug < 1)
            printf("chrise %d %d %d %g %g %g %g\n",
               n4[0],n4[1],n4[2],a,b,c,d);
      }
      if ((n4[1] > n4[j3]) || 
         ((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
      if ((n4[2] > n4[j3]) ||
         ((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
   }
   for (j = 0; j < n4[j3]; ++j)
         rts[j] = qrts[j][j3];
   if (debug < 1)
       printf("chris chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
   ++nchr[30+n4[j3]];
   ++nchr[35+j3];
   return(n4[j3]);
}; /* chris */