home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / cephes22 / polyn / polrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.4 KB  |  219 lines
  1. /*                            polrt.c
  2.  *
  3.  *    Find roots of a polynomial
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * typedef struct
  10.  *    {
  11.  *    double r;
  12.  *    double i;
  13.  *    }cmplx;
  14.  *
  15.  * double xcof[], cof[];
  16.  * int m;
  17.  * cmplx root[];
  18.  *
  19.  * polrt( xcof, cof, m, root )
  20.  *
  21.  *
  22.  *
  23.  * DESCRIPTION:
  24.  *
  25.  * Iterative determination of the roots of a polynomial of
  26.  * degree m whose coefficient vector is xcof[].  The
  27.  * coefficients are arranged in ascending order; i.e., the
  28.  * coefficient of x**m is xcof[m].
  29.  *
  30.  * The array cof[] is working storage the same size as xcof[].
  31.  * root[] is the output array containing the complex roots.
  32.  *
  33.  *
  34.  * ACCURACY:
  35.  *
  36.  * Termination depends on evaluation of the polynomial at
  37.  * the trial values of the roots.  The values of multiple roots
  38.  * or of roots that are nearly equal may have poor relative
  39.  * accuracy after the first root in the neighborhood has been
  40.  * found.
  41.  *
  42.  */
  43.  
  44. /*                            polrt    */
  45. /* Complex roots of real polynomial */
  46. /* number of coefficients is m + 1 ( i.e., m is degree of polynomial) */
  47.  
  48. typedef struct
  49.     {
  50.     double r;
  51.     double i;
  52.     }cmplx;
  53.  
  54. int polrt( xcof, cof, m, root )
  55. double xcof[], cof[];
  56. int m;
  57. cmplx root[];
  58. {
  59. register double *p, *q;
  60. int i, j, nsav, n, n1, n2, nroot, iter, retry;
  61. int final;
  62. double mag, cofj;
  63. cmplx x0, x, xsav, dx, t, t1, u, ud;
  64. double abs();
  65.  
  66. final = 0;
  67. n = m;
  68. if( n <= 0 )
  69.     return(1);
  70. if( n > 36 )
  71.     return(2);
  72. if( xcof[m] == 0.0 )
  73.     return(4);
  74.  
  75. n1 = n;
  76. n2 = n;
  77. nroot = 0;
  78. nsav = n;
  79. q = &xcof[0];
  80. p = &cof[n];
  81. for( j=0; j<=nsav; j++ )
  82.     *p-- = *q++;    /*    cof[ n-j ] = xcof[j];*/
  83.  
  84. nxtrut:
  85. x0.r = 0.00500101;
  86. x0.i = 0.01000101;
  87. retry = 0;
  88.  
  89. tryagn:
  90. retry += 1;
  91. x.r = x0.r;
  92.  
  93. x0.r = -10.0 * x0.i;
  94. x0.i = -10.0 * x.r;
  95.  
  96. x.r = x0.r;
  97. x.i = x0.i;
  98.  
  99. finitr:
  100. iter = 0;
  101.  
  102. while( iter < 500 )
  103. {
  104. u.r = cof[n];
  105. if( u.r == 0.0 )
  106.     {        /* this root is zero */
  107.     x.r = 0;
  108.     n1 -= 1;
  109.     n2 -= 1;
  110.     goto zerrut;
  111.     }
  112. u.i = 0;
  113. ud.r = 0;
  114. ud.i = 0;
  115. t.r = 1.0;
  116. t.i = 0;
  117. p = &cof[n-1];
  118. for( i=0; i<n; i++ )
  119.     {
  120.     t1.r = x.r * t.r  -  x.i * t.i;
  121.     t1.i = x.r * t.i  +  x.i * t.r;
  122.     cofj = *p--;        /* evaluate polynomial */
  123.     u.r += cofj * t1.r;
  124.     u.i += cofj * t1.i;
  125.     cofj = cofj * (i+1);    /* derivative */
  126.     ud.r += cofj * t.r;
  127.     ud.i -= cofj * t.i;
  128.     t.r = t1.r;
  129.     t.i = t1.i;
  130.     }
  131.  
  132. mag = ud.r * ud.r  +  ud.i * ud.i;
  133. if( mag == 0.0 )
  134.     {
  135.     if( !final )
  136.         goto tryagn;
  137.     x.r = xsav.r;
  138.     x.i = xsav.i;
  139.     goto findon;
  140.     }
  141. dx.r = (u.i * ud.i  -  u.r * ud.r)/mag;
  142. x.r += dx.r;
  143. dx.i = -(u.r * ud.i  +  u.i * ud.r)/mag;
  144. x.i += dx.i;
  145. if( (abs(dx.i) + abs(dx.r)) < 1.0e-6 )
  146.     goto lupdon;
  147. iter += 1;
  148. }    /* while iter < 500 */
  149.  
  150. if( final )
  151.     goto lupdon;
  152. if( retry < 5 )
  153.     goto tryagn;
  154. return(3);
  155.  
  156. lupdon:
  157. /* Swap original and reduced polynomials */
  158. q = &xcof[nsav];
  159. p = &cof[0];
  160. for( j=0; j<=n2; j++ )
  161.     {
  162.     cofj = *q;
  163.     *q-- = *p;
  164.     *p++ = cofj;
  165.     }
  166. i = n;
  167. n = n1;
  168. n1 = i;
  169.  
  170. if( !final )
  171.     {
  172.     final = 1;
  173.     if( abs(x.i/x.r) < 1.0e-4 )
  174.         x.i = 0.0;
  175.     xsav.r = x.r;
  176.     xsav.i = x.i;
  177.     goto finitr;    /* do final iteration on original polynomial */
  178.     }
  179.  
  180. findon:
  181. final = 0;
  182. if( abs(x.i/x.r) >= 1.0e-5 )
  183.     {
  184.     cofj = x.r + x.r;
  185.     mag = x.r * x.r  +  x.i * x.i;
  186.     n -= 2;
  187.     }
  188. else
  189.     {        /* root is real */
  190. zerrut:
  191.     x.i = 0;
  192.     cofj = x.r;
  193.     mag = 0;
  194.     n -= 1;
  195.     }
  196. /* divide working polynomial cof(z) by z - x */
  197. p = &cof[1];
  198. *p += cofj * *(p-1);
  199. for( j=1; j<n; j++ )
  200.     {
  201.     *(p+1) += cofj * *p  -  mag * *(p-1);
  202.     p++;
  203.     }
  204.  
  205. setrut:
  206. root[nroot].r = x.r;
  207. root[nroot].i = x.i;
  208. nroot += 1;
  209. if( mag != 0.0 )
  210.     {
  211.     x.i = -x.i;
  212.     mag = 0;
  213.     goto setrut;    /* fill in the complex conjugate root */
  214.     }
  215. if( n > 0 )
  216.     goto nxtrut;
  217. return(0);
  218. }
  219.