3 * Find roots of a polynomial
15 * double xcof[], cof[];
19 * polrt( xcof, cof, m, root )
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].
30 * The array cof[] is working storage the same size as xcof[].
31 * root[] is the output array containing the complex roots.
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
45 /* Complex roots of real polynomial */
46 /* number of coefficients is m + 1 ( i.e., m is degree of polynomial) */
57 extern double fabs ( double );
62 int polrt( xcof, cof, m, root )
67 register double *p, *q;
68 int i, j, nsav, n, n1, n2, nroot, iter, retry;
71 cmplx x0, x, xsav, dx, t, t1, u, ud;
88 for( j=0; j<=nsav; j++ )
89 *p-- = *q++; /* cof[ n-j ] = xcof[j];*/
115 { /* this root is zero */
129 t1.r = x.r * t.r - x.i * t.i;
130 t1.i = x.r * t.i + x.i * t.r;
131 cofj = *p--; /* evaluate polynomial */
134 cofj = cofj * (i+1); /* derivative */
141 mag = ud.r * ud.r + ud.i * ud.i;
150 dx.r = (u.i * ud.i - u.r * ud.r)/mag;
152 dx.i = -(u.r * ud.i + u.i * ud.r)/mag;
154 if( (fabs(dx.i) + fabs(dx.r)) < 1.0e-6 )
157 } /* while iter < 500 */
166 /* Swap original and reduced polynomials */
169 for( j=0; j<=n2; j++ )
182 if( fabs(x.i/x.r) < 1.0e-4 )
186 goto finitr; /* do final iteration on original polynomial */
191 if( fabs(x.i/x.r) >= 1.0e-5 )
194 mag = x.r * x.r + x.i * x.i;
205 /* divide working polynomial cof(z) by z - x */
210 *(p+1) += cofj * *p - mag * *(p-1);
222 goto setrut; /* fill in the complex conjugate root */