The following code is to find a root for the polynomial, but it seems
to have the limitation on the order sicne the code as the follows
enountered a run time crash. Would you kind to let me know the bug of
this program and how to solve it.
/*const FLOAT *a, Coefficients
FLOAT *u, Real component of each root
FLOAT *v, Imaginary component of each root
FLOAT *conv, Convergence constant associated with each root
register long n, Degree of polynomial (order-1)
long maxiter, Maximum number of iterations
long fig The number of decimal figures to be computed
*/
//#include <fp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#ifdef DOUBLE_PRECISION
//# define FLOAT double
# define PARAMFLOAT double_t
#else /* DOUBLE_PRECISION */
# define FLOAT float
# define PARAMFLOAT float_t
#endif /* DOUBLE_PRECISION */
void FindPolynomialRoots( const float *a,float *u,float
*v,float*conv,register long n,long maxiter,long fig);
#define MAXN 16
int main (){
float *Coeff;
float *Real;
float *Imaginary;
float *Convergence;
register long Degree;
long Max_iter;
long decifig;
int i;
int seg_number=25;//max 18
float s0Val=0.001;
float wireLeni=0.08636;
Coeff=malloc(sizeof(float)*seg_number);
Real=malloc(sizeof(float)*(seg_number));
Imaginary=malloc(sizeof(float)*(seg_number));
Convergence=malloc(sizeof(double)*(seg_number));
Degree=seg_number-1;
Max_iter=1000;
decifig=12;
for(i=0;i<(seg_number-1);i++){
Real[i]=Imaginary[i]=Convergence[i]=0;
}
for (i=0; i<seg_number; i++){
Coeff[i]=s0Val;
printf("Coeff[%d] is %f\n",i,Coeff[i]);
}
Coeff[0]=Coeff[0]-wireLeni;
printf("Coeff[0] is %f\n",Coeff[0]);
FindPolynomialRoots(Coeff,Real,Imaginary,Convergence,Degree,Max_iter,decifig);
for(i=0;i<(seg_number-1);i++){
printf("Real %d is %f",i,Real[i]);
printf("Imag %d is %f",i,Imaginary[i]);
printf("\n");
}
return 0;
}
void FindPolynomialRoots( const float *a,float *u,float
*v,float *conv,register long n,long maxiter,long fig){
int i;
register int j;
FLOAT h[MAXN + 3], b[MAXN + 3], c[MAXN + 3], d[MAXN + 3], e[MAXN +
3];
/* [-2 : n] */
FLOAT K, ps, qs, pt, qt, s, rev, r;
int t;
FLOAT p, q;
/* Zero elements with negative indices */
b[2 + -1] = b[2 + -2] =
c[2 + -1] = c[2 + -2] =
d[2 + -1] = d[2 + -2] =
e[2 + -1] = e[2 + -2] =
h[2 + -1] = h[2 + -2] = 0.0;
/* Copy polynomial coefficients to working storage */
for (j = n; j >= 0; j--)
h[2 + j] = *a++; /* Note reversal of coefficients */
t = 1;
K = pow(10.0, (double)(fig)); /* Relative accuracy */
for (; h[2 + n] == 0.0; n--) { /* Look for zero high-order coeff.
*/
*u++ = 0.0;
*v++ = 0.0;
*conv++ = K;
}
INIT:
if (n == 0)
return;
ps = qs = pt = qt = s = 0.0;
rev = 1.0;
K = pow(10.0, (double)(fig));
if (n == 1) {
r = -h[2 + 1] / h[2 + 0];
goto LINEAR;
}
for (j = n; j >= 0; j--) /* Find geometric mean of coeff's */
if (h[2 + j] != 0.0)
s += log(fabs(h[2 + j]));
s = exp(s / (n + 1));
for (j = n; j >= 0; j--) /* Normalize coeff's by mean */
h[2 + j] /= s;
if (fabs(h[2 + 1] / h[2 + 0]) < fabs(h[2 + n - 1] / h[2 + n])) {
REVERSE:
t = -t;
for (j = (n - 1) / 2; j >= 0; j--) {
s = h[2 + j];
h[2 + j] = h[2 + n - j];
h[2 + n - j] = s;
}
}
if (qs != 0.0) {
p = ps;
q = qs;
} else {
if (h[2 + n - 2] == 0.0) {
q = 1.0;
p = -2.0;
} else {
q = h[2 + n] / h[2 + n - 2];
p = (h[2 + n - 1] - q * h[2 + n - 3]) / h[2 + n - 2];
}
if (n == 2)
goto QADRTIC;
r = 0.0;
}
ITERATE:
for (i = maxiter; i > 0; i--) {
for (j = 0; j <= n; j++) { /* BAIRSTOW */
b[2 + j] = h[2 + j] - p * b[2 + j - 1] - q * b[2 + j - 2];
c[2 + j] = b[2 + j] - p * c[2 + j - 1] - q * c[2 + j - 2];
}
if ((h[2 + n - 1] != 0.0) && (b[2 + n - 1] != 0.0)) {
if (fabs(h[2 + n - 1] / b[2 + n - 1]) >= K) {
b[2 + n] = h[2 + n] - q * b[2 + n - 2];
}
if (b[2 + n] == 0.0)
goto QADRTIC;
if (K < fabs(h[2 + n] / b[2 + n]))
goto QADRTIC;
}
for (j = 0; j <= n; j++) { /* NEWTON */
d[2 + j] = h[2 + j] + r * d[2 + j - 1];/* Calculate polynomial at r
*/
e[2 + j] = d[2 + j] + r * e[2 + j - 1];/* Calculate derivative at r
*/
}
if (d[2 + n] == 0.0)
goto LINEAR;
if (K < fabs(h[2 + n] / d[2 + n]))
goto LINEAR;
c[2 + n - 1] = -p * c[2 + n - 2] - q * c[2 + n - 3];
s = c[2 + n - 2] * c[2 + n - 2] - c[2 + n - 1] * c[2 + n - 3];
if (s == 0.0) {
p -= 2.0;
q *= (q + 1.0);
} else {
p += (b[2 + n - 1] * c[2 + n - 2] - b[2 + n] * c[2 + n - 3]) / s;
q += (-b[2 + n - 1] * c[2 + n - 1] + b[2 + n] * c[2 + n - 2]) / s;
}
if (e[2 + n - 1] == 0.0)
r -= 1.0; /* Minimum step */
else
r -= d[2 + n] / e[2 + n - 1]; /* Newton's iteration */
}
ps = pt;
qs = qt;
pt = p;
qt = q;
if (rev < 0.0)
K /= 10.0;
rev = -rev;
goto REVERSE;
LINEAR:
if (t < 0)
r = 1.0 / r;
n--;
*u++ = r;
*v++ = 0.0;
*conv++ = K;
for (j = n; j >= 0; j--) { /* Polynomial deflation by lin-nomial
*/
if ((d[2 + j] != 0.0) && (fabs(h[2 + j] / d[2 + j]) < K))
h[2 + j] = d[2 + j];
else
h[2 + j] = 0.0;
}
if (n == 0)
return;
goto ITERATE;
QADRTIC:
if (t < 0) {
p /= q;
q = 1.0 / q;
}
n -= 2;
if (0.0 < (q - (p * p / 4.0))) { /* Two complex roots */
*(u + 1) = *u = -p / 2.0;
u += 2;
s = sqrt(q - (p * p / 4.0));
*v++ = s;
*v++ = -s;
} else { /* Two real roots */
s = sqrt(((p * p / 4.0)) - q);
if (p < 0.0)
*u++ = -p / 2.0 + s;
else
*u++ = -p / 2.0 - s;
*u++ = q / u[-1];
*v++ = 0.0;
*v++ = 0.0;
}
*conv++ = K;
*conv++ = K;
for (j = n; j >= 0; j--) { /* Polynomial deflation by quadratic
*/
if ((b[2 + j] != 0.0) && (fabs(h[2 + j] / b[2 + j]) < K))
h[2 + j] = b[2 + j];
else
h[2 + j] = 0.0;
}
goto INIT;
}
#undef MAXN |
Clarification of Question by
july2003-ga
on
22 Jun 2003 21:46 PDT
Hi mathtalk-ga
Thanks for helping.
Would you able to fix it and let me know where did you fix (Like
highlighting it) since I am not familiar with this alogrithm and this
program is not originally developed by me.
Below is some infomation the author has provided in his code.
Thanks
/*******************************************************************************
* FindPolynomialRoots
*
* The Bairstow and Newton correction formulae are used for a
simultaneous
* linear and quadratic iterated synthetic division. The coefficients
of
* a polynomial of degree n are given as a[i] (i=0,i,..., n) where
a[0] is
* the constant term. The coefficients are scaled by dividing them by
* their geometric mean. The Bairstow or Newton iteration method will
* nearly always converge to the number of figures carried, fig,
either to
* root values or to their reciprocals. If the simultaneous Newton
and
* Bairstow iteration fails to converge on root values or their
* reciprocals in maxiter iterations, the convergence requirement will
be
* successively reduced by one decimal figure. This program
anticipates
* and protects against loss of significance in the quadratic
synthetic
* division. (Refer to "On Programming the Numerical Solution of
* Polynomial Equations," by K. W. Ellenberger, Commun. ACM 3 (Dec.
1960),
* 644-647.) The real and imaginary part of each root is stated as
u[i]
* and v[i], respectively, together with the corresponding constant,
* conv[i], used in the convergence test. This program has been used
* successfully for over a year on the Bendix G15-D (Intercard System)
and
* has recently been coded for the IBM 709 (Fortran system).
*/
|
Request for Question Clarification by
mathtalk-ga
on
24 Jun 2003 11:09 PDT
Yes. In C "float" means the single precision floating point numeric
datatype, and "double" means the double precision floating point
numeric datatype.
The usual "Pentium" Intel microprocessors and most others found in
PC's implement those datatypes according to an IEEE floating point
standard:
[IEEE Standard 754 Floating-Point]
http://research.microsoft.com/~hollasch/cgindex/coding/ieeefloat.html
The single precision "float" takes up 32 bits (4 bytes), some of which
is used to hold the floating point exponent and sign, with 23 bits
left to store the mantissa (fractional part of a number after
normalizing for powers of two). The double precision "double" takes
up 64 bits (8 bytes), and despite allowing more space for the
exponent, there is still 52 bits for a mantissa or fractional part.
The (relative) precision of these numeric datatypes depends on the
length of the mantissa, so "double" actually is more than twice as
precise as a "float" in that sense (52 bits being more than twice 23
bits).
In my opinion the "float" datatype is useless for your computational
purposes. The savings of "space" for the handful of variables can
hardly have any real merit, nor is there economy of time in performing
computations with the less precise "float" over "double". The
floating point arithmetic is almost certain to be done with a "math
coprocessor" instruction that takes away essentially all of the
performance penalty of working with doubles. [Graphics processing is
done with a lot of single precision computations performed in
parallel.]
Note that while you state the purpose of the code as "to find a root
for the polynomial," the comments in the code itself indicate that it
tries to find "roots" (plural). In particular the section of code
which will attempt to perform "deflation" will not produce accurate
roots after an "inaccurate" one is removed from the polynomial by
(synthetic) division.
If you want code that approximates to limited accuracy to a single
root of a polynomial, then I would not use the approach outlined in
this code. A simple bisection method would be more than adequate for
approximating a real root of a real polynomial to (say) five digits.
Complex roots of real polynomials are a bit trickier, which is where
Bairstow's method comes into play.
If your goal is rather to make the FindPolynomialRoots subroutine
work, so that it locates all polynomial roots as reliably as possible,
then perhaps some work on this program makes sense. After working
with the code to this point I feel the $10 list price you have
initially offered is not sufficient for the time it will take to do
the debugging needed to get the code working in reliable fashion.
Although another researcher may take an interest in working it at the
offered price, you may want to review the Google Answer pricing
guidelines here:
http://answers.google.com/answers/pricing.html
Multi-part questions are normally priced at $50 and above, and I'm
inclined to think of debugging a program of this complexity as a
multi-part question.
regards, mathtalk-ga
|
