diff options
Diffstat (limited to 'libm/double/polrt.c')
-rw-r--r-- | libm/double/polrt.c | 227 |
1 files changed, 0 insertions, 227 deletions
diff --git a/libm/double/polrt.c b/libm/double/polrt.c deleted file mode 100644 index b1cd88087..000000000 --- a/libm/double/polrt.c +++ /dev/null @@ -1,227 +0,0 @@ -/* 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); -} |