home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / BIQUAD.C < prev    next >
Encoding:
Text File  |  1984-07-25  |  2.1 KB  |  87 lines
  1.    biquad(a,r1,r2,r3,r4,tol,numrt)
  2.  
  3.       /*solves fourth order equations of the form         */
  4.       /*a(1)x**4 + a(2)x**3 + a(3)x**2 + a(4)x + a(5) = 0 */
  5.       /*for real and/or imaginary roots.                  */
  6.  
  7.      int *numrt;
  8.      float a[],r1[],r2[],r3[],r4[],tol;
  9.  
  10.      {
  11.       int i,n;
  12.       float xcoefs[4],ycoefs[4],yroots[6],zero;
  13.       float aa,b,c,d,e,f,y;
  14.       extern double fabs(),sqrt();
  15.  
  16.       zero = tol/10.;
  17.       if ((fabs(a[0]) - zero) <= 0.)
  18.       {
  19.       for(i = 0; i <= 3; i++)
  20.        xcoefs[i] = a[i + 1];
  21.       cubic (xcoefs,r1,r2,r3,tol,&n);
  22.       *numrt = n;
  23.       return;
  24.       }
  25.       *numrt = 4;
  26.       for(i = 0; i <= 3; i++)
  27.        xcoefs[i] = a[i + 1]/a[0];
  28.  
  29.       ycoefs[0] = 1.0;
  30.       ycoefs[1] = -xcoefs[1];
  31.       ycoefs[2] = xcoefs[0]*xcoefs[2] - 4.0*xcoefs[3];
  32.       ycoefs[3] = 4.0*xcoefs[1]*xcoefs[3]-xcoefs[0]*xcoefs[0]*xcoefs[3];
  33.       ycoefs[3] = ycoefs[3] - xcoefs[2]*xcoefs[2];
  34.  
  35.       cubic(ycoefs,r1,r2,r3,tol,&n);
  36.  
  37.       yroots[0] = r1[0];
  38.       yroots[1] = r1[1];
  39.       yroots[2] = r2[0];
  40.       yroots[3] = r2[1];
  41.       yroots[4] = r3[0];
  42.       yroots[5] = r3[1];
  43.  
  44.       y = 0.0;
  45.       for(i = 1; i <= 5; i+=2)
  46.       {
  47.        if((fabs(yroots[i]) - zero) <= 0.)
  48.        {
  49.         if((y - yroots[i-1]) <0.)
  50.         y = yroots[i - 1];
  51.        }
  52.       }
  53.       aa = xcoefs[0] * xcoefs[0]*0.25 - xcoefs[1] + y;
  54.       c = -xcoefs[3] + y*y*0.25;
  55.       if((fabs(aa) - zero) <= 0.)
  56.       {
  57.        e = 0.0;
  58.        if((fabs(c) - zero) <= 0.)
  59.        {
  60.         f = 0.0;
  61.         goto r10;
  62.        }
  63.        if(c < 0.) goto r20;
  64.         f = sqrt(c);
  65.         goto r10;
  66.       }
  67.       if(aa < 0.) goto r20;
  68.       e = sqrt(aa);
  69.       b = -xcoefs[2] + 0.5 * xcoefs[0] * y;
  70.       f = b/e*0.5;
  71.  
  72. r10:  ycoefs[0] = 1.0;
  73.       ycoefs[1] = xcoefs[0]*0.5 - e;
  74.       ycoefs[2] = y*0.5 - f;
  75.       quadrt(ycoefs,r1,r2,tol,&n);
  76.       ycoefs[0] = 1.0;
  77.       ycoefs[1] = xcoefs[0]*0.5 + e;
  78.       ycoefs[2] = y*0.5 + f;
  79.       quadrt(ycoefs,r3,r4,tol,&n);
  80.       return;
  81. r20:  *numrt = -*numrt;
  82.       return;
  83.  
  84.  
  85.  
  86.     }
  87.