# NAG FL Interfacec02aaf (poly_​complex_​fpml)

## ▸▿ Contents

Settings help

FL Name Style:

FL Specification Language:

## 1Purpose

c02aaf finds all the roots of a complex polynomial equation, using a fourth-order convergent modification of Laguerre's method.

## 2Specification

Fortran Interface
 Subroutine c02aaf ( a, n, z, berr, cond, conv,
 Integer, Intent (In) :: n, itmax, polish Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: conv(n) Real (Kind=nag_wp), Intent (Out) :: berr(n), cond(n) Complex (Kind=nag_wp), Intent (In) :: a(0:n) Complex (Kind=nag_wp), Intent (Out) :: z(n)
#include <nag.h>
 void c02aaf_ (const Complex a[], const Integer *n, const Integer *itmax, const Integer *polish, Complex z[], double berr[], double cond[], Integer conv[], Integer *ifail)
The routine may be called by the names c02aaf or nagf_zeros_poly_complex_fpml.

## 3Description

c02aaf attempts to find all the roots of the $n$th degree complex polynomial equation
 $P(z) = a0zn + a1zn-1 + a2zn-2 + ⋯ + an-1 z + an = 0 .$
The roots are located using a modified form of Laguerre's method, as implemented by Cameron (2018). An implicit deflation strategy is employed, which allows for high accuracy even when solving high degree polynomials. Linear ($n=1$) and quadratic ($n=2$) problems are obtained via the 'standard' closed formulae.
First, initial estimates of the roots are made using a method proposed by Bini (1996), which selects complex numbers along circles of suitable radii. Updates to each root approximation ${z}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,n$, are then made using the iterative formula from Petković et al. (1997):
 $z^j = zj - n Gj ± (n-1) (nHj-Gj2) ,$
where
 $Gj = P′(zj) P(zj) - ∑ i=1; ​i≠j n 1 (zj-zi) and Hj = - ( P′(zj) P(zj) ) ′ - ∑ i=1; ​i≠j n 1 (zj-zi) 2 .$
The nearest ${\stackrel{^}{z}}_{j}$ to the current approximation is used, by selecting the sign that maximizes the denominator of the correction term. The subtraction of the sum terms when computing ${G}_{j}$ and ${H}_{j}$ determines the implicit deflation strategy of the modified Laguerre method.
The relative backward error of a root approximation ${z}_{j}$ is given by
 $η(zj) = |P(zj)| α(zj) ,$
where $\alpha$ is the perturbed polynomial
 $α(z) = ∑ k=0 n (3.8⁢k+1) |an-k| |z|k .$
A root approximation is deemed to have converged if $\eta \left({z}_{j}\right)\le 2\text{​ ​}×$ machine precision, at which point updates of that root cease. If the stopping criterion holds, then the computed root is the exact root of a polynomial whose coefficients are no more perturbed than the floating-point computation of $P\left({z}_{j}\right)$.
The condition number of each root is also computed, as a measure of sensitivity to changes in the coefficients of the polynomial. It is given by
 $κ(zj) = α(zj) |zj| |P′(zj)| .$
Root approximations can be further refined with optional 'polishing' processes. A simple polishing process is provided that carries out a single iteration of Newton's method, which proves quick yet often effective. Alternatively, a compensated polishing process from Cameron and O'Neill (2019) can be applied. This iterative method combines the implicit deflation of the modified Laguerre method, with the accuracy of evaluating polynomials and their derivatives using the compensated Horner's method from Graillat et al. (2005). Compensated polishing yields approximations with a limiting accuracy as if computed in twice the working precision.
It is recommended that you read Section 9.1 for advice on selecting an appropriate polishing process.
Bini D A (1996) Numerical computation of polynomial zeros by means of Aberth's method Numerical Algorithms 13 179–200 Springer US
Cameron T R (2018) An effective implementation of a modified Laguerre method for the roots of a polynomial Numerical Algorithms Springer US https://doi.org/10.1007/s11075-018-0641-9
Cameron T R and O'Neill A (2019) On a compensated polishing technique for polynomial root solvers To be published
Graillat S, Louvet N, and Langlois P (2005) Compensated Horner scheme Technical Report Université de Perpignan Via Domitia
Petković M, Ilić S, and Tričković S (1997) A family of simultaneous zero-finding methods Computers & Mathematics with Applications (Volume 34) 10 49–59 https://doi.org/10.1016/S0898-1221(97)00206-X
Wilkinson J H (1959) The evaluation of the zeros of ill-conditioned polynomials. Part I Numerische Mathematik (Volume 1) 1 150–166 Springer-Verlag

## 5Arguments

1: $\mathbf{a}\left(0:{\mathbf{n}}\right)$Complex (Kind=nag_wp) array Input
On entry: ${\mathbf{a}}\left(\mathit{i}\right)$ must contain the coefficient of ${z}^{n-\mathit{i}}$, for $\mathit{i}=0,1,\dots ,n$.
Constraint: ${\mathbf{a}}\left(0\right)\ne \left(0.0,0.0\right)$.
2: $\mathbf{n}$Integer Input
On entry: $n$, the degree of the polynomial.
Constraint: ${\mathbf{n}}\ge 1$.
3: $\mathbf{itmax}$Integer Input
On entry: the maximum number of iterations to be performed.
Suggested value: ${\mathbf{itmax}}=30$.
Constraint: ${\mathbf{itmax}}\ge 1$.
4: $\mathbf{polish}$Integer Input
On entry: specifies the polishing technique used to refine root approximations.
${\mathbf{polish}}=0$
No polishing.
${\mathbf{polish}}=1$
Single iteration of Newton's method.
${\mathbf{polish}}=2$
Iterative refinement using the compensated Horner's method.
Suggested value: ${\mathbf{polish}}=1$.
Constraint: ${\mathbf{polish}}=0$, $1$ or $2$.
5: $\mathbf{z}\left({\mathbf{n}}\right)$Complex (Kind=nag_wp) array Output
On exit: ${\mathbf{z}}\left(\mathit{j}\right)$ holds the approximation of the root ${z}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,n$.
6: $\mathbf{berr}\left({\mathbf{n}}\right)$Real (Kind=nag_wp) array Output
On exit: ${\mathbf{berr}}\left(\mathit{j}\right)$ holds the relative backward error, $\eta \left({z}_{\mathit{j}}\right)$, for $\mathit{j}=1,2,\dots ,n$.
7: $\mathbf{cond}\left({\mathbf{n}}\right)$Real (Kind=nag_wp) array Output
On exit: ${\mathbf{cond}}\left(\mathit{j}\right)$ holds the condition number, $\kappa \left({z}_{\mathit{j}}\right)$, for $\mathit{j}=1,2,\dots ,n$.
8: $\mathbf{conv}\left({\mathbf{n}}\right)$Integer array Output
On exit: ${\mathbf{conv}}\left(\mathit{j}\right)$ indicates the convergence status of the root approximation, ${z}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,n$.
${\mathbf{conv}}\left(j\right)\ge 0$
Successfully converged after ${\mathbf{conv}}\left(j\right)$ iterations.
${\mathbf{conv}}\left(j\right)=-1$
Failed to converge after itmax iterations.
${\mathbf{conv}}\left(j\right)=-2$
Overflow was encountered.
Note: if ${\mathbf{polish}}=2$, conv refers to convergence in the compensated polishing process.
9: $\mathbf{ifail}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $-1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value $-1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. When the value $-\mathbf{1}$ or $\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
On entry, the complex variable ${\mathbf{a}}\left(0\right)=\left(0.0,0.0\right)$.
${\mathbf{ifail}}=2$
On entry, ${\mathbf{n}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{n}}\ge 1$.
${\mathbf{ifail}}=3$
On entry, ${\mathbf{itmax}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{itmax}}\ge 1$.
${\mathbf{ifail}}=4$
On entry, ${\mathbf{polish}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{polish}}=0$, $1$ or $2$.
${\mathbf{ifail}}=5$
Convergence has failed for at least one root approximation.
Check conv and consider increasing itmax.
${\mathbf{ifail}}=6$
c02aaf encountered overflow during at least one root approximation.
Check conv and consider scaling the polynomial (see Section 9.2).
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed.

## 8Parallelism and Performance

c02aaf is not threaded in any implementation.

### 9.1Selecting a Polishing Process

The choice of polishing technique ultimately depends on two factors: how well conditioned the problem is, and a preference between run time and accuracy. For a detailed analysis of the polishing techniques, see Cameron and O'Neill (2019).
Well-conditioned Problems
Simple polishing is effective in reducing the error in approximations of well-conditioned roots, doing so with a negligible increase in run time. Compensated polishing has comparable accuracy, but it is approximately ten times slower than when using simple polishing.
Simple polishing (${\mathbf{polish}}=1$) is recommended for well-conditioned problems.
Ill-conditioned Problems
There is a dramatic difference in accuracy between the two polishing techniques for ill-conditioned polynomials. Unpolished approximations are inaccurate and simple polishing often proves ineffective. However, compensated polishing is able to reduce errors by several orders of magnitude.
Compensated polishing (${\mathbf{polish}}=2$) is highly recommended for ill-conditioned problems.

### 9.2Scaling the Polynomial

c02aaf attempts to avoid overflow conditions where possible. However, if the routine fails with ${\mathbf{ifail}}={\mathbf{6}}$, such conditions could not be avoided for the given polynomial. Use conv to identify the roots for which overflow occurred, as other approximations may still have succeeded.
Extremely large and/or small coefficients are likely to be the cause of overflow failures. 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 routine to locate the zeros of the polynomial $sP\left(cz\right)$ for some suitable values of $c$ and $s$. For example, if the original polynomial was $P\left(z\right)={2}^{-100}i+{2}^{100}{z}^{20}$, then choosing $c={2}^{-10}$ and $s={2}^{100}$, for instance, would yield the scaled polynomial $i+{z}^{20}$, which is well-behaved relative to overflow and has zeros which are ${2}^{10}$ times those of $P\left(z\right)$.

## 10Example

The example program for c02aaf demonstrates two problems, given in the subroutines ex1_basic and ex2_polishing. (Note that by default, the second example is switched off because the results may be machine dependent. Edit the program in the obvious way to switch it on.)
Example 1: Basic Problem
This example finds the roots of the polynomial
 $a0 z5 + a1 z4 + a2 z3 + a3 z2 + a4 z+a5 = 0 ,$
where ${a}_{0}=\left(5.0+6.0i\right)$, ${a}_{1}=\left(30.0+20.0i\right)$, ${a}_{2}=-\left(0.2+6.0i\right)$, ${a}_{3}=\left(50.0+100000.0i\right)$, ${a}_{4}=-\left(2.0-40.0i\right)$ and ${a}_{5}=\left(10.0+1.0i\right)$.
Example 2: Polishing Processes
This example finds the roots of a polynomial of the form
 $(z-1) (z-2) ⋯ (z-n) = 0 ,$
first proposed by Wilkinson (1959) as an example of polynomials with ill-conditioned roots, that are sensitive to small changes in the coefficients. A polishing mode is demonstrated with $n=10$, and the maximum forward and relative backward errors of the approximations displayed.

### 10.1Program Text

Program Text (c02aafe.f90)

### 10.2Program Data

Program Data (c02aafe.d)

### 10.3Program Results

Program Results (c02aafe.r)