home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / NABAIR.C < prev    next >
Encoding:
Text File  |  1984-07-30  |  2.6 KB  |  121 lines
  1.    nabair(p1,q1,a,itag,n,eps,root,b,c)
  2.  
  3.       /*Bairstow's method for finding roots of nth degree polynomial    */
  4.       /*with real coefficients. polynomial is of the form               */
  5.       /*a(0)x**n + a(1)x**(n - 1) + a(2)x**(n-2) + ... + a(n-1)x + a(n).*/
  6.  
  7.       int *itag,*n;
  8.       float a[],b[],c[],eps,p1,q1,root[];
  9.  
  10.     {
  11.  
  12.       int i,iitag,j,l,m;
  13.       float cbarl,d,deltp,deltq,den,e,f,p,q,sum,sum1;
  14.       extern double fabs(),sqrt();
  15.  
  16.       iitag = *itag;
  17.       *itag = 0;
  18.       d = a[0];
  19.  
  20.       for(i = 0; i <= *n-1; i++)
  21.        a[i] = a[i + 1]/d;
  22.  
  23. strt:if(*n <= 1)
  24.       {
  25.        *itag = *itag + 1;
  26.        l = *itag + *itag - 1;
  27.        root[l] = 0.;
  28.        root[l-1] = -a[0];
  29.        *n = *itag;
  30.        *itag = 4;
  31.        return;
  32.       }
  33.       if(*n <= 2)
  34.       {
  35.        p = a[0];
  36.        q = a[1];
  37.        goto intr;
  38.       }
  39.       p = p1;
  40.       q = q1;
  41.       m = 1;
  42. reit: b[0] = a[0] - p;
  43.       b[1] = a[1] - p * b[0] - q;
  44.       l = *n - 2;
  45.       c[0] = b[0] - p;
  46.       c[1] = b[1] - p * c[0] - q;
  47.       if(l == 1) l = 2;                            /* watch this!!! */
  48.       for(j = 2; j <= l; j++)
  49.       {
  50.        b[j] = a[j] - p * b[j-1] - q * b[j-2];
  51.        c[j] = b[j] - p * c[j-1] - q * c[j-2];
  52.       }
  53.       l = *n - 2;
  54.       cbarl = c[l] - b[l];
  55.       den = -cbarl;
  56.       if(*n != 3) den = den * c[*n-4];
  57.       den = den + c[*n-3]*c[*n-3];
  58.       if(den == 0.0)
  59.       {
  60.        *n = *itag;
  61.        *itag = 1;
  62.        return;
  63.       }
  64.  
  65.       b[*n-1] = a[*n-1] - p * b[*n-2] - q * b[*n-3];
  66.       deltp = -b[*n-1];
  67.       if(*n != 3) deltp = deltp * c[*n-4];
  68.       deltp = (b[*n-2] * c[*n-3] + deltp)/den;
  69.       deltq = (b[*n-1]*c[*n-3] - b[*n-2]*cbarl)/den;
  70.       p = p + deltp;
  71.       q = q + deltq;
  72.       sum = fabs(deltp) + fabs(deltq);
  73.       if(m == 1) sum1 = sum;
  74.       if(m == 5 && sum > sum1)
  75.       {
  76.       *n = *itag;
  77.       *itag = 2;
  78.       return;
  79.       }
  80.       if(sum > eps)
  81.       {
  82.        if(m >= iitag)
  83.        {
  84.         *n = *itag;
  85.         *itag = 3;
  86.         return;
  87.        }
  88.         m ++;
  89.         goto reit;
  90.       }
  91. intr: d = -p * 0.5;
  92.       *itag = *itag + 2;
  93.       f = q - d*d;
  94.       e = sqrt(fabs(f));
  95.       l = *itag + *itag - 1;
  96.       if(f <= 0.0)
  97.       {
  98.       root[l] = 0.0;
  99.       root[l - 1] = d - e;
  100.       root[l-2] = 0.0;
  101.       root[l-3] = d + e;
  102.       }
  103.        else
  104.        {
  105.       root[l] = -e;
  106.       root[l - 1] = d;
  107.       root[l - 2] = e;
  108.       root[l - 3] = d;
  109.        }
  110.       *n = *n - 2;
  111.       if(*n <= 0)
  112.       {
  113.       *n = *itag;
  114.       *itag = 4;
  115.       return;
  116.       }
  117.       for(i = 0; i <= *n-1; i++)
  118.         a[i] = b[i];
  119.       goto strt;
  120.    }
  121.