zeros_poly
Finds the zeros of a polynomial with real coefficients using the Jenkins-Traub, three-stage algorithm. Alternatively, the classical companion matrix method, or Aberth’s method according to an implementation by D. A. Bini, can be selected.
Synopsis
#include <imsl.h>
f_complex *imsl_f_zeros_poly (int ndeg, float coef[], …, 0)
The type d_complex function is imsl_d_zeros_poly.
Required Arguments
int ndeg (Input)
Degree of the polynomial. For the Jenkins-Traub method, ndeg <= 100 is required.
float coef[] (Input)
Array with ndeg + 1 components containing the coefficients of the polynomial in increasing order by degree. The polynomial is coef[n] zn + coef [n − 1] zn-1 + … + coef[0], where n = ndeg.
Return Value
A pointer to the complex array of zeros of the polynomial. To release this space, use imsl_free. If no zeros are computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
f_complex *imsl_f_zeros_poly (int ndeg, float coef[],
IMSL_COMPANION_METHOD,
IMSL_ABERTH_METHOD,
IMSL_MAX_ITERATIONS, int max_itn
IMSL_RETURN_USER, f_complex root[],
0)
Optional Arguments
IMSL_COMPANION_METHOD
Chooses the classical companion matrix method to find the polynomial roots.
By default, the Jenkins-Traub algorithm is selected.
IMSL_ABERTH_METHOD
Chooses Aberth’s method according to an implementation by Bini to find the polynomial roots. Note that this method is internally always used in double precision, even if zeros_poly is called in single precision.
By default, the Jenkins-Traub algorithm is selected.
IMSL_MAX_ITERATIONS, int max_itn (Input)
The maximum number of iterations allowed in Aberth’s method IMSL_ABERTH_METHOD.
Default: max_itn = 30.
IMSL_RETURN_USER, f_complex root[] (Output)
Array with ndeg components containing the zeros of the polynomial.
Description
The function imsl_f_zeros_poly computes the n zeros of the polynomial
where the coefficients ai for i = 0, 1, …, n are real and n is the degree of the polynomial.
The function imsl_f_zeros_poly uses the Jenkins-Traub, three-stage algorithm (Jenkins and Traub 1970; Jenkins 1975). The zeros are computed one at a time for real zeros or two at a time for a complex conjugate pair. As the zeros are found, the real zero, or quadratic factor, is removed by polynomial deflation.
Alternatively, the classical companion matrix method can be used through optional argument IMSL_COMPANION_METHOD. The companion matrix method is based on the fact that if Cp denotes the companion matrix associated with p(z), then det(zI - Cp) = p(z), where I is the n x n identity matrix. Thus, det(z0I - Cp) = 0 if and only if z0 is a zero of p(z). This implies that computing the eigenvalues of Cp will yield the roots of p(z). The companion matrix method is thought to be more robust than the Jenkins-Traub algorithm in most cases.
A second alternative is a variant of Aberth’s method based on an implementation by D. A. Bini (1996). This method can be used via optional argument IMSL_ABERTH_METHOD. It’s a numerically very stable Newton-type method with rapid convergence properties that can also be applied to high-degree polynomials.
For larger degree polynomials, use of the double-precision versions of the root-finding methods is recommended. This applies especially to the Jenkins-Traub method.
***********************************************************************
* All the software contained in this library is protected by copyright. *
* Permission to use, copy, modify, and distribute this software for any *
* purpose without fee is hereby granted, provided that this entire notice *
* is included in all copies of any software which is or includes a copy *
* or modification of this software and in all copies of the supporting *
* documentation for such software. *
***********************************************************************
* THIS SOFTWARE IS BEING PROVIDED “AS IS”, WITHOUT ANY EXPRESS OR IMPLIED *
* WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY *
* MEMBER OF THE EDITIORIAL BOARD OF THE JOURNAL “NUMERICAL ALGORITHMS” *
* NOR ITS EDITOR - IN - CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, *
* ANY MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK *
* OF USING THE SOFTWARE LIES WITH THE PARTY DOING SO. *
***********************************************************************
* ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE *
* ABOVE STATEMENT. *
***********************************************************************
Examples
Example 1
This example finds the zeros of the third-degree polynomial
p(z) = z3 − 3z2 + 4z – 2
where z is a complex variable.
#include <imsl.h>
#define NDEG 3
int main()
{
f_complex *zeros;
static float coeff[NDEG + 1] = {-2.0, 4.0, -3.0, 1.0};
zeros = imsl_f_zeros_poly(NDEG, coeff, 0);
imsl_c_write_matrix ("The complex zeros found are", 1, 3,
zeros, 0);
}
Output
The complex zeros found are
1 2 3
( 1, 0) ( 1, 1) ( 1, -1)
Example 2
The same problem is solved with the return option.
#include <imsl.h>
#define NDEG 3
int main()
{
f_complex zeros[3];
static float coeff[NDEG + 1] = {-2.0, 4.0, -3.0, 1.0};
imsl_f_zeros_poly(NDEG, coeff,
IMSL_RETURN_USER, zeros, 0);
imsl_c_write_matrix ("The complex zeros found are", 1, 3,
zeros, 0);
}
Output
The complex zeros found are
1 2 3
( 1, 0) ( 1, 1) ( 1, -1)
Warning Errors
IMSL_ZERO_COEFF | The # leading coefficients of the polynomial are equal to zero. The last # roots will be set to ("infinity",0.0) where "infinity" is the positive machine infinity. |
IMSL_ZERO_COEFF_1 | The leading coefficient of the polynomial is equal to zero. The last root will be set to ("infinity", 0.0) where "infinity" is the positive machine infinity. |
IMSL_FEWER_ZEROS_FOUND | Only # roots were found. The "root" vector will contain the value for machine infinity in the last # locations. |
IMSL_POLY_CONSTANT_ZERO | The polynomial is identically zero. |
IMSL_POLY_CONSTANT_NONZERO | The polynomial is a nonzero constant. |
IMSL_NO_CONVERGE_MAX_ITER | Failure to converge within "max_itn" = # iterations for at least one of the "nroot" = # roots. |
IMSL_POLY_COEFFS_TOO_LARGE | At least one of the coefficients of the polynomial is too large, overflow is likely. |
IMSL_POLY_ROOT_BOUNDS_FAIL | The computation of some inclusion radii for the roots of the polynomial may fail. |
IMSL_INV_POLY_ZEROS_TOO_SMALL | # zero(s) are too small to represent their inverses as complex numbers. They are replaced by small numbers. |
IMSL_POLY_ZEROS_TOO_SMALL | # zero(s) are too small to be represented by complex numbers. They are set to 0. |
IMSL_POLY_ZEROS_TOO_LARGE | # zero(s) are too big to be represented by complex numbers. They are set to a large number instead. |
