/*   quartic.c
     version 54

     a test of relative accuracy of various quartic routines.

     use:
       quartic [-a] [-c n] [-q n] [-d n]
       (no parameters)  : a loop through 10,000 prechosen quartic coefficients
       -a               : improve cubic roots (default: no iteration)
       -c  n            : keep reading n cubic roots from standard input
                          (if n=0 read coefficients), terminate with 'q'.
       -d  n            : set debug value to 'n' (default: 5)
       -q  n            : keep reading n quartic roots from standard input
                          (if n=0 read coefficients), terminate with 'q'.

     debug <= 1 : print cases used

     method-
      http://linus.socs.uts.edu.au/~don/pubs/solving.html

      4 Apr 2004 bringing yacfraid notation into line with solving.html
      4 Apr 2004 fixed bug in final root calculation in yacfraid
      8 Mar 2004 corrected modified yacfraid algorithm
      8 Mar 2004 modified yacfraid algorithm
     31 Jan 2004 printing results when number of roots vary 
     30 Jan 2004 printing error of chosen cubic if debug<1
     27 Jan 2004 seeking worst coefficients for each algorithm
     27 Jan 2004 choosing best route for k in chris
     26 Jan 2004 conforming variables to solving.html
     26 Jan 2004 fixed bug in yacfraid for multiplicity 3
     25 Jan 2004 speeding up chris
     24 Jan 2004 fixed chris error (A <-> B)
     23 Jan 2004 use debug<1 for diagnostics
     21 Jan 2004 solve cubic in neumark, yacfraid and chris if d=0
     20 Jan 2004 3 roots to cubic in chris if e1=0
     19 Jan 2004 debugging Christianson routine
     14 Jan 2004 adding Christianson, cut determinant in quadratic call
      5 Jan 2004 improving cubic roots by Newton-Raphson iteration
     23 Dec 2003 seeking snag in yacfraid
     21 Dec 2003 putting diagnostic printout in yacfraid
     18 Dec 2003 allowing input of equation coefficients
     18 Dec 2003 recording cubic root giving least worst error
     17 Dec 2003 recording cubic root giving the most quartic roots
     16 Dec 2003 using the cubic root giving the most quartic roots
     16 Dec 2003 trying consistency of 3 cubic roots where available
     15 Dec 2003 trying all 3 cubic roots where available
     13 Dec 2003 trying to fix Neumark
     13 Dec 2003 cleaning up diagnostic format
     12 Dec 2003 initialising n,m,po3 in cubic
     12 Dec 2003 allow cubic to return 3 zero roots if p=q=r=0
     10 Dec 2003 added setargs
      2 Dec 2003 finding worst cases
      2 Dec 2003 negating j if p>0 in cubic
      1 Dec 2003 changing v in cubic from (sinsqk > doub0) to (sinsqk >= doub0)
      1 Dec 2003 test changing v in cubic from po3sq+po3sq to doub2*po3sq
     30 Nov 2003 counting cases in all solving routines
     29 Nov 2003 testing wsq >= doub0
     29 Nov 2003 mult by doub2 for v in cubic
     29 Nov 2003 cut po3cu from cubic
     29 Nov 2003 better quadratic
     29 Nov 2003 count agreements
     17 Nov 2003 count special cases
     16 Nov 2003 option of loop or read coefficients from input
     15 Nov 2003 fixed cubic() bug
     11 Nov 2003 added Brown and Yacoub & Fraidenraich's algorithms
     21 Jan 1989 quartic selecting Descartes, Ferrari, or Neumark algorithms 
     16 Jul 1981 Don Herbison-Evans

"Solving Quartics and Cubics for Graphics", D. Herbison-Evans,
Graphics Gems V (ed.: A. Paeth) Academic Press (Chesnut Hill), pp. 3-15 (1995). 

"Solving Quartics and Cubics for Graphics", D. Herbison-Evans,
Research Report CS-86-56, Department of Computer Science, University of Waterloo (1986)
 
"Caterpillars and the Inaccurate Solution of Cubic and Quartic Equations", 
D. Herbison-Evans, Australian Computer Science Communications, 
Vol. 5, No. 1, pp. 80-90 (1983)

   subroutines:
      setcns      - set constants
      setargs     - determine which equations to test
      looptest    - test 10000 coefficient combinations
      cases       - print statistics
      docoeffs    - read coefficients of a particular equation to solve
      cubictest   - test cubic solutions
      quartictest - test quartic solutions
      compare     - check accuracy of results
      errors      - find errors in a set of roots
      acos3       - find arcos(x/3)
      curoot      - find cube root
      quadratic   - solve a quadratic
      cubic       - solve a cubic
      cubnewton   - improve cubic roots by iteration
      quartic     - solve a quartic
      descartes   - use Descartes' algorithm to solve a quartic  
      ferrari     - use Ferrari's algorithm to solve a quartic 
      neumark     - use Neumark's algorithm to solve a quartic 
      yacfraid    - use Yacoub & Fraidenraich's algorithm to solve a quartic 
      chris       - use Christianson's algorithm to solve a quartic
*/

#include <stdio.h>
#define TRUE 1
#define FALSE 0
#define NCASES 40

double a,b,c,d;           /* quartic coefficients */
double cc[4],cd[4],cf[4],cn[4],cy[4];
double d3o8,d3o256;       /* double precision constants */
double doub0;
double doub1,doub2;
double doub3,doub4;
double doub5,doub6;
double doub8,doub9,doub12;
double doub16,doub24;
double doub27,doub64;
double doubmax;           /* approx square root of max double number */
double doubmin;           /* smallest double number */
double doubtol;           /* tolerance of double numbers */
double errq;
double inv2,inv3,inv4;
double inv8,inv16,inv32,inv64,inv128;
double maxc,maxd,maxf,maxn,maxy;
double p,q,r;             /* cubic coefficients */
double qrts[4][3];        /* quartic roots for each cubic root */
double rt3;               /* square root of 3 */
double rtsc[4],rtsd[4],rtsf[4],rtsn[4],rtsq[4],rtsy[4]; /* roots given by each algorithm */
double rterc[4],rterd[4],rterf[4],rtern[4],rterq[4],rtery[4]; /* errors of roots */
double worstc,worstd,worstf,worstn,worsty; /* worst error for each algotithm */
double worst3[3];         /* worst error for a given cubic root */
double x0,x1,x2,x3;

int agr,dis,zer;
int agreecd,agreecf,agreecn,agreecy;
int agreedf,agreedn,agreedy;
int agreefn,agreefy,agreeny;
int debug;                /* <1 for lots of diagnostics, >5 for none */   
int docubic,doquartic;
int iterate;
int j3;
int n1,n2,n3,n4[3];
int n,nc,nd,nf,nn,ny;
int nrtsc,nrtsd,nrtsf,nrtsn,nrtsq,nrtsy;
int nt3,nt2,nt2sq,nden,nm4,ndet1,ndet2,ndet3,ndet4,ndet5,nn;
int nqud[NCASES];
int ncub[NCASES];
int nchr[NCASES];
int ndes[NCASES];
int nfer[NCASES];
int nneu[NCASES];
int nyac[NCASES];
int nzc,nzd,nzf,nzn,nzy;
int zerodf,zerodn,zerody,zerofn,zerofy,zerony;
int zerodfn,zerodfy,zerodny,zerofny;
int tot;

double exp(),log(),sqrt(),cos(),acos(),fabs(),errors();
int cubic() ;
int descartes(),ferrari(),neumark(),yacfraid();

/***********************************/

main(argc,argv)
int argc;
char *argv[];
{
   double fabs();

   printf("quartic version 54\n");
   debug = 5;
   setcns();
   setargs(argc,argv);

more: if (argc > 1)
   {
      docoeffs();
      if ((fabs(x0)+fabs(x1)+fabs(x2)+fabs(x3)) != doub0)
         goto more;
      exit(0);
   }
   else
   {
      looptest();
      exit(0);
   }
} /* main */
/****************************************/

setcns()
/* 
     set up constants.

     called by main, quartictest.
*/
{
      int j;

      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 (j = 1; doub1+doubtol > doub1; ++ j)
      {
          doubtol *= inv2;
      }
      doubtol = sqrt(doubtol);

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

} /* setcns */
/***********************************/

setargs(argc,argv)
int argc;
char *argv[];
/*
   determine whether to do looptest or just individual equation

   10 Dec 2003 Don Herbison-Evans

   called by main
*/
{
   int j;

   docubic = FALSE;
   doquartic = FALSE;
   iterate = FALSE;
   for (j = 1; j < argc; ++j)
   {
      if (argv[j][0] == '-')
      {
         if (argv[j][1] == 'a')
         {
            iterate = TRUE;
            printf("iterate cubics\n");
         }
         else
         if (argv[j][1] == 'c')
         {
            docubic = TRUE;
            sscanf(argv[++j],"%d",&n);
            printf("do cubic %d\n",n);
         }
         else
         if (argv[j][1] == 'q') 
         {
            doquartic = TRUE;
            sscanf(argv[++j],"%d",&n);
            printf("do quartic %d\n",n);
         }
         else
         if (argv[j][1] == 'd')
         {
            sscanf(argv[++j],"%d",&debug);
            printf("debug %d\n",debug);
         }
      }
   }
} /* setargs */
/********************************************************/

looptest()
/*
   loop through a range of coefficient values 
   comparing the algorithms
  
   called by main.
   calls compare.
*/
{
   double ten4,ten8;
   double v[10];
   int j,k,l,m;

   printf("quartic loop test\n");
   nc = 0;
   nd = 0;
   nf = 0;
   nn = 0;
   ny = 0;
   maxc = doub0;
   maxd = doub0;
   maxf = doub0;
   maxn = doub0;
   maxy = doub0;
   tot = 0;
   zer = 0;
   agr = 0;
   dis = 0;
   agreecd = 0;
   agreecf = 0;
   agreecn = 0;
   agreecy = 0;
   agreedf = 0;
   agreedn = 0;
   agreedy = 0;
   agreefn = 0;
   agreefy = 0;
   agreeny = 0;
   for (j = 0; j < NCASES; ++j)
   {
      nqud[j] = 0;
      ncub[j] = 0;
      ndes[j] = 0;
      nfer[j] = 0;
      nneu[j] = 0;
      nyac[j] = 0;
   }
   ten4 = (double)10000;
   ten8 = ten4*ten4;
   j = 0;
   for (a = ten8; a > inv2/ten8; a /= ten4)
   {
      v[j] = a;
      ++j;
   }
   for (j = 0; j < 5; ++j)
      v[j+5] = -v[j];
   for (j = 0; j < 10; ++j)
   {
      a = v[j];
      for (k = 0; k < 10; ++k)
      {
         b = v[k];
         for (l = 0; l < 10; ++l)
         {
            c = v[l];
            for (m = 0; m < 10; ++m)
            {
               d = v[m];
               compare();
            } /* for m */
         } /* for l */
      } /* for k */
   } /* for j */
   printf("total cases: %d\n",
          tot);
   printf("number of real roots: five agree %d, disagree %d\n",
          agr,dis);
   printf("two agree: cd %d, cf %d, cn %d, cy %d\n",
          agreecd,agreecf,agreecn,agreecy);
   printf("           df %d, dn %d, dy %d, fn %d, fy %d, ny %d\n",
          agreedf,agreedn,agreedy,agreefn,agreefy,agreeny);
   printf("5 agree on no real roots: %d\n",
          zer);
   printf("no real roots: chris %d, desc %d, ferr %d, neum %d, yacf %d\n",
          nzc,nzd,nzf,nzn,nzy);
   printf("Christianson : best %d, worst %g %g %g %g, error %g\n",
          nc,cc[0],cc[1],cc[2],cc[3],maxc);
   printf("Descartes : best %d, worst %g %g %g %g, error %g\n",
          nd,cd[0],cd[1],cd[2],cd[3],maxd);
   printf("Ferrari : best %d, worst %g %g %g %g, error %g\n",
          nf,cf[0],cf[1],cf[2],cf[3],maxf);
   printf("Neumark : best %d, worst %g %g %g %g, error %g\n",
          nn,cn[0],cn[1],cn[2],cn[3],maxn);
   printf("Yacoub  : best %d, worst %g %g %g %g, error %g\n",
          ny,cy[0],cy[1],cy[2],cy[3],maxy);
   cases();
} /* looptest */
/********************************************/

cases()
/*
   print statistics

    2 Dec 2003 Don Herbison-Evans

   called by looptest, quartictest.
*/
{
   int j;

   printf("quadratic cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,nqud[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

   printf("cubic cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,ncub[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

   printf("descartes cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,ndes[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

   printf("ferrari cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,nfer[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

   printf("neumark cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,nneu[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

   printf("yacfraid cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,nyac[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

   printf("chris cases:\n");
   for (j = 0; j < NCASES; ++j)
   {
      printf("  %2d %d,",j,nchr[j]);
      if (j%5 == 4) printf("\n");
   }
   printf("\n");

} /* cases */
/************************************/

docoeffs()
/*
    read a set of coefficients from standard input
    and compare results of the different algorithms.

    16 Nov 2003 Don Herbison-Evans

    called by main.
    calls cubictest, quartictest.
*/
{
   char buf[256];
   int ix0,ix1,ix2,ix3;
   int descartes(),ferrari(),neumark(),yacfraid(),chris();
   float fa,fb,fc,fd;
   float fp,fq,fr;
   double fabs();

   printf("docoeffsa %d\n",n);
   x0 = doub0;
   x1 = doub0;
   x2 = doub0;
   x3 = doub0;  
   if (NULL != gets(buf))
   {
      if (buf[0] == 'q')
         exit(0);
      if (n == 0)
      {
         x0 = doub1;
         if (docubic == TRUE)
         {
            sscanf(buf,"%g %g %g",&fp,&fq,&fr);
            p = (double)fp;
            q = (double)fq;
            r = (double)fr;
         }
         else
         {
            sscanf(buf,"%g %g %g %g",&fa,&fb,&fc,&fd);
            printf("\n%g %g %g %g\n",fa,fb,fc,fd);
            a = (double)fa;
            b = (double)fb;
            c = (double)fc;
            d = (double)fd;
            printf("\n%g %g %g %g\n",a,b,c,d);
         }
      }
      else
      if (n == 1) 
      {
         sscanf(buf,"%d\n",&ix0);
         x0 = (double)ix0;
      }
      else
      if (n == 2)
      {
         sscanf(buf,"%d %d\n",&ix0,&ix1);
         x0 = (double)ix0;
         x1 = (double)ix1;
      }
      else
      if (n == 3)
      {
         sscanf(buf,"%d %d %d\n",&ix0,&ix1,&ix2);
         x0 = (double)ix0;
         x1 = (double)ix1;
         x2 = (double)ix2;
      }
      else
      if (n == 4)
      {
         sscanf(buf,"%d %d %d %d\n",&ix0,&ix1,&ix2,&ix3);
         x0 = (double)ix0;
         x1 = (double)ix1;
         x2 = (double)ix2;
         x3 = (double)ix3;
      }
   }
   if ((fabs(x0)+fabs(x1)+fabs(x2)+fabs(x3)) != doub0)
   {
      if (docubic == TRUE) cubictest();
      if (doquartic == TRUE) quartictest();
   }
} /* docoeffs */
/*****************************************/

cubictest()
/*
   called by docoeffs.
*/
{
   int j;
   int nrtsc;
   double rtsc[4];

   if (n == 1)
   {
      p = -x0;
      q = doub1;
      r = -x0;
   }
   else
   if (n == 3)
   {
      p = -(x0+x1+x2);
      q = x0*x1 + x0*x2 + x1*x2;
      r = -x0*x1*x2;
   }
   printf("\ncubic test: %g %g %g\n",
         x0,x1,x2);
   printf(" x^3 + %gx^2 + %gx + %g\n",
         p,q,r);
   nrtsc = cubic(p,q,r,rtsc);
   printf("%d roots\n",nrtsc);
   for (j = 0; j < nrtsc; ++j)
      printf("%g  ",rtsc[j]);
   printf("\n");
} /* cubictest */
/******************************************/

quartictest()
/*
   called by docoeffs.
   calls descartes, ferrari, neumark, yacfraid, 
         errors, cases, setcns.
*/
{
   int j;

   printf("quartictest %d\n",n);
   if (n == 4)
   {
      a = -(x0 + x1 + x2 + x3);
      b = x0*x1 + x0*x2 + x0*x3 + x1*x2 + x1*x3 + x2*x3;
      c = -(x0*x1*x2 + x0*x1*x3 + x0*x2*x3 + x1*x2*x3);
      d = x0*x1*x2*x3;
   }
   else
   if (n == 2)
   {
      a = -(x0 + x1);
      b = x0*x1;
      c = a;
      d = b;
   }
   printf("\nquartic test: %g %g %g %g\n",
         x0,x1,x2,x3);
   printf("x^4 + %gx^3 + %gx^2 + %gx + %g\n",
         a,b,c,d);

   setcns();
   printf("\nDescartes\n");
   nrtsd = descartes(a,b,c,d,rtsd);
   worstd = errors(a,b,c,d,rtsd,rterd,nrtsd);
   if (nrtsd > 0)
         for (j = 0; j < nrtsd; ++j)
            printf(" %g %g\n",rtsd[j],rterd[j]);
   else 
         printf("no real roots found\n");

   setcns();
   printf("\nFerrari\n");
   nrtsf = ferrari(a,b,c,d,rtsf);
   worstf = errors(a,b,c,d,rtsf,rterf,nrtsf);
   if (nrtsf > 0)
         for (j = 0; j < nrtsf; ++j)
            printf(" %g %g\n",rtsf[j],rterf[j]);
   else 
         printf("no real roots found\n");

   setcns();
   printf("\nNeumark\n");
   nrtsn = neumark(a,b,c,d,rtsn);
   worstn = errors(a,b,c,d,rtsn,rtern,nrtsn);
   if (nrtsn > 0)
         for (j = 0; j < nrtsn; ++j)
            printf(" %g %g\n",rtsn[j],rtern[j]);
   else 
         printf("no real roots found\n");

   setcns();
   printf("\nYacoub and Fraidenraich\n");
   nrtsy = yacfraid(a,b,c,d,rtsy);
   worsty = errors(a,b,c,d,rtsy,rtery,nrtsy);
   if (nrtsy > 0)
        for (j = 0; j < nrtsy; ++j)
            printf(" %g %g\n",rtsy[j],rtery[j]);
   else 
         printf("no real roots found\n");

   setcns();
   printf("\nChristianson\n");
   nrtsc = chris(a,b,c,d,rtsc);
   worstc = errors(a,b,c,d,rtsc,rterc,nrtsc);
   if (nrtsc > 0)
         for (j = 0; j < nrtsc; ++j)
            printf(" %g %g\n",rtsc[j],rterc[j]);
   else 
         printf("no real roots found\n");
   printf("\n x^4 + %gx^3 + %gx^2 + %gx + %g\n",
         a,b,c,d);
   if (debug < 2) cases();
} /* quartictest */
/***************************************/

compare()
/*
    called by looptest.
    calls errors, descartes, ferrari, neumark, yacfraid, chris.
*/
{
   int n;
   double errors();

   ++tot;
   nrtsd = descartes(a,b,c,d,rtsd);
   worstd = errors(a,b,c,d,rtsd,rterd,nrtsd);
   nrtsf = ferrari(a,b,c,d,rtsf);
   worstf = errors(a,b,c,d,rtsf,rterf,nrtsf);
   nrtsn = neumark(a,b,c,d,rtsn,worstn);
   worstn = errors(a,b,c,d,rtsn,rtern,nrtsn);
   nrtsy = yacfraid(a,b,c,d,rtsy);
   worsty = errors(a,b,c,d,rtsy,rtery,nrtsy);
   nrtsc = chris(a,b,c,d,rtsc);
   worstc = errors(a,b,c,d,rtsc,rterc,nrtsc);

   if (nrtsd == 0) ++nzd;
   if (nrtsf == 0) ++nzf;
   if (nrtsn == 0) ++nzn;
   if (nrtsy == 0) ++nzy;
   if (nrtsc == 0) ++nzc;

   if (nrtsc == nrtsd) ++agreecd;
   if (nrtsc == nrtsf) ++agreecf;
   if (nrtsc == nrtsn) ++agreecn;
   if (nrtsc == nrtsy) ++agreecy;
   if (nrtsd == nrtsf) ++agreedf;
   if (nrtsd == nrtsn) ++agreedn;
   if (nrtsd == nrtsy) ++agreedy;
   if (nrtsf == nrtsn) ++agreefn;
   if (nrtsf == nrtsy) ++agreefy;
   if (nrtsn == nrtsy) ++agreeny;

   if ((nrtsc == nrtsd) && (nrtsc == nrtsf)
    && (nrtsc == nrtsn) && (nrtsc == nrtsy))
   {
      if (nrtsf == 0) ++zer;
      ++agr;
   }
   else
      ++dis;

     
   if ((nrtsc == nrtsd) && (nrtsc == nrtsf)
    && (nrtsc == nrtsn) && (nrtsc == nrtsy) && (nrtsc != 0))
   {
       if ((worstd < worstf) && (worstd < worstn)
        && (worstd < worsty) && (worstd < worstc)) ++nd;

       if ((worstf < worstd) && (worstf < worstn)
        && (worstf < worsty) && (worstf < worstc)) ++nf;

       if ((worstn < worstd) && (worstn < worstf) 
        && (worstn < worsty) && (worstn < worstc)) ++nn;

       if ((worsty < worstd) && (worsty < worstf)
        && (worsty < worstn) && (worsty < worstc)) ++ny;

       if ((worstc < worstd) && (worstc < worstf) 
        && (worstc < worstn) && (worstc < worsty)) ++nc;

       if (maxd < worstd)
       {
          maxd = worstd;
          cd[0]=a; cd[1]=b; cd[2]=c; cd[3]=d;
       }
       if (maxf < worstf)
       {
          maxf = worstf;
          cf[0]=a; cf[1]=b; cf[2]=c; cf[3]=d;
       }
       if (maxn < worstn)
       {
          maxn = worstn;
          cn[0]=a; cn[1]=b; cn[2]=c; cn[3]=d;
       }
       if (maxy < worsty)
       {
          maxy = worsty;
          cy[0]=a; cy[1]=b; cy[2]=c; cy[3]=d;
       }
       if (maxc < worstc)
       {
           maxc = worstc;
           cc[0]=a; cc[1]=b; cc[2]=c; cc[3]=d;
       }

   } /* nrts agree */
} /* compare */
/****************************************************/

double errors(a,b,c,d,rts,rterr,nrts)
double a,b,c,d,rts[4],rterr[4];
int nrts;
/*
   find the errors

   called by quartictest, docoeff, compare,
             chris, descartes, ferrari, neumark, yacfraid.
*/
{
   int k;
   double deriv,test,worst;
   double fabs(),sqrt(),curoot();

   worst = doubmax;
   if (nrts > 0)
   {
      worst =  doub0;
      for (  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] = fabs(test/deriv);
            else
            {
               deriv = (doub12*rts[k]+doub6*a)*rts[k]+doub2*b;
               if (deriv != doub0)
                   rterr[k] = sqrt(fabs(test/deriv));
               else
               {
                  deriv = doub24*rts[k]+doub6*a;
                  if (deriv != doub0)
                     rterr[k] = curoot(fabs(test/deriv));
                  else
                     rterr[k] = sqrt(sqrt(fabs(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 */
/**********************************************/

double acos3(x)
   double x ;
/* 
     find cos(acos(x)/3) 
    
     16 Jul 1981   Don Herbison-Evans 

     called by cubic . 
*/
{
   double value;
   double acos(),cos();

   value = cos(acos(x)*inv3);
   return(value);
} /* acos3 */
/***************************************/

double curoot(x)
   double x ;
/* 
     find cube root of x.

     30 Jan 1989   Don Herbison-Evans 

     called by cubic . 
*/
{
   double exp(),log();
   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 */
/****************************************************/

int quadratic(b,c,rts)
   double b,c,rts[4];
/* 
     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 j,nquad;
   double dis,rtdis ;
   double ans[2] ;

   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 */
/**************************************************/

int cubic(p,q,r,v3)
double p,q,r,v3[4];
/* 
   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    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;
   double ans;
   double curoot();
   double acos3();
   double sqrt(),fabs();
   int quadratic();

   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 */
/****************************************************/

cubnewton(p,q,r,n3,v3)
int n3;
double p,q,r,v3[4];
/*
    improve roots of cubic by Newton-Raphson iteration

    5 Jan 2004  Don Herbison-Evans

    called by cubic.
*/
{
   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 */
/****************************************************/

int quartic(a,b,c,d,rts)
double a,b,c,d,rts[4];
/*
   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 quadratic(),ferrari(),neumark(),yacfraid();
   int j,k,nq,nr;
   double odd, even;
   double roots[4];
   double fabs();

   if (fabs(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 */
/*****************************************/

int descartes(a,b,c,d,rts)
double a,b,c,d,rts[4];
/*
   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
*/
{
   double errors();
   int quadratic(),cubic();

   int j;
   double v1[4],v2[4],v3[4];
   double k,y;
   double dis;
   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 */
/****************************************************/

int ferrari(a,b,c,d,rts)
   double a,b,c,d,rts[4];
/* 
     solve the quartic equation - 

     x**4 + a*x**3 + b*x**2 + c*x + d = 0 

     called by quartic, compare, quartictest.
     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 
*/
{
   double errors();
   int cubic(),quadratic();

   int j,k;
   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 */
/*****************************************/

int neumark(a,b,c,d,rts)
   double a,b,c,d,rts[4];
/* 
     solve the quartic equation - 

     x**4 + a*x**3 + b*x**2 + c*x + d = 0 

     called by quartic, compare, quartictest.
     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 j,k;
   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;
   double errors();
   int cubic();
   int quadratic();

   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 */
/****************************************************/

int yacfraid(a,b,c,d,rts)
   double a,b,c,d,rts[4];
/* 
     solve the quartic equation - 

     x**4 + a*x**3 + b*x**2 + c*x + d = 0 

     called by quartic, compare, quartictest.
     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 i,j;
   double y;
   double v1[4],v2[4],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;
   double errors();
   int cubic();
   int quadratic();

   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 */
/*****************************************/

int chris(a,b,c,d,rts)
double a,b,c,d,rts[4];
/*
   Solve quartic equation using
   Christianson's algorithm

   called by compare, quartictest.
   calls errors, quadratic, cubic.

   Bruce Christianson, Solving Quartics Using Palindromes,
      Mathematical Gazette, Vol. 75, pp. 327-328 (1991)

   14 Jan 2004  Don Herbison-Evans
*/
{
   double errors();
   int quadratic(),cubic();

   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(fabs(y*ycu + e2*ysq + e1*y + e0)));
         if (y == doub0)
             k1 = doub0;
         else
             k1 = sqrt(fabs(ysq + inv2*e2 + inv4*e1/y));
         g1 = (doub4*ycu + doub2* y*e2 + e1)*k;
         printf("chrisc %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 */