home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / CUBIC.C < prev    next >
Encoding:
Text File  |  1984-09-01  |  2.0 KB  |  89 lines
  1.     cubic(a,rt1,rt2,rt3,tol,numrt)
  2.  
  3.       /*this function finds the roots (real or   */
  4.       /*real and complex) to any cubic equation. */
  5.  
  6.       int *numrt;
  7.       float a[],rt1[],rt2[],rt3[],tol;
  8.  
  9.     {
  10.  
  11.       int it,n,nda,i;
  12.       float ang,b[2],c[4],work[3],x,y,zz,zero;
  13.       extern double fabs(),sqrt(),sin(),cos(),acos(),exp(),log(),pow();
  14.  
  15.       zero = tol/10.;
  16.  
  17.       if((fabs(a[0]) - zero) <= 0)
  18.       {
  19.        for(i = 0;i <= 2; i++)
  20.         work[i] = a[i+1];
  21.       quadrt(work, rt1,rt2,tol,&n);
  22.       *numrt = n;
  23.       return;
  24.       }
  25.       *numrt = 3;
  26.       it = 0;
  27.       for(i = 1; i <= 3; i++)
  28.        a[i] = a[i]/a[0];
  29.       a[0] = 1.0;
  30.  
  31.       if(fabs(a[1]) > zero)
  32.       {
  33.        nda = 3;
  34.        b[0] = 1.0;
  35.        b[1] = -(a[1]*.333333);
  36.        lincng(a,nda,b,c);
  37.        it = 1;
  38.       }
  39.        else
  40.        {
  41.         for(i = 0; i <= 3; i++)
  42.          c[i] = a[i];
  43.        }
  44.       x = c[3] * c[3] * .25 + c[2]*c[2]*c[2] * .037037;
  45.       if(x >= 0.0)
  46.       {
  47.        x = sqrt(x);
  48.        y = -(c[3]*.5);
  49.        i = 0;
  50.        rt1[i] = y + x;
  51. l5:     n = 0;
  52.         if(rt1[0] <  0.0) n = 1;
  53.         if(fabs(rt1[0]) > zero)
  54.         {
  55.          rt1[i] = exp((log(fabs(rt1[i])))*.333333);
  56.          if(n == 1) rt1[i] = -rt1[i];
  57.         }
  58.       if(i == 1) goto l8;
  59.       i = 1;
  60.       rt1[i] = y - x;
  61.       goto l5;
  62.  
  63. l8:   rt2[1] = ((rt1[0] - rt1[1])*.5)*1.732051;
  64.       rt1[0] = rt1[0] + rt1[1];
  65.       rt1[1] = 0.0;
  66.       rt2[0] = - rt1[0] * .5;
  67.       rt3[0] = rt2[0];
  68.       rt3[1] = -rt2[1];
  69.       }
  70.        else
  71.        {
  72.         zz = fabs(c[2]);
  73.         x = -(c[3]*.5)/sqrt(pow(zz,3.)*.037037);
  74.         ang= acos(x);
  75.         y = 2.0 * (sqrt(zz*.333333));
  76.         ang = ang*.333333;
  77.         rt1[0] = y * cos(ang);
  78.         rt2[0] = y * cos(ang + 2.09440);
  79.         rt3[0] = y * cos(ang + 4.18879);
  80.         rt1[1] = 0.0;
  81.         rt2[1] = 0.0;
  82.         rt3[1] = 0.0;
  83.        }
  84.       if(it == 0) return;
  85.       rt1[0] = rt1[0] + b[1];
  86.       rt2[0] = rt2[0] + b[1];
  87.       rt3[0] = rt3[0] + b[1];
  88.     }
  89.