/* polrt.c * * Find roots of a polynomial * * * * SYNOPSIS: * * typedef struct * { * double r; * double i; * }cmplx; * * double xcof[], cof[]; * int m; * cmplx root[]; * * polrt( xcof, cof, m, root ) * * * * DESCRIPTION: * * Iterative determination of the roots of a polynomial of * degree m whose coefficient vector is xcof[]. The * coefficients are arranged in ascending order; i.e., the * coefficient of x**m is xcof[m]. * * The array cof[] is working storage the same size as xcof[]. * root[] is the output array containing the complex roots. * * * ACCURACY: * * Termination depends on evaluation of the polynomial at * the trial values of the roots. The values of multiple roots * or of roots that are nearly equal may have poor relative * accuracy after the first root in the neighborhood has been * found. * */ /* polrt */ /* Complex roots of real polynomial */ /* number of coefficients is m + 1 ( i.e., m is degree of polynomial) */ #include "mconf.h" /* typedef struct { double r; double i; }cmplx; */ #ifdef ANSIPROT extern double fabs ( double ); #else double fabs(); #endif int polrt( xcof, cof, m, root ) double xcof[], cof[]; int m; cmplx root[]; { register double *p, *q; int i, j, nsav, n, n1, n2, nroot, iter, retry; int final; double mag, cofj; cmplx x0, x, xsav, dx, t, t1, u, ud; final = 0; n = m; if( n <= 0 ) return(1); if( n > 36 ) return(2); if( xcof[m] == 0.0 ) return(4); n1 = n; n2 = n; nroot = 0; nsav = n; q = &xcof[0]; p = &cof[n]; for( j=0; j<=nsav; j++ ) *p-- = *q++; /* cof[ n-j ] = xcof[j];*/ xsav.r = 0.0; xsav.i = 0.0; nxtrut: x0.r = 0.00500101; x0.i = 0.01000101; retry = 0; tryagn: retry += 1; x.r = x0.r; x0.r = -10.0 * x0.i; x0.i = -10.0 * x.r; x.r = x0.r; x.i = x0.i; finitr: iter = 0; while( iter < 500 ) { u.r = cof[n]; if( u.r == 0.0 ) { /* this root is zero */ x.r = 0; n1 -= 1; n2 -= 1; goto zerrut; } u.i = 0; ud.r = 0; ud.i = 0; t.r = 1.0; t.i = 0; p = &cof[n-1]; for( i=0; i= 1.0e-5 ) { cofj = x.r + x.r; mag = x.r * x.r + x.i * x.i; n -= 2; } else { /* root is real */ zerrut: x.i = 0; cofj = x.r; mag = 0; n -= 1; } /* divide working polynomial cof(z) by z - x */ p = &cof[1]; *p += cofj * *(p-1); for( j=1; j 0 ) goto nxtrut; return(0); }