summaryrefslogtreecommitdiff
path: root/libm/double/polrt.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/polrt.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/double/polrt.c')
-rw-r--r--libm/double/polrt.c227
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);
-}