quarcube.c

[Game.Programming].Academic - Graphics Gems (6 books source code)

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

   called by quartic.
*/
{
   int k;
   double deriv,test;
   double fabs(),sqrt(),curoot();

   if (nrts > 0)
   {
      for (  k = 0 ; k < nrts ; ++ k )
      {
         test = (((rts[k]+a)*rts[k]+b)*rts[k]+c)*rts[k]+d ;
         if (test == nought) rterr[k] = nought;
         else
         {
            deriv =
               ((doub4*rts[k]+doub3*a)*rts[k]+doub2*b)*rts[k]+c ;
            if (deriv != nought)
               rterr[k] = fabs(test/deriv);
            else
            {
               deriv = (doub12*rts[k]+doub6*a)*rts[k]+doub2*b ;
               if (deriv != nought)
                   rterr[k] = sqrt(fabs(test/deriv)) ;
               else
               {
                  deriv = doub24*rts[k]+doub6*a ;
                  if (deriv != nought)
                     rterr[k] = curoot(fabs(test/deriv));
                  else
                     rterr[k] = sqrt(sqrt(fabs(test)/doub24));
               }
            }
         }
         fprintf(stderr,"errorsa  %d %9g %9g\n",
               k,rts[k],rterr[k]);
      }
   }
   else fprintf(stderr,"errors ans: none\n");
} /* errors */
/**********************************************/

int qudrtc(b,c,rts,dis)
   double b,c,rts[4],dis ;
/* 
     solve the quadratic equation - 

         x**2+b*x+c = 0 

     called by  quartic, ferrari, neumark, ellcut 
*/
{
   int nquad;
   double rtdis ;

   if (dis >= nought)
   {
      nquad = 2 ;
      rtdis = sqrt(dis) ;
      if (b > nought) rts[0] = ( -b - rtdis)*inv2 ;
         else rts[0] = ( -b + rtdis)*inv2 ;
      if (rts[0] == nought) rts[1] =  -b ;
      else rts[1] = c/rts[0] ;
   }
   else
   {
      nquad = 0;
      rts[0] = nought ;
      rts[1] = nought ;
   }
   return(nquad);
} /* qudrtc */
/**************************************************/

double cubic(p,q,r)
double p,q,r;
/* 
     find the lowest real root of the cubic - 
       x**3 + p*x**2 + q*x + r = 0 

   input parameters - 
     p,q,r - coeffs of cubic equation. 

   output- 
     cubic - a real root. 

   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 

     Charles Prineas   April 1981 

     called by  neumark.
     calls  acos3 
*/
{
   int nrts;
   double po3,po3sq,qo3;
   double uo3,u2o3,uo3sq4,uo3cu4 ;
   double v,vsq,wsq ;
   double m,mcube,n;
   double muo3,s,scube,t,cosk,sinsqk ;
   double root;
   double curoot();
   double acos3();
   double sqrt(),fabs();

   m = nought;
   nrts =0;
   if ((p > doubmax) || (p <  -doubmax)) root = -p;
   else
   if ((q > doubmax) || (q <  -doubmax))
   {
      if (q > nought) root = -r/q;
      else root = -sqrt(-q);
   }
   else
   if ((r > doubmax)|| (r <  -doubmax)) root =  -curoot(r) ;
   else
   {
      po3 = p*inv3 ;
      po3sq = po3*po3 ;
      if (po3sq > doubmax) root =  -p ;
      else
      {
         v = r + po3*(po3sq + po3sq - q) ;
         if ((v > doubmax) || (v < -doubmax)) root = -p ;
         else
         {
            vsq = v*v ;
            qo3 = q*inv3 ;
            uo3 = qo3 - po3sq ;
            u2o3 = uo3 + uo3 ;
            if ((u2o3 > doubmax) || (u2o3 < -doubmax))
            {
               if (p == nought)
               {
                  if (q > nought) root =  -r/q ;
                  else root =  -sqrt(-q) ;
               }
               else root =  -q/p ;
            }
            uo3sq4 = u2o3*u2o3 ;
            if (uo3sq4 > doubmax)
            {
               if (p == nought)
               {
                  if (q > nought) root = -r/q ;
                  else root = -sqrt(fabs(q)) ;
               }
               else root = -q/p ;
            }
            uo3cu4 = uo3sq4*uo3 ;
            wsq = uo3cu4 + vsq ;
            if (wsq >= nought)
            {
/* 
     cubic has one real root 
*/
               nrts = 1;
               if (v <= nought) mcube = ( -v + sqrt(wsq))*inv2 ;
               if (v  > nought) mcube = ( -v - sqrt(wsq))*inv2 ;
               m = curoot(mcube) ;
               if (m != nought) n = -uo3/m ;
                  else n = nought;
               root = m + n - po3 ;
            }
            else
            {
               nrts = 3;
/* 
     cubic has three real roots 
*/
               if (uo3 < nought)
               {
                  muo3 = -uo3;
                  s = sqrt(muo3) ;
                  scube = s*muo3;
                  t =  -v/(scube+scube) ;
                  cosk = acos3(t) ;
                  if (po3 < nought)
                     root = (s+s)*cosk - po3;
                  else
                  {
                     sinsqk = doub1 - cosk*cosk ;
                     if (sinsqk < nought) sinsqk = nought ;
                     root = s*( -cosk - rt3*sqrt(sinsqk)) - po3 ;
                  }
               }
               else
/* 
     cubic has multiple root -  
*/
               root = curoot(v) - po3 ;
            }
         }
      }
   }
   fprintf(stderr,"cubic %g %g %g %d %g\n",p,q,r,nrts,root);
   return(root);
} /* cubic */
/***************************************/

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

     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.
 
     Don Herbison-Evans   30/1/89 

     called by cubic . 
*/
{
   double exp(),log();
   double value;
   double absx;
   int neg;

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

int simple(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
   calls     cubic, qudrtc.

     input - 
   a,b,c,e - coeffs of equation. 

     output - 
   nquar - number of real roots. 
   rts - array of root values. 

     method :  unstabilized Ferrari-Lagrange
     Abramowitz,M. & Stegun I.A.
     Handbook of Mathematical Functions
     Dover 1972 (ninth printing), pp. 17-18

     calls  cubic, qudrtc 
*/
{
   double cubic();
   int qudrtc();

   int nquar,n1,n2 ;
   double asq,y;
   double v1[4],v2[4] ;
   double p,q,r ;
   double e,f,esq,fsq ;
   double g,gg,h,hh;

   fprintf(stderr,"\nsimple %g %g %g %g\n",a,b,c,d);
   asq = a*a;

   p = -b ;
   q = a*c-doub4*d ;
   r = -asq*d - c*c + doub4*b*d ;
   y = cubic(p,q,r) ;

   esq = inv4*asq - b + y;
   fsq = inv4*y*y - d;
   if (esq < nought) return(0);
   else
   if (fsq < nought) return(0);
   else
   {
      e = sqrt(esq);
      f = sqrt(fsq);
      g = inv2*a - e;
      h = inv2*y - f;
      gg = inv2*a + e;
      hh = inv2*y + f;
      n1 = qudrtc(gg,hh,v1, gg*gg - doub4*hh) ;
      n2 = qudrtc(g,h,v2, g*g - doub4*h) ;
      nquar = n1+n2 ;
      rts[0] = v1[0] ;
      rts[1] = v1[1] ;
      rts[n1+0] = v2[0] ;
      rts[n1+1] = v2[1] ;
      return(nquar);
   }
} /* simple */
/*****************************************/
 
int descartes(a,b,c,d,rts)
double a,b,c,d,rts[4];
/*
   Solve quartic equation using
   Descartes-Euler-Cardano algorithm

   Strong, T. "Elemementary and Higher Algebra"
      Pratt and Oakley, p. 469 (1859)

     29 Jun 1994  Don Herbison-Evans
*/
{
   int qudrtc();
   double cubic();

   int nrts;
   int r1,r2;
   double v1[4],v2[4];
   double y;
   double p,q,r;
   double A,B,C;
   double m,n1,n2;
   double d3o8,d3o256;
   double inv8,inv16;
   double asq;
   double Binvm;

   d3o8 = (double)3/(double)8;
   inv8 = doub1/(double)8;
   inv16 = doub1/(double)16;
   d3o256 = (double)3/(double)256;

   fprintf(stderr,"\nDescartes %f %f %f %f\n",a,b,c,d);
   asq = a*a;

   A = b - asq*d3o8;
   B = c + a*(asq*inv8 - b*inv2);
   C = d + asq*(b*inv16 - asq*d3o256) - a*c*inv4;

   p = doub2*A;
   q = A*A - doub4*C;
   r = -B*B;

/***
   inv64 = doub1/(double)64;
   p = doub2*b - doub3*a*a*inv4 ;
   q = b*b - a*a*b - doub4*d + doub3*a*a*a*a*inv16 + a*c;
   r = -c*c - a*a*a*a*a*a*inv64 - a*a*b*b*inv4
       -a*a*a*c*inv4 + a*b*c + a*a*a*a*b*inv8;
***/

   y = cubic(p,q,r) ;
   if (y <= nought) 
      nrts = 0;
   else
   {
      m = sqrt(y);
      Binvm = B/m;
      n1 = (y + A + Binvm)*inv2;
      n2 = (y + A - Binvm)*inv2;
      r1 = qudrtc(-m, n1, v1, y-doub4*n1);
      r2 = qudrtc( m, n2, v2, y-doub4*n2);
      rts[0] = v1[0]-a*inv4;
      rts[1] = v1[1]-a*inv4;
      rts[r1] = v2[0]-a*inv4;
      rts[r1+1] = v2[1]-a*inv4;
      nrts = r1+r2;
   } 
   return(nrts);
} /* descartes */
/****************************************************/