home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / NAMULL.C < prev    next >
Encoding:
Text File  |  1984-08-30  |  3.6 KB  |  153 lines
  1.    namull(root,n,eps,max)
  2.  
  3.       /* subroutine finds n complex roots of the equation f(x) = 0 */
  4.       /* using Muller's method. f(x) may be transcendental and     */
  5.       /* complex.                                                  */
  6.  
  7.       int max,*n;
  8.       float eps,root[];
  9.  
  10.     {
  11.       int i1=1,i2=2,i3=3,i4=4,inf,j,j1,j2,k,k1,k2,m,mm;
  12.       float alam[2],b1[2],b2[2],b3[2],d[2],delta[2],f1[2],f2[2],f3[2];
  13.       float fprt[2],g[2],h[2],temp[2],rt[2],xn[2],s;
  14.       static float one[2]={1.,0.},two[2]={-2.,0.},four[2]={4.,0.};
  15.       extern double fabs(),sqrt();
  16.  
  17.       k = 0;
  18. reit: k = k + 1;
  19.       mm = 0;
  20.       k2 = k + k;
  21.       k1 = k2 - 1;
  22.       if(root[k1-1] != 0.0 || root[k2-1] != 0.0)
  23.       {
  24.       rt[0] = 0.9 * root[k1-1];
  25.       rt[1] = 0.9 * root[k2-1];
  26.       m = 1;
  27.       }
  28.        else
  29.        {
  30.         rt[0] = -1.0;
  31.         rt[1] = 0.0;
  32.         h[0] = -1.0;
  33.         h[1] = 0.0;
  34.         alam[0] = -0.5;
  35.         alam[1] = 0.0;
  36.         m = 3;
  37.        }
  38.         inf = 1;
  39.       while(inf == 1)
  40.       {
  41.        value(rt,fprt);
  42.        s = sqrt(fprt[0]*fprt[0] + fprt[1]*fprt[1]);
  43.        if(s < eps)
  44.        {
  45.         root[k1-1] = xn[0];
  46.         root[k2-1] = xn[1];
  47.         if(k != *n) goto reit;
  48.         return;
  49.        }
  50.        mm = mm + 1;
  51.        if(mm >= max)
  52.        {
  53.         *n = k - 1;
  54.         return;
  55.        }
  56.        if(k != 1)
  57.        {
  58.         for(j = 1; j <= k-1; j++)
  59.         {
  60.          j2 = j + j;
  61.          j1 = j2 - 1;
  62.          temp[0] = rt[0] - root[j1-1];
  63.          temp[1] = rt[1] - root[j2-1];
  64.          cxarth(fprt,temp,fprt,i4);
  65.         }
  66.        }
  67.        switch(m){
  68.  case 1:
  69.        f1[0] = fprt[0];
  70.        f1[1] = fprt[1];
  71.        rt[0] = 1.1 * root[k1-1];
  72.        rt[1] = 1.1 * root[k2-1];
  73.        m = 2;
  74.        break;
  75.  case 2:
  76.        f2[0] = fprt[0];
  77.        f2[1] = fprt[1];
  78.        h[0] = root[k1-1] - rt[0];
  79.        h[1] = root[k2-1] - rt[1];
  80.        alam[0] = .2 * root[k1-1];
  81.        alam[1] = .2 * root[k2-1];
  82.        rt[0] = root[k1-1];
  83.        rt[1] = root[k2-1];
  84.        m = 5;
  85.        break;
  86.  case 3:
  87.        f1[0] = fprt[0];
  88.        f1[1] = fprt[1];
  89.        rt[0] = 1.0;
  90.        rt[1] = 0.0;
  91.        m = 4;
  92.        break;
  93.  case 4:
  94.        f2[0] = fprt[0];
  95.        f2[1] = fprt[1];
  96.        rt[0] = 0.0;
  97.        rt[1] = 0.0;
  98.        m = 5;
  99.        break;
  100.  case 6:
  101.        f1[0] = f2[0];
  102.        f1[1] = f2[1];
  103.        f2[0] = f3[0];
  104.        f2[1] = f3[1];
  105.  case 5:
  106.        f3[0] = fprt[0];
  107.        f3[1] = fprt[1];
  108.        cxarth(one,alam,delta,i1);
  109.        cxarth(f1,alam,b1,i3);
  110.        cxarth(b1,alam,b1,i3);
  111.        cxarth (delta,delta,b2,i3);
  112.        cxarth(b2,f2,b2,i3);
  113.        cxarth(alam,delta,b3,i1);
  114.        cxarth(b3,f3,b3,i3);
  115.        cxarth(b1,b2,g,i2);
  116.        cxarth(g,b3,g,i1);
  117.        cxarth(f1,alam,b1,i3);
  118.        cxarth(f2,delta,b2,i3);
  119.        cxarth(b1,b2,b2,i2);
  120.        cxarth(b2,f3,b2,i1);
  121.        cxarth(b2,alam,b2,i3);
  122.        cxarth(b2,delta,b2,i3);
  123.        cxarth(b2,f3,b2,i3);
  124.        cxarth(b2,four,b2,i3);
  125.        cxarth(g,g,b1,i3);
  126.        cxarth(b1,b2,d,i2);
  127.        cxsqrt(d,d);
  128.        cxarth(g,d,b1,i1);
  129.        cxarth(g,d,b2,i2);
  130.        d[0] = (b1[0] - b1[1])*(b1[0] + b1[1]);
  131.        d[1] = (b2[0] - b2[1])*(b2[0] + b2[1]);
  132.        if(d[0] >= d[1])
  133.        {
  134.        d[0] = b1[0];
  135.        d[1] = b1[1];
  136.        }
  137.         else
  138.         {
  139.          d[0] = b2[0];
  140.          d[1] = b2[1];
  141.         }
  142.        cxarth(f3,delta,b1,i3);
  143.        cxarth(b1,two,b1,i3);
  144.        cxarth(b1,d,alam,i4);
  145.        cxarth(h,alam,h,i3);
  146.        cxarth(rt,h,xn,i1);
  147.        rt[0] = xn[0];
  148.        rt[1] = xn[1];
  149.        m = 6;
  150.        }
  151.       }
  152.      }
  153.