/* 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 <math.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<n; i++ ) { t1.r = x.r * t.r - x.i * t.i; t1.i = x.r * t.i + x.i * t.r; cofj = *p--; /* evaluate polynomial */ u.r += cofj * t1.r; u.i += cofj * t1.i; cofj = cofj * (i+1); /* derivative */ ud.r += cofj * t.r; ud.i -= cofj * t.i; t.r = t1.r; t.i = t1.i; } mag = ud.r * ud.r + ud.i * ud.i; if( mag == 0.0 ) { if( !final ) goto tryagn; x.r = xsav.r; x.i = xsav.i; goto findon; } dx.r = (u.i * ud.i - u.r * ud.r)/mag; x.r += dx.r; dx.i = -(u.r * ud.i + u.i * ud.r)/mag; x.i += dx.i; if( (fabs(dx.i) + fabs(dx.r)) < 1.0e-6 ) goto lupdon; iter += 1; } /* while iter < 500 */ if( final ) goto lupdon; if( retry < 5 ) goto tryagn; return(3); lupdon: /* Swap original and reduced polynomials */ q = &xcof[nsav]; p = &cof[0]; for( j=0; j<=n2; j++ ) { cofj = *q; *q-- = *p; *p++ = cofj; } i = n; n = n1; n1 = i; if( !final ) { final = 1; if( fabs(x.i/x.r) < 1.0e-4 ) x.i = 0.0; xsav.r = x.r; xsav.i = x.i; goto finitr; /* do final iteration on original polynomial */ } findon: final = 0; if( fabs(x.i/x.r) >= 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<n; j++ ) { *(p+1) += cofj * *p - mag * *(p-1); p++; } setrut: root[nroot].r = x.r; root[nroot].i = x.i; nroot += 1; if( mag != 0.0 ) { x.i = -x.i; mag = 0; goto setrut; /* fill in the complex conjugate root */ } if( n > 0 ) goto nxtrut; return(0); }