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[],

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_FEWER_ZEROS_FOUND_1	Only # roots were found after "max_itn" = # iterations. The "root" vector will contain the value for machine infinity in the last # locations. A larger value for "max_itn" may help to find more roots.
IMSL_POLY_CONSTANT_ZERO	The polynomial is identically zero.
IMSL_POLY_CONSTANT_NONZERO	The polynomial is a nonzero constant.
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.