home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch1_1 / quarcube.c
Encoding:
C/C++ Source or Header  |  1995-04-04  |  18.3 KB  |  796 lines
  1. /*
  2.    doquartic.c
  3.  
  4.    a test of a quartic solving routine
  5.  
  6.    Don Herbison-Evans   24 June 1994
  7. */
  8.  
  9. #include <stdio.h>
  10.  
  11. double nought;
  12. double doub1,doub2;
  13. double doub3,doub4;
  14. double doub6,doub12;
  15. double doub24;
  16. double doubmax;       /* approx square root of max double number */
  17. double doubmin;       /* smallest double number */
  18. double doubtol;       /* tolerance of double numbers */
  19. double rt3;
  20. double inv2,inv3,inv4;
  21.  
  22. void setcns();
  23. int descartes();
  24. double exp(),log(),sqrt(),cos(),acos();
  25. double cubic();
  26. int qudrtc();
  27. int quartic();
  28. int ferrari();
  29. int neumark();
  30. void errors();
  31.  
  32. main()
  33. {
  34.    double a,b,c,d;
  35.    double rts[4];
  36.    double rterr[4];
  37.    int nr;
  38.  
  39.    setcns();
  40.    a = -(double)10;
  41.    b = (double)35;
  42.    c = -(double)50;
  43.    d = (double)24;
  44.    nr = descartes(a,b,c,d,rts);
  45.    errors(a,b,c,d,rts,rterr,nr);
  46.    nr = ferrari(a,b,c,d,rts);
  47.    errors(a,b,c,d,rts,rterr,nr);
  48.    nr = neumark(a,b,c,d,rts);
  49.    errors(a,b,c,d,rts,rterr,nr);
  50.    exit(0);
  51. }
  52. /****************************************************/
  53.  
  54. void setcns()
  55. /* 
  56.      set up constants 
  57. */
  58. {
  59.       int j;
  60.  
  61.       nought = (double)0;
  62.       doub1 = (double)1;
  63.       doub2 = (double)2;
  64.       doub3 = (double)3;
  65.       doub4 = (double)4;
  66.       doub6 = (double)6;
  67.       doub12 = (double)12;
  68.       doub24 = (double)24;
  69.       inv2 = doub1/(double)2;
  70.       inv3 = doub1/(double)3;
  71.       inv4 = doub1/(double)4;
  72.       rt3 = sqrt(doub3) ;
  73.  
  74.       doubtol = doub1;
  75.       for (  j = 1 ; doub1+doubtol > doub1 ; ++ j )
  76.       {
  77.           doubtol *= inv2;
  78.       }
  79.       doubtol = sqrt(doubtol);
  80.  
  81.       doubmin = inv2 ;
  82.       for (  j = 1 ; j <= 100 ; ++ j )
  83.       {
  84.           doubmin=doubmin*doubmin ;
  85.           if ((doubmin*doubmin) <= (doubmin*doubmin*inv2))
  86.               break;
  87.       }
  88.       doubmax=0.7/sqrt(doubmin) ;
  89.  
  90. } /* setcns */
  91. /***********************************/
  92.  
  93. int quartic(a,b,c,d,rts,rterr)
  94. double a,b,c,d,rts[4],rterr[4];
  95. /*
  96.    Solve quartic equation using either
  97.    quadratic, Ferrari's or Neumark's algorithm.
  98.  
  99.    called by 
  100.    calls  qudrtc, ferrari, neumark.
  101.  
  102.      21 Jan 1989  Don Herbison-Evans
  103. */
  104. {
  105.    int qudrtc(),ferrari(),neumark();
  106.    int j,k,nq,nr;
  107.    double odd, even;
  108.    double roots[4];
  109.  
  110.    if (a < nought) odd = -a; else odd = a;
  111.    if (c < nought) odd -= c; else odd += c;
  112.    if (b < nought) even = -b; else even = b;
  113.    if (d < nought) even -= d; else even += d;
  114.    if (odd < even*doubtol)
  115.    {
  116.       nq = qudrtc(b,d,roots,b*b-doub4*d);
  117.       j = 0;
  118.       for (k = 0; k < nq; ++k)
  119.       {
  120.          if (roots[k] > nought)
  121.          {
  122.             rts[j] = sqrt(roots[k]);
  123.             rts[j+1] = -rts[j];
  124.             ++j; ++j;
  125.          }
  126.       }
  127.       nr = j;
  128.    }
  129.    else
  130.    {
  131.       if (a < nought) k = 1; else k = 0;
  132.       if (b < nought) k += k+1; else k +=k; 
  133.       if (c < nought) k += k+1; else k +=k; 
  134.       if (d < nought) k += k+1; else k +=k; 
  135.       switch (k)
  136.       {
  137.               case 0 : nr = ferrari(a,b,c,d,rts) ; break;
  138.               case 1 : nr = neumark(a,b,c,d,rts) ; break;
  139.               case 2 : nr = neumark(a,b,c,d,rts) ; break;
  140.               case 3 : nr = ferrari(a,b,c,d,rts) ; break;
  141.               case 4 : nr = ferrari(a,b,c,d,rts) ; break;
  142.               case 5 : nr = neumark(a,b,c,d,rts) ; break;
  143.               case 6 : nr = ferrari(a,b,c,d,rts) ; break;
  144.               case 7 : nr = ferrari(a,b,c,d,rts) ; break;
  145.               case 8 : nr = neumark(a,b,c,d,rts) ; break;
  146.               case 9 : nr = ferrari(a,b,c,d,rts) ; break;
  147.               case 10 : nr = ferrari(a,b,c,d,rts) ; break;
  148.               case 11 : nr = neumark(a,b,c,d,rts) ; break;
  149.               case 12 : nr = ferrari(a,b,c,d,rts) ; break;
  150.               case 13 : nr = ferrari(a,b,c,d,rts) ; break;
  151.               case 14 : nr = ferrari(a,b,c,d,rts) ; break;
  152.               case 15 : nr = ferrari(a,b,c,d,rts) ; break;
  153.       }
  154.    }
  155.    errors(a,b,c,d,rts,rterr,nr);
  156.    return(nr);
  157. } /* quartic */
  158. /*****************************************/
  159.  
  160. int ferrari(a,b,c,d,rts)
  161.    double a,b,c,d,rts[4];
  162. /* 
  163.      solve the quartic equation - 
  164.  
  165.    x**4 + a*x**3 + b*x**2 + c*x + d = 0 
  166.  
  167.    called by quartic
  168.    calls     cubic, qudrtc.
  169.  
  170.      input - 
  171.    a,b,c,e - coeffs of equation. 
  172.  
  173.      output - 
  174.    nquar - number of real roots. 
  175.    rts - array of root values. 
  176.  
  177.      method :  Ferrari - Lagrange
  178.      Theory of Equations, H.W. Turnbull p. 140 (1947)
  179.  
  180.      calls  cubic, qudrtc 
  181. */
  182. {
  183.    double cubic();
  184.    int qudrtc();
  185.  
  186.    int nquar,n1,n2 ;
  187.    double asq,ainv2;
  188.    double v1[4],v2[4] ;
  189.    double p,q,r ;
  190.    double y;
  191.    double e,f,esq,fsq,ef ;
  192.    double g,gg,h,hh;
  193.  
  194.    fprintf(stderr,"\nFerrari %g %g %g %g\n",a,b,c,d);
  195.    asq = a*a;
  196.  
  197.    p = b ;
  198.    q = a*c-doub4*d ;
  199.    r = (asq - doub4*b)*d + c*c ;
  200.    y = cubic(p,q,r) ;
  201.  
  202.    esq = inv4*asq - b - y;
  203.    if (esq < nought) return(0);
  204.    else
  205.    {
  206.       fsq = inv4*y*y - d;
  207.       if (fsq < nought) return(0);
  208.       else
  209.       {
  210.          ef = -(inv4*a*y + inv2*c);
  211.          if ( ((a > nought)&&(y > nought)&&(c > nought))
  212.            || ((a > nought)&&(y < nought)&&(c < nought))
  213.            || ((a < nought)&&(y > nought)&&(c < nought))
  214.            || ((a < nought)&&(y < nought)&&(c > nought))
  215.            ||  (a == nought)||(y == nought)||(c == nought)
  216.          )
  217. /* use ef - */
  218.          {
  219.             if ((b < nought)&&(y < nought)&&(esq > nought))
  220.             {
  221.                e = sqrt(esq);
  222.                f = ef/e;
  223.             }
  224.             else if ((d < nought) && (fsq > nought))
  225.             {
  226.                f = sqrt(fsq);
  227.                e = ef/f;
  228.             }
  229.             else
  230.             {
  231.                e = sqrt(esq);
  232.                f = sqrt(fsq);
  233.                if (ef < nought) f = -f;
  234.             }
  235.          }
  236.          else
  237.          {
  238.             e = sqrt(esq);
  239.             f = sqrt(fsq);
  240.             if (ef < nought) f = -f;
  241.          }
  242. /* note that e >= nought */
  243.          ainv2 = a*inv2;
  244.          g = ainv2 - e;
  245.          gg = ainv2 + e;
  246.          if ( ((b > nought)&&(y > nought))
  247.            || ((b < nought)&&(y < nought)) )
  248.          {
  249.             if (( a > nought) && (e != nought)) g = (b + y)/gg;
  250.                else if (e != nought) gg = (b + y)/g;
  251.          }
  252.          if ((y == nought)&&(f == nought))
  253.          {
  254.             h = nought;
  255.             hh = nought;
  256.          }
  257.          else if ( ((f > nought)&&(y < nought))
  258.                 || ((f < nought)&&(y > nought)) )
  259.          {
  260.             hh = -inv2*y + f;
  261.             h = d/hh;
  262.          }
  263.          else
  264.          {
  265.             h = -inv2*y - f;
  266.             hh = d/h;
  267.          }
  268.          n1 = qudrtc(gg,hh,v1, gg*gg - doub4*hh) ;
  269.          n2 = qudrtc(g,h,v2, g*g - doub4*h) ;
  270.          nquar = n1+n2 ;
  271.          rts[0] = v1[0] ;
  272.          rts[1] = v1[1] ;
  273.          rts[n1+0] = v2[0] ;
  274.          rts[n1+1] = v2[1] ;
  275.          return(nquar);
  276.       }
  277.    }
  278. } /* ferrari */
  279. /*****************************************/
  280.  
  281. int neumark(a,b,c,d,rts)
  282.    double a,b,c,d,rts[4];
  283. /* 
  284.      solve the quartic equation - 
  285.  
  286.    x**4 + a*x**3 + b*x**2 + c*x + d = 0 
  287.  
  288.    called by quartic
  289.    calls     cubic, qudrtc.
  290.  
  291.      input parameters - 
  292.    a,b,c,e - coeffs of equation. 
  293.  
  294.      output parameters - 
  295.    nquar - number of real roots. 
  296.    rts - array of root values. 
  297.  
  298.      method -  S. Neumark 
  299.  
  300.      Solution of Cubic and Quartic Equations - Pergamon 1965 
  301.         translated to C with help of Shawn Neely
  302.  
  303. */
  304. {
  305.    int nquar,n1,n2 ;
  306.    double y,g,gg,h,hh,gdis,gdisrt,hdis,hdisrt,g1,g2,h1,h2 ;
  307.    double bmy,gerr,herr,y4,d4,bmysq ;
  308.    double v1[4],v2[4] ;
  309.    double asq ;
  310.    double p,q,r ;
  311.    double hmax,gmax ;
  312.    double cubic();
  313.    int qudrtc();
  314.  
  315.    fprintf(stderr,"\nNeumark %g %g %g %g\n",a,b,c,d);
  316.    asq = a*a ;
  317.  
  318.    p =  -b*doub2 ;
  319.    q = b*b + a*c - doub4*d ;
  320.    r = (c - a*b)*c + asq*d ;
  321.    y = cubic(p,q,r) ;
  322.  
  323.    bmy = b - y ;
  324.    y4 = y*doub4 ;
  325.    d4 = d*doub4 ;
  326.    bmysq = bmy*bmy ;
  327.    gdis = asq - y4 ;
  328.    hdis = bmysq - d4 ;
  329.    if ((gdis < nought) || (hdis < nought)) return(0);
  330.    else
  331.    {
  332.       g1 = a*inv2 ;
  333.       h1 = bmy*inv2 ;
  334.       gerr = asq + y4 ;
  335.       herr = hdis ;
  336.       if (d > nought) herr = bmysq + d4 ;
  337.       if ((y < nought) || (herr*gdis > gerr*hdis))
  338.       {
  339.          gdisrt = sqrt(gdis) ;
  340.          g2 = gdisrt*inv2 ;
  341.          if (gdisrt != nought) h2 = (a*h1 - c)/gdisrt ;
  342.             else h2 = nought;
  343.       }
  344.       else
  345.       {
  346.          hdisrt = sqrt(hdis) ;
  347.          h2 = hdisrt*inv2 ;
  348.          if (hdisrt != nought) g2 = (a*h1 - c)/hdisrt ;
  349.             else g2 = nought;
  350.       }
  351. /* 
  352.      note that in the following, the tests ensure non-zero 
  353.      denominators -  
  354. */
  355.       h = h1 - h2 ;
  356.       hh = h1 + h2 ;
  357.       hmax = hh ;
  358.       if (hmax < nought) hmax =  -hmax ;
  359.       if (hmax < h) hmax = h ;
  360.       if (hmax <  -h) hmax =  -h ;
  361.       if ((h1 > nought)&&(h2 > nought)) h = d/hh ;
  362.       if ((h1 < nought)&&(h2 < nought)) h = d/hh ;
  363.       if ((h1 > nought)&&(h2 < nought)) hh = d/h ;
  364.       if ((h1 < nought)&&(h2 > nought)) hh = d/h ;
  365.       if (h > hmax) h = hmax ;
  366.       if (h <  -hmax) h =  -hmax ;
  367.       if (hh > hmax) hh = hmax ;
  368.       if (hh <  -hmax) hh =  -hmax ;
  369.  
  370.       g = g1 - g2 ;
  371.       gg = g1 + g2 ;
  372.       gmax = gg ;
  373.       if (gmax < nought) gmax =  -gmax ;
  374.       if (gmax < g) gmax = g ;
  375.       if (gmax <  -g) gmax =  -g ;
  376.       if ((g1 > nought)&&(g2 > nought)) g = y/gg ;
  377.       if ((g1 < nought)&&(g2 < nought)) g = y/gg ;
  378.       if ((g1 > nought)&&(g2 < nought)) gg = y/g ;
  379.       if ((g1 < nought)&&(g2 > nought)) gg = y/g ;
  380.       if (g > gmax) g = gmax ;
  381.       if (g <  -gmax) g =  -gmax ;
  382.       if (gg > gmax) gg = gmax ;
  383.       if (gg <  -gmax) gg =  -gmax ;
  384.  
  385.       n1 = qudrtc(gg,hh,v1, gg*gg - doub4*hh) ;
  386.       n2 = qudrtc(g,h,v2, g*g - doub4*h) ;
  387.       nquar = n1+n2 ;
  388.       rts[0] = v1[0] ;
  389.       rts[1] = v1[1] ;
  390.       rts[n1+0] = v2[0] ;
  391.       rts[n1+1] = v2[1] ;
  392.  
  393.       return(nquar);
  394.    }
  395. } /* neumark */
  396. /****************************************************/
  397.  
  398. void errors(a,b,c,d,rts,rterr,nrts)
  399. double a,b,c,d,rts[4],rterr[4];
  400. int nrts;
  401. /*
  402.    find the errors
  403.  
  404.    called by quartic.
  405. */
  406. {
  407.    int k;
  408.    double deriv,test;
  409.    double fabs(),sqrt(),curoot();
  410.  
  411.    if (nrts > 0)
  412.    {
  413.       for (  k = 0 ; k < nrts ; ++ k )
  414.       {
  415.          test = (((rts[k]+a)*rts[k]+b)*rts[k]+c)*rts[k]+d ;
  416.          if (test == nought) rterr[k] = nought;
  417.          else
  418.          {
  419.             deriv =
  420.                ((doub4*rts[k]+doub3*a)*rts[k]+doub2*b)*rts[k]+c ;
  421.             if (deriv != nought)
  422.                rterr[k] = fabs(test/deriv);
  423.             else
  424.             {
  425.                deriv = (doub12*rts[k]+doub6*a)*rts[k]+doub2*b ;
  426.                if (deriv != nought)
  427.                    rterr[k] = sqrt(fabs(test/deriv)) ;
  428.                else
  429.                {
  430.                   deriv = doub24*rts[k]+doub6*a ;
  431.                   if (deriv != nought)
  432.                      rterr[k] = curoot(fabs(test/deriv));
  433.                   else
  434.                      rterr[k] = sqrt(sqrt(fabs(test)/doub24));
  435.                }
  436.             }
  437.          }
  438.          fprintf(stderr,"errorsa  %d %9g %9g\n",
  439.                k,rts[k],rterr[k]);
  440.       }
  441.    }
  442.    else fprintf(stderr,"errors ans: none\n");
  443. } /* errors */
  444. /**********************************************/
  445.  
  446. int qudrtc(b,c,rts,dis)
  447.    double b,c,rts[4],dis ;
  448. /* 
  449.      solve the quadratic equation - 
  450.  
  451.          x**2+b*x+c = 0 
  452.  
  453.      called by  quartic, ferrari, neumark, ellcut 
  454. */
  455. {
  456.    int nquad;
  457.    double rtdis ;
  458.  
  459.    if (dis >= nought)
  460.    {
  461.       nquad = 2 ;
  462.       rtdis = sqrt(dis) ;
  463.       if (b > nought) rts[0] = ( -b - rtdis)*inv2 ;
  464.          else rts[0] = ( -b + rtdis)*inv2 ;
  465.       if (rts[0] == nought) rts[1] =  -b ;
  466.       else rts[1] = c/rts[0] ;
  467.    }
  468.    else
  469.    {
  470.       nquad = 0;
  471.       rts[0] = nought ;
  472.       rts[1] = nought ;
  473.    }
  474.    return(nquad);
  475. } /* qudrtc */
  476. /**************************************************/
  477.  
  478. double cubic(p,q,r)
  479. double p,q,r;
  480. /* 
  481.      find the lowest real root of the cubic - 
  482.        x**3 + p*x**2 + q*x + r = 0 
  483.  
  484.    input parameters - 
  485.      p,q,r - coeffs of cubic equation. 
  486.  
  487.    output- 
  488.      cubic - a real root. 
  489.  
  490.    global constants -
  491.      rt3 - sqrt(3) 
  492.      inv3 - 1/3 
  493.      doubmax - square root of largest number held by machine 
  494.  
  495.      method - 
  496.      see D.E. Littlewood, "A University Algebra" pp.173 - 6 
  497.  
  498.      Charles Prineas   April 1981 
  499.  
  500.      called by  neumark.
  501.      calls  acos3 
  502. */
  503. {
  504.    int nrts;
  505.    double po3,po3sq,qo3;
  506.    double uo3,u2o3,uo3sq4,uo3cu4 ;
  507.    double v,vsq,wsq ;
  508.    double m,mcube,n;
  509.    double muo3,s,scube,t,cosk,sinsqk ;
  510.    double root;
  511.    double curoot();
  512.    double acos3();
  513.    double sqrt(),fabs();
  514.  
  515.    m = nought;
  516.    nrts =0;
  517.    if ((p > doubmax) || (p <  -doubmax)) root = -p;
  518.    else
  519.    if ((q > doubmax) || (q <  -doubmax))
  520.    {
  521.       if (q > nought) root = -r/q;
  522.       else root = -sqrt(-q);
  523.    }
  524.    else
  525.    if ((r > doubmax)|| (r <  -doubmax)) root =  -curoot(r) ;
  526.    else
  527.    {
  528.       po3 = p*inv3 ;
  529.       po3sq = po3*po3 ;
  530.       if (po3sq > doubmax) root =  -p ;
  531.       else
  532.       {
  533.          v = r + po3*(po3sq + po3sq - q) ;
  534.          if ((v > doubmax) || (v < -doubmax)) root = -p ;
  535.          else
  536.          {
  537.             vsq = v*v ;
  538.             qo3 = q*inv3 ;
  539.             uo3 = qo3 - po3sq ;
  540.             u2o3 = uo3 + uo3 ;
  541.             if ((u2o3 > doubmax) || (u2o3 < -doubmax))
  542.             {
  543.                if (p == nought)
  544.                {
  545.                   if (q > nought) root =  -r/q ;
  546.                   else root =  -sqrt(-q) ;
  547.                }
  548.                else root =  -q/p ;
  549.             }
  550.             uo3sq4 = u2o3*u2o3 ;
  551.             if (uo3sq4 > doubmax)
  552.             {
  553.                if (p == nought)
  554.                {
  555.                   if (q > nought) root = -r/q ;
  556.                   else root = -sqrt(fabs(q)) ;
  557.                }
  558.                else root = -q/p ;
  559.             }
  560.             uo3cu4 = uo3sq4*uo3 ;
  561.             wsq = uo3cu4 + vsq ;
  562.             if (wsq >= nought)
  563.             {
  564. /* 
  565.      cubic has one real root 
  566. */
  567.                nrts = 1;
  568.                if (v <= nought) mcube = ( -v + sqrt(wsq))*inv2 ;
  569.                if (v  > nought) mcube = ( -v - sqrt(wsq))*inv2 ;
  570.                m = curoot(mcube) ;
  571.                if (m != nought) n = -uo3/m ;
  572.                   else n = nought;
  573.                root = m + n - po3 ;
  574.             }
  575.             else
  576.             {
  577.                nrts = 3;
  578. /* 
  579.      cubic has three real roots 
  580. */
  581.                if (uo3 < nought)
  582.                {
  583.                   muo3 = -uo3;
  584.                   s = sqrt(muo3) ;
  585.                   scube = s*muo3;
  586.                   t =  -v/(scube+scube) ;
  587.                   cosk = acos3(t) ;
  588.                   if (po3 < nought)
  589.                      root = (s+s)*cosk - po3;
  590.                   else
  591.                   {
  592.                      sinsqk = doub1 - cosk*cosk ;
  593.                      if (sinsqk < nought) sinsqk = nought ;
  594.                      root = s*( -cosk - rt3*sqrt(sinsqk)) - po3 ;
  595.                   }
  596.                }
  597.                else
  598. /* 
  599.      cubic has multiple root -  
  600. */
  601.                root = curoot(v) - po3 ;
  602.             }
  603.          }
  604.       }
  605.    }
  606.    fprintf(stderr,"cubic %g %g %g %d %g\n",p,q,r,nrts,root);
  607.    return(root);
  608. } /* cubic */
  609. /***************************************/
  610.  
  611. double acos3(x)
  612.    double x ;
  613. /* 
  614.      find cos(acos(x)/3) 
  615.     
  616.      Don Herbison-Evans   16/7/81 
  617.  
  618.      called by cubic . 
  619. */
  620. {
  621.    double value;
  622.    double acos(),cos();
  623.  
  624.    value = cos(acos(x)*inv3);
  625.    return(value);
  626. } /* acos3 */
  627. /***************************************/
  628.  
  629. double curoot(x)
  630.    double x ;
  631. /* 
  632.      find cube root of x.
  633.  
  634.      Don Herbison-Evans   30/1/89 
  635.  
  636.      called by cubic . 
  637. */
  638. {
  639.    double exp(),log();
  640.    double value;
  641.    double absx;
  642.    int neg;
  643.  
  644.    neg = 0;
  645.    absx = x;
  646.    if (x < nought)
  647.    {
  648.       absx = -x;
  649.       neg = 1;
  650.    }
  651.    value = exp( log(absx)*inv3 );
  652.    if (neg == 1) value = -value;
  653.    return(value);
  654. } /* curoot */
  655. /****************************************************/
  656.  
  657. int simple(a,b,c,d,rts)
  658.    double a,b,c,d,rts[4];
  659. /* 
  660.      solve the quartic equation - 
  661.  
  662.    x**4 + a*x**3 + b*x**2 + c*x + d = 0 
  663.  
  664.    called by quartic
  665.    calls     cubic, qudrtc.
  666.  
  667.      input - 
  668.    a,b,c,e - coeffs of equation. 
  669.  
  670.      output - 
  671.    nquar - number of real roots. 
  672.    rts - array of root values. 
  673.  
  674.      method :  unstabilized Ferrari-Lagrange
  675.      Abramowitz,M. & Stegun I.A.
  676.      Handbook of Mathematical Functions
  677.      Dover 1972 (ninth printing), pp. 17-18
  678.  
  679.      calls  cubic, qudrtc 
  680. */
  681. {
  682.    double cubic();
  683.    int qudrtc();
  684.  
  685.    int nquar,n1,n2 ;
  686.    double asq,y;
  687.    double v1[4],v2[4] ;
  688.    double p,q,r ;
  689.    double e,f,esq,fsq ;
  690.    double g,gg,h,hh;
  691.  
  692.    fprintf(stderr,"\nsimple %g %g %g %g\n",a,b,c,d);
  693.    asq = a*a;
  694.  
  695.    p = -b ;
  696.    q = a*c-doub4*d ;
  697.    r = -asq*d - c*c + doub4*b*d ;
  698.    y = cubic(p,q,r) ;
  699.  
  700.    esq = inv4*asq - b + y;
  701.    fsq = inv4*y*y - d;
  702.    if (esq < nought) return(0);
  703.    else
  704.    if (fsq < nought) return(0);
  705.    else
  706.    {
  707.       e = sqrt(esq);
  708.       f = sqrt(fsq);
  709.       g = inv2*a - e;
  710.       h = inv2*y - f;
  711.       gg = inv2*a + e;
  712.       hh = inv2*y + f;
  713.       n1 = qudrtc(gg,hh,v1, gg*gg - doub4*hh) ;
  714.       n2 = qudrtc(g,h,v2, g*g - doub4*h) ;
  715.       nquar = n1+n2 ;
  716.       rts[0] = v1[0] ;
  717.       rts[1] = v1[1] ;
  718.       rts[n1+0] = v2[0] ;
  719.       rts[n1+1] = v2[1] ;
  720.       return(nquar);
  721.    }
  722. } /* simple */
  723. /*****************************************/
  724.  
  725. int descartes(a,b,c,d,rts)
  726. double a,b,c,d,rts[4];
  727. /*
  728.    Solve quartic equation using
  729.    Descartes-Euler-Cardano algorithm
  730.  
  731.    Strong, T. "Elemementary and Higher Algebra"
  732.       Pratt and Oakley, p. 469 (1859)
  733.  
  734.      29 Jun 1994  Don Herbison-Evans
  735. */
  736. {
  737.    int qudrtc();
  738.    double cubic();
  739.  
  740.    int nrts;
  741.    int r1,r2;
  742.    double v1[4],v2[4];
  743.    double y;
  744.    double p,q,r;
  745.    double A,B,C;
  746.    double m,n1,n2;
  747.    double d3o8,d3o256;
  748.    double inv8,inv16;
  749.    double asq;
  750.    double Binvm;
  751.  
  752.    d3o8 = (double)3/(double)8;
  753.    inv8 = doub1/(double)8;
  754.    inv16 = doub1/(double)16;
  755.    d3o256 = (double)3/(double)256;
  756.  
  757.    fprintf(stderr,"\nDescartes %f %f %f %f\n",a,b,c,d);
  758.    asq = a*a;
  759.  
  760.    A = b - asq*d3o8;
  761.    B = c + a*(asq*inv8 - b*inv2);
  762.    C = d + asq*(b*inv16 - asq*d3o256) - a*c*inv4;
  763.  
  764.    p = doub2*A;
  765.    q = A*A - doub4*C;
  766.    r = -B*B;
  767.  
  768. /***
  769.    inv64 = doub1/(double)64;
  770.    p = doub2*b - doub3*a*a*inv4 ;
  771.    q = b*b - a*a*b - doub4*d + doub3*a*a*a*a*inv16 + a*c;
  772.    r = -c*c - a*a*a*a*a*a*inv64 - a*a*b*b*inv4
  773.        -a*a*a*c*inv4 + a*b*c + a*a*a*a*b*inv8;
  774. ***/
  775.  
  776.    y = cubic(p,q,r) ;
  777.    if (y <= nought) 
  778.       nrts = 0;
  779.    else
  780.    {
  781.       m = sqrt(y);
  782.       Binvm = B/m;
  783.       n1 = (y + A + Binvm)*inv2;
  784.       n2 = (y + A - Binvm)*inv2;
  785.       r1 = qudrtc(-m, n1, v1, y-doub4*n1);
  786.       r2 = qudrtc( m, n2, v2, y-doub4*n2);
  787.       rts[0] = v1[0]-a*inv4;
  788.       rts[1] = v1[1]-a*inv4;
  789.       rts[r1] = v2[0]-a*inv4;
  790.       rts[r1+1] = v2[1]-a*inv4;
  791.       nrts = r1+r2;
  792.    } 
  793.    return(nrts);
  794. } /* descartes */
  795. /****************************************************/
  796.