quarcube.c

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

/*
   doquartic.c

   a test of a quartic solving routine

   Don Herbison-Evans   24 June 1994
*/

#include <stdio.h>

double nought;
double doub1,doub2;
double doub3,doub4;
double doub6,doub12;
double doub24;
double doubmax;       /* approx square root of max double number */
double doubmin;       /* smallest double number */
double doubtol;       /* tolerance of double numbers */
double rt3;
double inv2,inv3,inv4;

void setcns();
int descartes();
double exp(),log(),sqrt(),cos(),acos();
double cubic();
int qudrtc();
int quartic();
int ferrari();
int neumark();
void errors();

main()
{
   double a,b,c,d;
   double rts[4];
   double rterr[4];
   int nr;

   setcns();
   a = -(double)10;
   b = (double)35;
   c = -(double)50;
   d = (double)24;
   nr = descartes(a,b,c,d,rts);
   errors(a,b,c,d,rts,rterr,nr);
   nr = ferrari(a,b,c,d,rts);
   errors(a,b,c,d,rts,rterr,nr);
   nr = neumark(a,b,c,d,rts);
   errors(a,b,c,d,rts,rterr,nr);
   exit(0);
}
/****************************************************/

void setcns()
/* 
     set up constants 
*/
{
      int j;

      nought = (double)0;
      doub1 = (double)1;
      doub2 = (double)2;
      doub3 = (double)3;
      doub4 = (double)4;
      doub6 = (double)6;
      doub12 = (double)12;
      doub24 = (double)24;
      inv2 = doub1/(double)2;
      inv3 = doub1/(double)3;
      inv4 = doub1/(double)4;
      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=0.7/sqrt(doubmin) ;

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

int quartic(a,b,c,d,rts,rterr)
double a,b,c,d,rts[4],rterr[4];
/*
   Solve quartic equation using either
   quadratic, Ferrari's or Neumark's algorithm.

   called by 
   calls  qudrtc, ferrari, neumark.

     21 Jan 1989  Don Herbison-Evans
*/
{
   int qudrtc(),ferrari(),neumark();
   int j,k,nq,nr;
   double odd, even;
   double roots[4];

   if (a < nought) odd = -a; else odd = a;
   if (c < nought) odd -= c; else odd += c;
   if (b < nought) even = -b; else even = b;
   if (d < nought) even -= d; else even += d;
   if (odd < even*doubtol)
   {
      nq = qudrtc(b,d,roots,b*b-doub4*d);
      j = 0;
      for (k = 0; k < nq; ++k)
      {
         if (roots[k] > nought)
         {
            rts[j] = sqrt(roots[k]);
            rts[j+1] = -rts[j];
            ++j; ++j;
         }
      }
      nr = j;
   }
   else
   {
      if (a < nought) k = 1; else k = 0;
      if (b < nought) k += k+1; else k +=k; 
      if (c < nought) k += k+1; else k +=k; 
      if (d < nought) k += k+1; else k +=k; 
      switch (k)
      {
              case 0 : nr = ferrari(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 = ferrari(a,b,c,d,rts) ; break;
              case 5 : nr = neumark(a,b,c,d,rts) ; break;
              case 6 : nr = ferrari(a,b,c,d,rts) ; break;
              case 7 : nr = ferrari(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 = ferrari(a,b,c,d,rts) ; break;
              case 11 : nr = neumark(a,b,c,d,rts) ; break;
              case 12 : nr = ferrari(a,b,c,d,rts) ; break;
              case 13 : nr = ferrari(a,b,c,d,rts) ; break;
              case 14 : nr = ferrari(a,b,c,d,rts) ; break;
              case 15 : nr = ferrari(a,b,c,d,rts) ; break;
      }
   }
   errors(a,b,c,d,rts,rterr,nr);
   return(nr);
} /* quartic */
/*****************************************/

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

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

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

     method :  Ferrari - Lagrange
     Theory of Equations, H.W. Turnbull p. 140 (1947)

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

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

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

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

   esq = inv4*asq - b - y;
   if (esq < nought) return(0);
   else
   {
      fsq = inv4*y*y - d;
      if (fsq < nought) return(0);
      else
      {
         ef = -(inv4*a*y + inv2*c);
         if ( ((a > nought)&&(y > nought)&&(c > nought))
           || ((a > nought)&&(y < nought)&&(c < nought))
           || ((a < nought)&&(y > nought)&&(c < nought))
           || ((a < nought)&&(y < nought)&&(c > nought))
           ||  (a == nought)||(y == nought)||(c == nought)
         )
/* use ef - */
         {
            if ((b < nought)&&(y < nought)&&(esq > nought))
            {
               e = sqrt(esq);
               f = ef/e;
            }
            else if ((d < nought) && (fsq > nought))
            {
               f = sqrt(fsq);
               e = ef/f;
            }
            else
            {
               e = sqrt(esq);
               f = sqrt(fsq);
               if (ef < nought) f = -f;
            }
         }
         else
         {
            e = sqrt(esq);
            f = sqrt(fsq);
            if (ef < nought) f = -f;
         }
/* note that e >= nought */
         ainv2 = a*inv2;
         g = ainv2 - e;
         gg = ainv2 + e;
         if ( ((b > nought)&&(y > nought))
           || ((b < nought)&&(y < nought)) )
         {
            if (( a > nought) && (e != nought)) g = (b + y)/gg;
               else if (e != nought) gg = (b + y)/g;
         }
         if ((y == nought)&&(f == nought))
         {
            h = nought;
            hh = nought;
         }
         else if ( ((f > nought)&&(y < nought))
                || ((f < nought)&&(y > nought)) )
         {
            hh = -inv2*y + f;
            h = d/hh;
         }
         else
         {
            h = -inv2*y - f;
            hh = d/h;
         }
         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);
      }
   }
} /* 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
   calls     cubic, qudrtc.

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

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

     method -  S. Neumark 

     Solution of Cubic and Quartic Equations - Pergamon 1965 
        translated to C with help of Shawn Neely

*/
{
   int nquar,n1,n2 ;
   double y,g,gg,h,hh,gdis,gdisrt,hdis,hdisrt,g1,g2,h1,h2 ;
   double bmy,gerr,herr,y4,d4,bmysq ;
   double v1[4],v2[4] ;
   double asq ;
   double p,q,r ;
   double hmax,gmax ;
   double cubic();
   int qudrtc();

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

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

   bmy = b - y ;
   y4 = y*doub4 ;
   d4 = d*doub4 ;
   bmysq = bmy*bmy ;
   gdis = asq - y4 ;
   hdis = bmysq - d4 ;
   if ((gdis < nought) || (hdis < nought)) return(0);
   else
   {
      g1 = a*inv2 ;
      h1 = bmy*inv2 ;
      gerr = asq + y4 ;
      herr = hdis ;
      if (d > nought) herr = bmysq + d4 ;
      if ((y < nought) || (herr*gdis > gerr*hdis))
      {
         gdisrt = sqrt(gdis) ;
         g2 = gdisrt*inv2 ;
         if (gdisrt != nought) h2 = (a*h1 - c)/gdisrt ;
            else h2 = nought;
      }
      else
      {
         hdisrt = sqrt(hdis) ;
         h2 = hdisrt*inv2 ;
         if (hdisrt != nought) g2 = (a*h1 - c)/hdisrt ;
            else g2 = nought;
      }
/* 
     note that in the following, the tests ensure non-zero 
     denominators -  
*/
      h = h1 - h2 ;
      hh = h1 + h2 ;
      hmax = hh ;
      if (hmax < nought) hmax =  -hmax ;
      if (hmax < h) hmax = h ;
      if (hmax <  -h) hmax =  -h ;
      if ((h1 > nought)&&(h2 > nought)) h = d/hh ;
      if ((h1 < nought)&&(h2 < nought)) h = d/hh ;
      if ((h1 > nought)&&(h2 < nought)) hh = d/h ;
      if ((h1 < nought)&&(h2 > nought)) hh = d/h ;
      if (h > hmax) h = hmax ;
      if (h <  -hmax) h =  -hmax ;
      if (hh > hmax) hh = hmax ;
      if (hh <  -hmax) hh =  -hmax ;

      g = g1 - g2 ;
      gg = g1 + g2 ;
      gmax = gg ;
      if (gmax < nought) gmax =  -gmax ;
      if (gmax < g) gmax = g ;
      if (gmax <  -g) gmax =  -g ;
      if ((g1 > nought)&&(g2 > nought)) g = y/gg ;
      if ((g1 < nought)&&(g2 < nought)) g = y/gg ;
      if ((g1 > nought)&&(g2 < nought)) gg = y/g ;
      if ((g1 < nought)&&(g2 > nought)) gg = y/g ;
      if (g > gmax) g = gmax ;
      if (g <  -gmax) g =  -gmax ;
      if (gg > gmax) gg = gmax ;
      if (gg <  -gmax) gg =  -gmax ;
 
      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);
   }
} /* neumark */
/****************************************************/

void errors(a,b,c,d,rts,rterr,nrts)