NAG Library Function Document
nag_zeros_real_poly (c02agc)
1 Purpose
nag_zeros_real_poly (c02agc) finds all the roots of a real polynomial equation, using a variant of Laguerre's method.
2 Specification
#include <nag.h> 
#include <nagc02.h> 
void 
nag_zeros_real_poly (Integer n,
const double a[],
Nag_Boolean scale,
Complex z[],
NagError *fail) 

3 Description
nag_zeros_real_poly (c02agc) attempts to find all the roots of the
$n$th degree real polynomial equation
The roots are located using a modified form of Laguerre's method, originally proposed by
Smith (1967).
The method of Laguerre (see
Wilkinson (1965)) can be described by the iterative scheme
where
$H\left({z}_{k}\right)=\left(n1\right)\left[\left(n1\right){\left({P}^{\prime}\left({z}_{k}\right)\right)}^{2}nP\left({z}_{k}\right){P}^{\prime \prime}\left({z}_{k}\right)\right]$, and
${z}_{0}$ is specified.
The sign in the denominator is chosen so that the modulus of the Laguerre step at ${z}_{k}$, viz. $\leftL\left({z}_{k}\right)\right$, is as small as possible. The method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots.
The function generates a sequence of iterates
${z}_{1},{z}_{2},{z}_{3},\dots ,$ such that
$\leftP\left({z}_{k+1}\right)\right<\leftP\left({z}_{k}\right)\right$ and ensures that
${z}_{k+1}+L\left({z}_{k+1}\right)$ ‘roughly’ lies inside a circular region of radius
$\leftF\right$ about
${z}_{k}$ known to contain a zero of
$P\left(z\right)$; that is,
$\leftL\left({z}_{k+1}\right)\right\le \leftF\right$, where
$F$ denotes the Fejér bound (see
Marden (1966)) at the point
${z}_{k}$. Following
Smith (1967),
$F$ is taken to be
$\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(B,1.445nR\right)$, where
$B$ is an upper bound for the magnitude of the smallest zero given by
${r}_{1}$ is the zero
$X$ of smaller magnitude of the quadratic equation
and the Cauchy lower bound
$R$ for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation
Starting from the origin, successive iterates are generated according to the rule
${z}_{k+1}={z}_{k}+L\left({z}_{k}\right)$, for
$k=1,2,3,\dots $, and
$L\left({z}_{k}\right)$ is ‘adjusted’ so that
$\leftP\left({z}_{k+1}\right)\right<\leftP\left({z}_{k}\right)\right$ and
$\leftL\left({z}_{k+1}\right)\right\le \leftF\right$. The iterative procedure terminates if
$P\left({z}_{k+1}\right)$ is smaller in absolute value than the bound on the rounding error in
$P\left({z}_{k+1}\right)$ and the current iterate
${z}_{p}={z}_{k+1}$ is taken to be a zero of
$P\left(z\right)$ (as is its conjugate
${\stackrel{}{z}}_{p}$ if
${z}_{p}$ is complex). The deflated polynomial
$\stackrel{~}{P}\left(z\right)=P\left(z\right)/\left({zz}_{p}\right)$ of degree
$n1$ if
${z}_{p}$ is real (
$\stackrel{~}{P}\left(z\right)=P\left(z\right)/\left(\left({zz}_{p}\right)\left(z{\stackrel{}{z}}_{p}\right)\right)$ of degree
$n2$ if
${z}_{p}$ is complex) is then formed, and the above procedure is repeated on the deflated polynomial until
$n<3$, whereupon the remaining roots are obtained via the ‘standard’ closed formulae for a linear (
$n=1$) or quadratic (
$n=2$) equation.
4 References
Marden M (1966) Geometry of polynomials Mathematical Surveys 3 American Mathematical Society, Providence, RI
Smith B T (1967) ZERPOL: a zero finding algorithm for polynomials using Laguerre's method Technical Report Department of Computer Science, University of Toronto, Canada
Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford
5 Arguments
 1:
n – IntegerInput
On entry: $n$, the degree of the polynomial.
Constraint:
${\mathbf{n}}\ge 1$.
 2:
a[${\mathbf{n}}+1$] – const doubleInput

On entry: ${\mathbf{a}}\left[\mathit{i}\right]$ must contain ${a}_{\mathit{i}}$ (i.e., the coefficient of ${z}^{n\mathit{i}}$), for $\mathit{i}=0,1,\dots ,n$.
Constraint:
${\mathbf{a}}\left[0\right]\ne 0.0$.
 3:
scale – Nag_BooleanInput

On entry: indicates whether or not the polynomial is to be scaled. See
Section 9 for advice on when it may be preferable to set
${\mathbf{scale}}=\mathrm{Nag\_FALSE}$ and for a description of the scaling strategy.
Suggested value:
${\mathbf{scale}}=\mathrm{Nag\_TRUE}$.
 4:
z[n] – ComplexOutput

On exit: the real and imaginary parts of the roots are stored in
${\mathbf{z}}\left[\mathit{i}\right]\mathbf{.}\mathbf{re}$ and
${\mathbf{z}}\left[\mathit{i}\right]\mathbf{.}\mathbf{im}$ respectively, for
$\mathit{i}=0,1,\dots ,n1$. Complex conjugate pairs of roots are stored in consecutive pairs of
z; that is,
${\mathbf{z}}\left[i+1\right]\mathbf{.}\mathbf{re}={\mathbf{z}}\left[i\right]\mathbf{.}\mathbf{re}$ and
${\mathbf{z}}\left[i+1\right]\mathbf{.}\mathbf{im}={\mathbf{z}}\left[i\right]\mathbf{.}\mathbf{im}$
 5:
fail – NagError *Input/Output

The NAG error argument (see
Section 3.6 in the Essential Introduction).
6 Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
 NE_BAD_PARAM

On entry, argument $\u27e8\mathit{\text{value}}\u27e9$ had an illegal value.
 NE_INT_ARG_LT

On entry, ${\mathbf{n}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{n}}\ge 1$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call
and any array sizes. If the call is correct then please contact
NAG for
assistance.
 NE_POLY_NOT_CONV

The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact
NAG immediately, as some basic assumption for the arithmetic has been violated.
 NE_POLY_OVFLOW

The function cannot evaluate
$P\left(z\right)$ near some of its zeros without overflow. Please contact
NAG immediately.
 NE_POLY_UNFLOW

The function cannot evaluate
$P\left(z\right)$ near some of its zeros without underflow. Please contact
NAG immediately.
 NE_REAL_ARG_EQ

On entry, ${\mathbf{a}}\left[0\right]=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{a}}\left[0\right]\ne 0.0$.
7 Accuracy
All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed.
8 Parallelism and Performance
Not applicable.
If
${\mathbf{scale}}=\mathrm{Nag\_TRUE}$, then a scaling factor for the coefficients is chosen as a power of the base
$b$ of the machine so that the largest coefficient in magnitude approaches
$\mathit{thresh}={b}^{{e}_{\mathrm{max}}p}$. You should note that no scaling is performed if the largest coefficient in magnitude exceeds
$\mathit{thresh}$, even if
${\mathbf{scale}}=\mathrm{Nag\_TRUE}$. (
$b$,
${e}_{\mathrm{max}}$ and
$p$ are defined in
Chapter x02.)
However, with
${\mathbf{scale}}=\mathrm{Nag\_TRUE}$, overflow may be encountered when the input coefficients
${a}_{0},{a}_{1},{a}_{2},\dots ,{a}_{n}$ vary widely in magnitude, particularly on those machines for which
${b}^{4p}$ overflows. In such cases,
scale should be set to Nag_FALSE and the coefficients scaled so that the largest coefficient in magnitude does not exceed
${b}^{{e}_{\mathrm{max}}2p}$.
Even so, the scaling strategy used in nag_zeros_real_poly (c02agc) is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, you are recommended to scale the independent variable $\left(z\right)$ so that the disparity between the largest and smallest coefficient in magnitude is reduced. That is, use the function to locate the zeros of the polynomial $d\times P\left(cz\right)$ for some suitable values of $c$ and $d$. For example, if the original polynomial was $P\left(z\right)={2}^{100}+{2}^{100}{z}^{20}$, then choosing $c={2}^{10}$ and $d={2}^{100}$, for instance, would yield the scaled polynomial $1+{z}^{20}$, which is wellbehaved relative to overflow and underflow and has zeros which are ${2}^{10}$ times those of $P\left(z\right)$.
If the function fails with
NE_POLY_NOT_CONV,
NE_POLY_UNFLOW or
NE_POLY_OVFLOW, then the real and imaginary parts of any roots obtained before the failure occurred are stored in
z in the reverse order in which they were found. More precisely,
${\mathbf{z}}\left[{\mathbf{n}}1\right]\mathbf{.}\mathbf{re}$ and
${\mathbf{z}}\left[{\mathbf{n}}1\right]\mathbf{.}\mathbf{im}$ contain the real and imaginary parts of the 1st root found,
${\mathbf{z}}\left[{\mathbf{n}}2\right]\mathbf{.}\mathbf{re}$ and
${\mathbf{z}}\left[{\mathbf{n}}2\right]\mathbf{.}\mathbf{im}$ contain the real and imaginary parts of the 2nd root found, and so on. The real and imaginary parts of any roots not found will be set to a large negative number, specifically
$1.0/\left(\sqrt{2.0}\times {\mathbf{nag\_real\_safe\_small\_number}}\right)$.
10 Example
To find the roots of the 5th degree polynomial ${z}^{5}+2{z}^{4}+3{z}^{3}+4{z}^{2}+5z+6=0$.
10.1 Program Text
Program Text (c02agce.c)
10.2 Program Data
Program Data (c02agce.d)
10.3 Program Results
Program Results (c02agce.r)