Noname manuscript No.
(will be inserted by the editor)

An Effective Implementation of a Modified Laguerre
Method for the Roots of a Polynomial
Thomas R. Cameron

Received: date / Accepted: date

Abstract Two common strategies for computing all roots of a polynomialwith Laguerre’s method are explicit deflation and Maehly’s procedure. The
former is only a semi-stable process and is not suitable for solving large degree
polynomial equations. In contrast, the latter implicitly deflates the polynomialusing previously accepted roots and is, therefore, a more practical strategy
for solving large degree polynomial equations. However, since the roots of apolynomial are computed sequentially, this method cannot take advantage ofparallel systems. In this article, we present an implementation of a modified
Laguerre method for the simultaneous approximation of all roots of a polynomial. We provide a derivation of this method along with a detailed analysis
of our algorithm’s initial estimates, stopping criterion, and stability. Finally,
the results of several numerical experiments are provided to verify our analysis
and the effectiveness of our algorithm.
Keywords Laguerre’s method · polynomial roots · mathematical software
Mathematics Subject Classification (2010) 26C10 · 65H04 · 65Y20
1 Introduction
In the mid-1800’s, Laguerre first published his root-finding method in Racinesd’une Équation Algébrique and it was later republished in Oeuvres de Laguerre
in 1898 [21]. Since then, research on modifications, generalizations, and analysis of Laguerre’s method has permeated the mathematical literature [10, 12,
17, 20, 24, 34, 36, 38]. In addition, Laguerre’s method has been applied to thematrix and polynomial eigenvalue problem [15, 22, 23, 30].
Thomas R. CameronDavidson College
Mathematics and Computer Science Department
Tel.: +1 (704)-894-2241
E-mail: thcameron@davidson.edu

2Thomas R. Cameron

Despite a significant amount of attention in the literature, we are only
aware of two popular software packages that employ Laguerre’s method tosolve for the roots of a polynomial: zroots from Numerical Recipes and C02AFF
from the NAG library [13, 26, 28, 35]. Both zroots and C02AFF use forward
deflation to ensure that subsequent iterations tend to a different root. However, zroots does nothing to ensure that the roots are computed in ascendingorder of absolute value, which is vital for the backward stability of forward
deflation [40]. In contrast, the NAG library routine C02AFF uses a modified
Laguerre method that was first proposed by Smith and guarantees the roots
are computed in ascending order of absolute value [38]. For this reason, theroutine C02AFF is more reliable than zroots.
Modern-day polynomial root problems, such as those that stem from computer algebra applications or digital signal processing, require software thatcan solve for the roots of a polynomial of degree well above 100 and sometimes on the order of several thousand or higher [29, 37]. Therefore, an explicit
deflation strategy such as forward deflation used by zroots and CO2AFF isnot suitable for solving today’s polynomial root problems. Indeed, explicit deflation is at best a semi-stable process since the cumulative effect after many
steps may result in a significant loss of accuracy [19]. We tested this notion on
the C02AFF routine using polynomials with random coefficients; for polynomials of degree higher than 500, the routine would often fail to converge andcontinue to run for minutes on end. This test serves as the impetus for our
research on implicit deflation strategies that help make Laguerre’s method asuitable option for solving large degree polynomial equations.
There have been several implicit deflation strategies proposed for Laguerre’s
method [17, 25, 33, 34]. They differ from the explicit deflation strategies since
they result in a modification of Laguerre’s method rather than a change to theunderlying polynomial. The modifications proposed in [17, 34] are mathematically interesting: They have local fourth-order, or higher, convergence and canbe run in parallel. However, these methods are not reliable in practice. In [24],it was noted that convergence could be quite slow when two approximations
are closer to each other than to a root. This situation frequently arises when
dealing with high degree polynomials. Therefore, these methods may only besuitable for iterative refinement. Modifying Laguerre’s method using Maehly’s
procedure [25] is far more practical and has been used to solve matrix andpolynomial eigenvalue problems [22, 30]. However, this method computes theroots of a polynomial sequentially and, therefore, cannot take advantage ofparallel systems.
In contrast, the family of methods derived in [33] allow for the simultaneous
approximation of all roots of a polynomial. These methods have local fourthorder convergence and computationally verifiable conditions that guarantee
their convergence for simple roots [32]. We are particularly interested in themember of this family that is analogous to Laguerre’s method. We will demonstrate that this method has a workload that is well-suited for data-parallelism
and use it to develop an algorithm that is effective for solving large degree
polynomial equations.

Implementation of a Modified Laguerre Method3

The outline of this article is as follows: In Section 2 we derive an instance ofthe modified Laguerre method from [33] that we use to form an algorithm forthe computation of all roots of a polynomial. In Sections 2.1–2.3, we provide a
detailed analysis of the algorithm’s initial approximations, stopping criterion,and stability. Finally, in Section 3, the results of several numerical experiments
are provided to verify our analysis and the effectiveness of our algorithm.

2 Derivation of the Method
Consider the polynomial p in variable z defined byp(z) = a0 + a1 z + · · · + am z m ,(1)

where a0 am 6= 0. If a0 = 0, then deflate the polynomial and apply the following method to the resulting polynomial. Denote by {z1 , . . . , zm } the current
approximations to the roots {ζ1 , . . . , ζm } of p. We are interested in updatingeach approximation zj until it is “close enough” to the root ζj .
Laguerre’s method performs this update by solving a quadratic equation,
which we denote by Qj (z) = 0. For a derivation of Qj , see [36]. Note that if
the roots of p are real, then both solutions of Qj (z) = 0 are guaranteed to
be closer to ζj than the current approximation zj . This is what accounts forthe global convergence of Laguerre’s method when all roots are real. While
global convergence is not guaranteed when the roots are complex, it is our
experience that Laguerre’s method often performs just as well in this case. A
similar sentiment was noted by Parlett in [30].
We denote the roots of Qj by ẑj and note thatẑj = zj −

m
Gj ±

q

(m − 1)(mHj − G2j )

,

(2)

where
0
m
1
p (zj ) X=
and Hj = −Gj =p(zj )
(z
−
ζi )
j
i=1

0

p (zj )p(zj )

!0
=

m
X

1
. (3)
(z
−
ζi )2
j
i=1

The update to zj is given by the nearest ẑj ; that is, we select the sign thatmaximizes the denominator of the fraction in (2). We reference this update as
the correction term of zj .
In [33], two one-parameter families of iteration functions are derived fromthe classical Weierstrass’ method and the Hansen-Patrick formula [16]. For aparticular choice of the parameter, a modified Laguerre method is obtained forthe simultaneous approximation of all roots of a polynomial. The derivationn [33] allows for repeated roots, but we assume that all roots are simple since
we cannot expect to know the multiplicity a priori.

4Thomas R. Cameron

In the case of simple roots, the modified Laguerre method from [33] isequivalent to making the following changes to the terms in (3):0m
p (zj ) X
1
Gj =−and Hj = −p(zj )
(z
−
zi )
j
i=1

0

p (zj )p(zj )

i6=j

!0
−

m
X

1
. (4)
(z
−
zi )2
j
i=1
i6=j

Moreover, this is equivalent to applying Laguerre’s method to the functionp(z),i=1 (z − zi )fj (z) = Qm

(5)

i6=j

when computing the correction term of zj . Note that this is a classical approach
for implementing implicit deflation and it is at the core of the Aberth-EhrlichBörsch-Supan method [1, 7, 11]. Essentially, we are creating poles at all otherroot approximations and therefore avoiding unnecessary multiple convergences
to the same root.
As noted in Section 1, the method from [33] has strong virtues including
local fourth-order convergence and a parallel workload. We use this method
to form an algorithm for the computation of all roots of a polynomial.
The outline of our algorithm is as follows: We begin with the initial approximations whose computation is discussed in Section 2.1. Then, for a fixed number of iterations, we loop through the array of approximations (z1 , . . . , zm ). If
zj is not close enough to ζj , which is quantified in Section 2.3, then we update
the approximation using (2) and (4), as described in Section 2.2. The loop
through the array of approximations is considered one iteration of our algorithm. Each iteration can be implemented sequentially or in parallel and bothstyles are outlined in Algorithm 1 and Algorithm 2, respectively. The former
has the benefit of a slightly improved convergence rate, whereas the latter will
have the obvious benefit of speedup.
Algorithm 1 Sequential Style(z1 , . . . , zm ) ← initial approximations via Algorithm 3
while i < itmax dofor j = 1 to m do
if zj is not close enough to ζj thenCompute roots of Qj (λ) via (2), using (4)zj ← root that maximizes denominator of (2)end ifend fori←i+1
end while

The parfor loops in Algorithm 2 are well-suited for data-parallelism. Thatis, the parallel tasks are executing the same function across different partsof the dataset. The actual implementation and the effectiveness thereof depends on the appropriate system support, such as a multiprocessor and anImplementation of a Modified Laguerre Method5

Algorithm 2 Parallel Style(z1 , . . . , zm ) ← initial approximations via Algorithm 3
while i < itmax doparfor j = 1 to m do
if zj is not close enough to ζj thenCompute Gj and Hj via (4) and store
end ifend parfor
parfor j = 1 to m do
if zj is not close enough to ζj thenUse Gj and Hj to compute roots of Qj (λ) via (2)zj ← root that maximizes denominator of (2)end ifend parfori←i+1
end while

operating system that supports multithreading. For now, our focus is on theimplementation of the modified Laguerre method from [33] for solving large
degree polynomial equations. We consider the parallel implementation a future step in our work, at which point we will consider other features such as
adaptivity and multi-precision as is done with MPSolve [6].

2.1 Initial Approximations
As with all iterative methods, the performance of the modified Laguerre methodfrom [33] depends considerably upon the quality of initial approximations. We
note that this method has a computationally verifiable condition for the guaranteed convergence to simple roots [32, Theorem 5.2]. However, this does little
for us in practice since if the test condition fails we are given no information
on how to improve our initial approximations. In this section, we outline a procedure originally due to Bini [4] for the computation of initial approximations
that, in practice, almost always converge to the roots of a polynomial.
In essence, this procedure selects complex numbers along circles of suitableradii; what constitutes suitable radii can be formalized by applying Pellet’sTheorem [31]. We state the theorem below and note that an elementary proof
can be found in [39].
Theorem 1 Let p be a polynomial as defined in (1). For each k such thatak 6= 0, consider the equation|ak |z k =

m
X
i=0
i6=k

|ai |z i .

(6)(i.) If k = 0 there exists one real positive solution s0 , and p has no roots ofmoduli less than s0 .

6Thomas R. Cameron(ii.) If 0 < k < m there are either no real positive solutions or two realpositive solutions tk ≤ sk . In the latter case, p has no roots in the open
annulus A(tk , sk ) and exactly k roots of moduli less than or equal to tk .
(iii.) If k = m there exists one real positive solution tm , and p has no rootsof moduli greater than tm .
Let 0 = k1 < k2 < · · · < kq = m be the values of k for which (6) has
positive solution(s) and let
s0 = sk1 ≤ tk2 ≤ sk2 ≤ · · · ≤ tkq−1 ≤ skq−1 ≤ tkq = tm
be these solutions. Then, by [4, Theorem 5], any polynomial in the class(m)X
i
P(p) =
bi z : |bi | = |ai |i=0

has (ki+1 − ki ) zeros in the closed annulus Ā(ski , tki+1 ), for i = 1, . . . , q − 1,and no zeros in the open annulus A(tki , ski ), for i = 1, . . . , q.
While these inclusion results can be used to determine initial approximations for our iterative method, forming positive solutions of (6) requires solving
several polynomial equations. Therefore, we seek a strategy with a smaller timecomplexity. To this end, defineaiak

1/(k−i)uk = maxaiak

1/(k−i)vk = mini<ki>kai 6=0,k = 1, . . . , m,,

k = 0, . . . , m − 1,

(7)



ai
u0 = 1/ 1 + max,
i>0 a0ai.
v0 = 1 + max
i<m amThen it follows from [4, Theorem 4] that the positive solutions of (6) satisfy
uki ≤ tki ≤ ski ≤ vki ,i = 1, . . . , q.

(8)Consider the set C = {(i, log |ai |), i = 0, 1, . . . , m} and denote byγ(C) = {(k̂i , log |ak̂i |), 0 = k̂1 < k̂2 < · · · < k̂q̂ = m}
the upper envelope of the convex hull of C. By [4, Theorem 6] we know that thevertices of γ(C) satisfy ul ≤ vl if and only if l ∈ {k̂1 , . . . , k̂q̂ }. Therefore, by (8),we have {k1 , . . . , kq } ⊆ {k̂1 , . . . , k̂q̂ }, where the containing set can be computed
efficiently. In fact, a divide and conquer method is used in [4] to compute theconvex hull γ(C) in O(m log m) time. Since the set C is ordered with respect tothe first coordinate, we use Andrew’s monotone chain algorithm to computeγ(C), which, in this case, has time complexity O(m) [2]. The pseudo-code forcomputing the initial approximations is given in Algorithm 3, where i is theimaginary unit and σ is any nonzero number [4].

Implementation of a Modified Laguerre Method

7

Algorithm 3 initial approximations
Compute k̂1 , . . . , k̂q̂ via Andrew’s monotone chain algorithm.
for i = 1 to q̂ − 1 don = k̂i+1 − k̂i
uk̂

i+1

=

ak̂

1/n

i

ak̂

i+1

for j = 1 to n do
2π2πie( n j+ m +σ)izk̂ +j = uk̂
ii+1
end forend for2.2 Computing the Correction Term
The update of the root approximation zj requires the computation of thecorrection term ẑj via (2), where Gj and Hj are defined in (4). In this section,
we outline a numerically stable method for computing the correction term.
Both summations in (4) represent a backward stable computation. Indeed,each summand can be computed in at most 3 flops, each of which is exactly
rounded assuming the IEEE 754 standard [14]. Moreover, the summation iswell-known to be backward stable [18, Section 4.2]. Finally, polynomial evaluation via the Ruffini-Horner rule is a backward stable computation. We outline
the latter result as it will be used when deriving our stopping criterion.
Let fl(p(ξ)) denote the computed value of p evaluated at ξ ∈ C using
floating-point arithmetic with machine precision µ. Then, if the Ruffini-Horner
rule is used, it follows from [4, Theorem 7] that we havefl(p(ξ)) =mXai (1 + i )ξ i ,

(9)

i=0

√
where |i | < ((2 2 + 1)i + 1)µ + O(µ2 ). A similar result holds for the first andsecond derivative of p.
The Ruffini-Horner method is prone to overflow, for instance, when a large
degree polynomial with positive coefficients is evaluated at ξ where |ξ| > 1. For
this reason, we introduce the reversal polynomial pR in the variable z definedbypR (z) = a0 z m + · · · + am−1 z + am .
Let z 6= 0 and ρ = 1/z, then it is easy to verify that the following holds true
!
0
0
p (z)
pR (ρ)=ρ m−ρ,p(z)pR (ρ)

(10)
!0
!0 
0
0
0
p (z)
pR (ρ) pR (ρ)
22
−= ρ m − 2ρ
−ρ
.
p(z)
pR (ρ)pR (ρ)
When |ξ| ≤ 1, we apply the Ruffini-Horner rule to the polynomial p, and its
first and second derivative. Otherwise, we apply the Ruffini-Horner rule to8Thomas R. Cameron

the reversal polynomial pR , and its first and second derivative, evaluated atρ = 1/ξ. In either case, it follows from (10) that we can easily obtain thedesired values needed for computing the correction term.
Now that the main components of Algorithm 1 have been established, weanalyze our algorithm’s time and space complexity. The initial approximations are computed via Algorithm 3 in O(m) time. Then, each approximation
is updated in O(m) time using the Ruffini-Horner rule as described above.
Assuming that itmax is fixed, i.e., does not depend on the degree of the polynomial, it follows that Algorithm 1 can approximate all roots of a polynomialin O(m2 ) time. Finally, we note that our algorithm has O(m) space complexity
since the required storage is the coefficients of the polynomial, its approximateroots, and the corresponding backward error and condition number.
2.3 Stopping Criterion
Let ξ be an approximation of the root ζ for the polynomial p as definedin (1), where all computations have been done in floating-point arithmetic withmachine precision µ. In this section, we derive the backward error of ξ andthe condition number of ζ, which are used to establish our stopping criterion
and provide a first-order bound on the forward error in our approximation.
Both results in Theorem 2 and Theorem 3 can be found in [8]. We provide
elementary proofs for completeness and clarity.
For i = 0, 1, . . . , m, denote by ei an arbitrary representation of the tolerances against which the perturbations ∆ai to ai will be measured. Then theperturbed polynomial in the variable z is defined by∆p(z) = ∆a0 + ∆a1 z + · · · + ∆am z m ,(11)and for convenience we define the sign of a complex number z by(z̄/|z|, z 6= 0,
sgn(z) =0,
z = 0.
The backward error of ξ represents the smallest perturbations ∆ai forwhich (p + ∆p)(ξ) = 0. Therefore, a natural definition of the backward errorof ξ is given byη(ξ) = min{ : (p + ∆p)(ξ) = 0, |∆ai | ≤ |ei |, i = 0, 1, . . . , m}.
Note that when ei = 1 we are measuring the absolute backward error andwhen ei = ai we are measuring the relative backward error.
Theorem 2 The backward error η(ξ) is given byη(ξ) =
where α(ξ) =Pmii=0 |ei ||ξ| .

|p(ξ)|
,
α(ξ)

(12)

Implementation of a Modified Laguerre Method

9

Proof First we show that the right hand side is a lower bound for η(ξ). To thisend, note that|p(ξ)| = |∆p(ξ)|
m
m
X
X
≤
|ξ|i |∆ai | ≤ |ξ|i |ei |.
i=0

i=0

Next, we show that this lower bound is attained by the perturbations∆ai = −1sgn(ξ i )|ei |p(ξ).
α(ξ)

To this end, note that∆p(ξ) = −mX
1
sgn(ξ i )|ei |p(ξ)ξ i
α(ξ)
i=0

=−

p(ξ) X|ei ||ξ|i = −p(ξ).
α(ξ) i=0m

Therefore, (p + ∆p)(ξ) = 0 and |∆ai |/|ei | = |p(ξ)|/α(ξ) for all i = 0, 1, . . . , m.
tuMotivated by
thebackward
error
analysisofthe
Ruffini-Horner
rule
in
(9),
√
we let ei = ((2 2 + 1)i + 1)ai . Then, if the approximate root satisfies
η(ξ) ≤ µ(13)

it follows that ξ is a root of (p + ∆p), where ∆p is defined in (11) and its
coefficients ∆ai are no bigger than i from (9), for i = 0, 1, . . . , m. Conversely,if η(ξ) > µ, then ξ being a root of (p+∆p) implies that ∆ai is bigger than i , forsome i = 0, 1, . . . , m. Therefore, (13) denotes a viable stopping criterion thatguarantees iterations do not terminate until fl(p(ξ)) is no longer a reliable
computation for updating the approximate root ξ. Furthermore, if ξ is anapproximate root that satisfies (13), then the computation of ξ was backward
stable. There is no guarantee of backward stability for all polynomials sinceconvergence is not guaranteed. However, in practice, our algorithm almost
always performs in a backward stable manner, see Section 3.
As noted in Section 2.2, if |ξ| > 1 then we apply the Ruffin-Horner ruleto the reversal polynomial pR at ρ = 1/ξ in order to compute the correction
term. In this case, the backward error can be computed as follows
η(ξ) =
where αR (ρ) =|pR (ρ)|,αR (ρ)Pm(14)i
i=0 |em−i ||ρ| . This result is easily verified by noting that

|ρ|m |p(1/ρ)||pR (ρ)|=αR (ρ)|ρ|m α(1/ρ)|p(ξ)|
=
.
α(ξ)

10

Thomas R. Cameron

The condition number of a root is a measurement of its sensitivity tochanges in the coefficients of p. Therefore, a natural definition of the condition
number of a nonzero simple root ζ is given by|∆ζ|: (p + ∆p)(ζ + ∆ζ) = 0, |∆ai | ≤ |ei |, i = 0, 1, . . . , m .
κ(ζ, p) = lim sup
|ζ|→0Theorem 3 The condition number κ(ζ, p) is given byκ(ζ, p) =where α(ζ) =α(ζ),|ζ||p0 (ζ)|

(15)

Pm

i
i=0 |ei ||ζ| .

Proof Note the following Taylor series expansions
0p(ζ + ∆ζ) = p(ζ) + ∆ζp (ζ) + O(∆ζ 2 )and0∆p(ζ + ∆ζ) = ∆p(ζ) + ∆ζ∆p (ζ) + O(∆ζ 2 ).
Then the constraint (p + ∆p)(ζ + ∆ζ) = 0 implies that∆p(ζ)+ O(∆ζ 2 ).
∆ζ = − 0
p (ζ) + ∆p0 (ζ)
Therefore,
|∆ζ|
α(ζ)
O(∆ζ 2 )
≤
+
0
0
|ζ|
|ζ||p (ζ) + ∆p (ζ)|
and by taking the limit as  → 0 it follows that
κ(ζ, p) ≤α(ζ).
|ζ||p0 (ζ)|

To show that this upper bound is attained, define ∆ai = −sgn(ζ i )|ei |. Then
we have −∆p(ζ) = α(ζ) and it follows that|∆ζ|αζO(∆ζ 2 )
=
+
.
0
0
|ζ|
|ζ||(p (ζ) + ∆p (ζ)|

Taking the limit as  → 0 gives the desired equality.

tuSuppose that ζ > 1 and define ρ = 1/ζ. Then, we compute the condition
number as follows
κ(ζ, p) =αR (ρ),
0|m · pR (ρ) − ρ · pR (ρ)|(16)

Implementation of a Modified Laguerre Method

where αR (ρ) =

11

Pm

i
i=0 |em−i ||ρ| . This result is easily verified by noting thatα(ζ)|ζ|m αR (ρ)
=
0
0
m−1
|z||p (ζ)||ζ||mζ
pR (ρ) − ζ m−2 pR (ρ)||ζ|m αR (ρ)αR (ρ)= m=.
00|ζ| |m · pR (ρ) − ρ · pR (ρ)||m · pR (ρ) − ρ · pR (ρ)|
In summary of the stopping criterion, we compute the backward error
η(ξ) via (12), if ξ ≤ 1, or via (14) otherwise. If η(ξ) ≤ µ, then updatesof the root approximation ξ cease and we compute an approximation of thecondition number κ(ζ, p). This is done by replacing ζ with the approximate
root ξ in (15), if ξ ≤ 1, or in (16) otherwise. Both values are returned andcan be used to determine if the computation of ξ was backward stable andthe approximate sensitivity of ζ to changes in the coefficients of p. Finally, theproduct η(ξ)κ(ζ, p) can be used to provide a first-order bound on the forward
error in the approximate root ξ.

3 Numerical Experiments
In this section, we present the results of several numerical experiments toverify our analysis in Sections 2.1–2.3 and the effectiveness of Algorithm 1.
The results that follow are from tests run on an Intel Core i5 CPU running
at 2.7 GHz with 16GB of memory. All code was compiled using gfortran 8.2.0
with the -O2 optimization flag in the IEEE double-precision floating-point
standard. Fortran 90 source code of Algorithm 1, denoted FPML, and all tests
is available online at https://github.com/trcameron/FPML.
We provide comparisons with two open source Fortran 90 root solvers:
AMVW from [3] and Polzeros from [4]. We are interested in comparing against
AMVW as it is a recent advancement in eigenvalue methods for computing
the roots of a polynomial. It has guaranteed backward stability, assuming themethod converges, and quadratic time complexity. Note that we are comparing
against the most recent single-shift version of AMVW which is maintained
at https://github.com/eiscor/eiscor.
We are interested in comparing against Polzeros as it strongly motivatedthe development of Algorithm 1 from the initial approximations to the stopping
criterion. Furthermore, Polzeros is driven by the Aberth-Ehrlich-Börsh-Supanmethod, which is similar to the modified Laguerre method that we employ.
We recognize that Polzeros has been superseded by MPSolve [5, 6] but bothmake use of the Aberth-Ehrlich-Börsch-Supan method. Furthermore, MPSolveis an adaptive multi-precision software written in C, and as of version 3.0, it is
implemented in parallel [6]. In the future, we are interested in writing sourcecode for Algorithm 2 in C and applying other features such as adaptivity andmulti-precision at which point we will compare against MPSolve. For now,
we are comparing against a Fortran 90 version of Polzeros, which is available
at https://jblevins.org/mirror/amiller/.

12

Thomas R. Cameron

Table 1 Convergence Rate Tests

Iteration

Error-1

Error-2Error-3Error-4Error-5Error-6

1
2
3
4
5
6
7
8
9

0.71
0.4
1.06
2.45 · 10−4
8.07 · 10−15
2.7 · 10−17
0
0
0

1.75
1.27
0.11
9.98 · 10−2
0.94
1.13 · 10−96.84 · 10−15
0
0

1.76
0.42
0.84
0.4
9.96 · 10−2
2.24 · 10−4
3.41 · 10−16
0
0

0.83
0.65
0.72
0.38
0.78
3.95 · 10−2
2.96 · 10−3
2.94 · 10−15
0

7.13
16.6
0.75
0.23
0.34
8.4 · 10−4
1.08 · 10−13
0
0

3.47
4.74
0.75
0.82
0.16
0.52
9.68 · 10−9
5.18 · 10−140

We also ran tests on the C02AFF routine from the NAG library [28]. For
large degree polynomials, the performance of C02AFF is not comparable withthe other root solvers tested. For this reason, we do not include the results of
C02AFF in the CPU benchmark tests. However, on small degree polynomials,
the C02AFF routine performs quite well, so we include the results of C02AFF
in the accuracy benchmark tests.
In Section 3.1, we demonstrate that the fourth-order convergence of Algorithm 1 can be observed in practice. In Section 3.2, we provide several benchmark tests for comparing the elapsed CPU time of FPML with AMVW and
Polzeros when solving an array of polynomial equations. Finally, in Section 3.3,
we use accuracy benchmark tests that rely on special polynomials to test for
convergence difficulties and the deflation stability of FPML, AMVW, Polzeros,
and C02AFF.

3.1 Convergence Rate Tests
Six polynomials are used to test the local convergence rate of FPML. Eachpolynomial has random roots in the unit circle of the complex plane. After
generating these random roots, we use the rootstocoeffs function from thetesting in [3] to form the coefficients of the corresponding polynomial. Thisfunction relies on multi-precision software for Fortran [27].
Starting with the initial approximations outlined in Section 2.1, the error after each iteration is recorded as the maximum relative forward error in
our root approximations. The results are recorded in Table 1, where the column Error-i corresponds to the error in the root approximations for the ith
polynomial; the ith polynomial has degree 4 + 4i, for i = 0, 1, . . . , 5.
Note that fourth-order convergence (or better) can be observed in eachcolumn of Table 1, which indicates a surprising feature of the modified Laguerre method from [33]. Often, higher order methods do not display theirconvergence rate in practice. For instance, the method in [17] has fourth-orderconvergence, but the approximations must be so close to the roots before thisconvergence rate will set in that it will rarely be noticed in double-precision
floating-point arithmetic.

Implementation of a Modified Laguerre Method

13

3.2 CPU Time Benchmark Tests
Our CPU benchmark tests are split into three categories of increasing difficulty: The first category consists of polynomials with random coefficients andthe roots of unity. The second category consists of polynomials with wellconditioned roots along a circular curve in the complex plane. The third category consists of a polynomial with random roots in the unit circle and thetruncated exponential.
3.2.1 Category One
Each random polynomial has coefficients with real and imaginary parts uniformly distributed on the interval [−1, 1]. Figure 1 includes the elapsed timemeasured in seconds and the error measured as the average maximum backward error over all trials. The backward error is computed via (12), where ei
is defined as in (13). Iterations are run for polynomials of degree 80 to degree 20480, doubling the degree on each step. The number of trials performed
ranges from 512 to 2, cutting in half as the degree doubles.
2πi
The roots of the polynomial z n − 1 are the n roots of unity: e n j for
j = 1, . . . , n. Figure 1 includes the elapsed time measured in seconds and theerror measured as the average of the maximum relative forward error over all
trials. Again, iterations are run for polynomials of degree 80 to degree 20480,
doubling the degree on each step; the number of trials ranges from 512 to 2,
cutting in half as the degree doubles.
Note that the backward error of FPML and Polzeros is nearly identical;
similarly, the forward error is the same except for a slight spike from FPMLat degree 10240. In addition, the elapsed time verifies our time complexity
analysis of Algorithm 1, see Section 2.2, since both AMVW and Polzeros
are known to have a quadratic time complexity. Furthermore, the fact that
FPML is faster than Polzeros indicates that the fourth-order convergence ofthe modified Laguerre method from [33] is making up for its additional cost periteration over the Aberth-Ehrlich-Börsh-Supan method. Finally, we note thatthe growth in the backward error of AMVW is undesirable as it shows thatthe resulting approximations ξ are roots of a polynomial whose coefficients are
more perturbed than the worst case scenario for fl(p(ξ)), see Section 2.3.
3.2.2 Category TwoPn
Pn
1
Two polynomials are tested: p1 (z) = i=0 (i + 1)z i and p2 (z) = i=0 (i+1)
zi.
Both polynomials have well-conditioned roots that lie on a circular curve inthe complex plane. Figure 2 includes the elapsed time measured in seconds andthe error measured as the average maximum backward error over all trials. The
backward error is computed via (12), where ei is defined as in (13). Iterations
are run for polynomials of degree 80 to degree 10240, doubling the degree on
each step. The number of trials performed ranges from 512 to 4, cutting inhalf as the degree doubles.

14

Thomas R. Cameron
Roots of Unity

Random Coefficients

10

1

10

0

FPML
Polzeros
AMVW

102
10
Time (s)Time (s)

102

10−1

1

FPML
Polzeros
AMVW

100
10−1
10−2

10−2

10−3

10−3

10−4
10−13

Forward Error

Backward Error

10−13
10−14

10−15

10−14

10−15
10−16
102

103
Degree

104

102

103
Degree104

Fig. 1 Category One Testing

Note that the backward error for FPML, Polzeros, and AMVW range
between 10−16 and 10−15 when testing on p1 (z). However, for the testing on
p2 (z), there is a significant spike from Polzeros that indicates the failure toconverge for at least one root approximation. This somewhat surprising failurefrom Polzeros persisted even after increasing the number of iterations.

3.2.3 Category Three
Two polynomials are tested: The first has random roots in the unit circle of thecomplex plane; these polynomials arePcreated as outlined in Section 3.1. Themsecond is the truncated exponential: i=0 xi /i!. The truncated exponential iswell-known to have ill-conditioned roots, especially for larger m. Also, randomroots in the unit circle can be ill-conditioned as a result of clustering. In
Figure 3, we test on polynomials of degree 10 to degree 100, increasing thedegree by 2 on each step. The number of trials is 512 for each degree, andthe error is measured as the average maximum backward error over all trials.
Again, the backward error is computed via (12), where ei is defined as in (13).
Note that the backward error of FPML and Polzeros are both within unit
roundoff. However, the backward error of AMVW ranges between 10−16 and
100 as the degree increases from 10 to 100.

Implementation of a Modified Laguerre Method

15
Testing on p2 (z)Testing on p1 (z)FPMLPolzeros
AMVW

101
100

100
Time (s)Time (s)FPML
Polzeros
AMVW

101

10−1

10−1

10−2

10−2

10−3

10−3

10−15

Backward Error

Backward Error

10−1

10−16

10−5
10−9
10−13
10−17

102

103
Degree

104

102

103
Degree104

Fig. 2 Category Two Testing
Truncated Exponential

Random Roots in Unit Circle
FPML
Polzeros
AMVW

FPML
Polzeros
AMVW
10−3
Time (s)Time (s)

10−3

10−4

10−4

10−5
10−3
Backward Error

Backward Error

10−1
10−5
10

−9

10−13

10−7
10−11
10−15

10−17
101

102
DegreeFig. 3 Category Three Testing

10−19

101

102
Degree16

Thomas R. Cameron

Table 2 Special polynomials tested by Chandrasekaran et al in [9]Poly No.

Description

Deg.

Roots

1
2
3
4
5
6
7
8
9
10

Wilkinson polynomial
Wilkinson polynomial
Wilkinson polynomial
scale and shifted Wilkinson polynomial
reverse Wilkinson polynomial
reverse Wilkinson polynomial
reverse Wilkinson polynomialprescribed roots of varying scaleprescribed roots of varying scale -3
Chebyshev polynomial

10
15
20
20
10
15
20
20
20
20

1, . . . , 10
1, . . . , 15
1, . . . , 20
−2.1, −1.9, . . . , 1.7
1, 1/2, . . . , 1/10
1, 1/2, . . . , 1/15
1, 1/2, . . . , 1/20
2−10 , 2−9 , . . . , 29
2−10 − 3, . . . , 29 − 3π)cos( 2j−140

11

z 20 + z 19 + · · · + z + 120

ei 21 π2jTable 3 Special polynomials used to test MPSolve in [5]Poly No.

Description

Deg.

Roots

12
13
14

C. Traverso
Mandelbrot
Mandelbrot

24
31
63

known
known
known

3.3 Accuracy Benchmark Tests
The polynomials in Table 2 are tested in [3, 9]. Many of these polynomials are
infamous for their ill-conditioned roots and are suitable for testing convergencedifficulties, e.g., Wilkinson and Chebyshev polynomial.
Polynomials taken from MPSolve [5] and tested in [3] are given in Table 3.
We take the roots computed by MPSolve to be exact. The polynomial provided by Carlo Traverso arises from the symbolic processing of a system ofpolynomial equations and has multiple roots. The Mandelbrot polynomials
are defined recursively as follows: p0 (z) = 1 and pj (z) = zpj−1 (z)2 + 1, forj = 1, 2, . . . , k, with m = 2k − 1.
The polynomials in Table 4 are from Jenkins and Traub [19] and are designed to test convergence difficulties and the deflation stability of a root
finding method. The cubic polynomial p1 (z) = (z − a)(z + a)(z − 1) examines
the stopping criterion. Underflow is likely in the polynomial p3 (z) =
Qm
−j
(z
−10
). Thus, this polynomial may be used to test whether a stopping
j=1
criterion based on roundoff error fails. To test the performance of our methodon multiple or nearly multiple zeros we use the polynomials
p4 (z) = (z − .1)3 (z − .6)(z − .7),

p5 (z) = (z − .1)4 (z − .2)3 (z − .3)2 (z − .4),

p6 (z) = (z − .1)(z − 1.001)(z − .998)(z − 1.00002)(z − .99999),

p7 (z) = (z − .001)(z − .01)(z − .1)(z − .1 + ia)(z − .1 − ia)(z − 1)(z − 10),p8 (z) = (z + 1)5 .

Implementation of a Modified Laguerre Method

17

Table 4 Special Polynomials to test root finding algorithms from [19]Poly No.

Description

Deg.

Roots

15
16
17
18
19
20
21
22
23
24
25
26
27

p1 (z) with a = 10−8
p1 (z) with a = 10−15
p1 (z) with a = 108
p1 (z) with a = 1015
p3 (z)
p3 (z)
p4 (z)
p5 (z)
p6 (z)
p7 (z) with a = 0
p7 (z) with a = 10−10
p7 (z) with a = 10−6
p8 (z)

3
3
3
3
10
20
6
10
5
7
7
7
5

a, −a, 1a, −a, 1a, −a, 1a, −a, 1
10−1 , . . . , 10−10
10−1 , . . . , 10−20
1/10 (m. 3), 5/10, 6/10, 7/10
1/10 (m. 4), 2/10 (m. 3), 3/10 (m. 2), 4/10
0.1, 0.998, 1.00002, 0.99999
10−3 , 10−2 , 10−1 (m. 2), 1, 10
10−3 , 10−2 , 10−1 − ia, 10−1 + ia, 1, 10
10−3 , 10−2 , 10−1 − ia, 10−1 + ia, 1, 10−1 (m. 5)

28
29
30
31

p9 (z)
p10 (z) with a = 103
p10 (z) with a = 106
p10 (z) with a = 109

20
3
3
3

102 e 10 j , 10−2 e 10 j
a, 1, 1/aa, 1, 1/aa, 1, 1/a32

p11 (z) with m = 1560

e 2m j , 0.9e 2m j

2πi

2πi

πi

πi

33

p11 (z) with m = 20

80

e

34

p11 (z) with m = 25

100

e

πi j
2m
πi j
2m

πi

, 0.9e 2m j
πi

, 0.9e 2m j

The last convergence test is on the polynomial p9 (z) = (z 10 −10−20 )(z 10 +1020 )whose equimodular roots are known to cause difficulties for algorithms thatuse powering techniques to separate zeros. Finally, we are interested in testing
the deflation stability of our method. To this we end, we use the polynomials

p10 (z) = (z − a)(z − 1)(z − a−1 ),p11 (z) =

m−1
Y

j=1−m

πi

(z − e 2m j )

3m
Y

πi

(z − .9e 2m j ).

j=m

The maximum relative forward error from FPML, AMVW, Polzeros, andC02AFF when solving for the roots of all polynomials listed in Tables 2–
4 is recorded in Table 5. It is clear from Table 5 that FPML and Polzeros
have comparable accuracy. This is not surprising since the deflation strategies of both these algorithms are similar and the stopping criterion is thesame. It is interesting to observe that there are several tests where the performance of AMVW is significantly worse than FPML, Polzeros, and C02AFF(e.g., Poly No. 6, 12, 16, 19, 20, 28). Finally, we note that there are several tests
where C02AFF is superior to FPML, AMVW, and Polzeros (e.g., Poly No.
1, 2, 3, 23). The superior performance of C02AFF is the effect of a sharper
stopping criterion.

18

Thomas R. Cameron

Table 5 Accuracy Benchmark Tests

Poly No.

FPML

Polzeros

AMVW

C02AFF

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

3.04 · 10−10
3.3 · 10−5
8.8 · 10−2
6.3 · 10−11
4.44 · 10−9
2.09 · 10−5
8.07 · 10−2
4.8 · 10−14
6.27 · 10−2
1.61 · 10−98.28 · 10−15
1.73 · 10−7
1.37 · 10−4
0.29
1.65 · 10−16
1.97 · 10−16
3.14 · 10−16
2.48 · 10−16
4.24 · 10−16
4.56 · 10−15
2.86 · 10−5
1.58 · 10−3
1.49 · 10−4
2.82 · 10−5
2.82 · 10−5
1.93 · 10−5
2.35 · 10−3
2.54 · 10−15
5.17 · 10−23
9.69 · 10−16
1.38 · 10−15
1.89 · 10−7
4.94 · 10−7
7.04 · 10−7

3.38 · 10−9
5.01 · 10−5
6.01 · 10−2
1.62 · 10−12
5.29 · 10−10
2.28 · 10−6
6.32 · 10−2
5.33 · 10−15
6.01 · 10−2
3.54 · 10−10
1 · 10−15
4 · 10−7
9.17 · 10−6
0.3
3.03 · 10−19
1.04 · 10−15
2.65 · 10−16
2.65 · 10−16
4.24 · 10−16
3.83 · 10−15
2.68 · 10−5
1.65 · 10−3
6.08 · 10−5
2.92 · 10−5
2.92 · 10−5
2.29 · 10−5
2.11 · 10−3
1.59 · 10−16
4.34 · 10−22
1.16 · 10−16
1.19 · 10−16
2.01 · 10−7
5.42 · 10−7
4.57 · 10−7

3.12 · 10−10
6.26 · 10−6
6.48
3.86 · 10−12
2.16 · 10−6
0.44
2.01
0.46
8.91 · 10−2
2.05 · 10−10
1.68 · 10−15
115
4.31 · 10−7
0.24
5.54 · 10−9
0.2
5 · 10−93.75 · 10−16
1.74 · 107
3.04 · 1016
5.98 · 10−5
1.4 · 10−2
2.52 · 10−5
1.41 · 10−4
1.41 · 10−4
8.32 · 10−5
1.19 · 10−3
1.01
2.27 · 10−16
2.33 · 10−16
4.14 · 10−16
5.87 · 10−8
1.98 · 10−7
2.97 · 10−7

4.25 · 10−11
4.39 · 10−8
4.44 · 10−4
1.17 · 10−12
1.03 · 10−10
9.88 · 10−7
1.97 · 10−3
1.15 · 10−14
2.93 · 10−2
2.4 · 10−11
2.22 · 10−15
2.9 · 10−7
1.02 · 10−70.68
01.97 · 10−16
1.49 · 10−16
1.25 · 10−16
4.24 · 10−16
6.35 · 10−16
2.05 · 10−5
1.35 · 10−3
4.46 · 10−7
1.6 · 10−5
1.6 · 10−5
1.28 · 10−503.29 · 10−150001.14 · 10−7
7.23 · 10−7
1.36 · 10−6

4 Conclusion
The modified Laguerre method from [33] has strong virtues including local fourth-order convergence and a workload that is well-suited for dataparallelism. The former virtue is proved in [33, Section 4] and we demonstrate
the latter in Algorithm 2.

Implementation of a Modified Laguerre Method

19

We derived an instance of the modified Laguerre method in 2 and used
it to form an algorithm for the approximation of all roots of a polynomial,
which we outline in Algorithm 1. Furthermore, we provide a detailed analysis
of our algorithm’s initial approximations, stopping criterion, and stability in
Sections 2.1–2.3.
In Section 3, numerical experiments are provided to verify our analysis andthe effectiveness of our algorithm for solving large degree polynomial equations.
The experiments in Section 3.1 demonstrate that the fourth-order convergenceof Algorithm 1 can be observed in practice. The experiments in Section 3.2
verify the quadratic complexity of our algorithm and highlight its superior
speed in comparison to the root solvers AMVW from [3] and Polzeros from [4].
The experiments in Section 3.3 indicate the reliability of our algorithm evenin the face of difficult polynomial equations.
Finally, we note that all computations from FPML in Section 3 were backward stable as the stopping criterion in (13) was satisfied by all root approximations. Future research includes the implementation of Algorithm 2 alongwith other features such as adaptivity and multi-precision.
Acknowledgements We are grateful to the NAG for providing a free trial that was usedin the testing of our work. In addition, we are thankful for numerous private conversations
with Dario Bini and David Watkins throughout our research and for the referees’ comments
and suggestions that significantly improved this article.

References
1. Aberth, O.: Iteration methods for finding all zeros of a polynomial simultaneously. Math.
Comp. 27(122), 339–344 (1973)
2. Andrew, A.M.: Another efficient algorithm for convex hulls in two dimensions. Info.
Proc. Letters 9(15), 216–219 (1979)
3. Aurentz, J.L., Mach, T., Vandebril, R., Watkins, D.S.: Fast and backward stable computation of roots of polynomials. SIAM J. Matrix Anal. Appl. 36(3), 942–973 (2015)
4. Bini, D.A.: Numerical computation of polynomial zeros by means of Aberth’s method.
Numer. Algorithms 13, 179–200 (1996)
5. Bini, D.A., Fiorentino, G.: Design, analysis, and implementation of a multiprecision
polynomial root finder. Numer. Algorithms 23, 127–173 (2000)
6. Bini, D.A., Robol, L.: Solving secular and polynomial equations: A multiprecision algorithm. J. Comput. Appl. Math. 272, 276–292 (2014)
7. Börsch-Supan, W.: A posteriori error bounds for the zeros of polynomials. Numer.
Math. 5, 380–398 (1963)
8. Chaitin-Chatelin, F., Frayssé, V.: Lectures on Finite Precision Computations. Software,
Environments and Tools. SIAM, Philadelphia, PA (1996)
9. Chandrasekaran, S., Gu, M., Xia, J., Zhu, J.: A fast QR algorithm for companion
matrices. In: Recent Advances in Matrix and Operator Theory, pp. 111–143. Birkhäuser
Basel, Basel, Switzerland (2008)
10. Chen, T.C.: Improving Laguerre’s method to cope with symmetry pitfalls in polynomialroot-finding. In: Proceedings of the Eigth International Colloquium on Differential
Equations, pp. 105–110. VSP, Utrecht, Netherlands (1998)
11. Ehrlich, L.W.: A modified Newton method for polynomials. Commun. ACM 10(2),107–108 (1967)
12. Foster, L.V.: Generalizations of Laguerre’s method: Higher order methods. SIAM J.
Numer. Anal. 18(6), 1004–1018 (1981)

20

Thomas R. Cameron

13. Goedecker, S.: Remarks on algorithms to find roots of polynomials. SIAM J. Sci. Comput. 15(5), 1059–1063 (1994)
14. Goldberg, D.: What every computer scientist should know about floating-point arithmetic. ACM Comput. Surv. 23, 5–48 (1991)
15. Grad, J., Zakrajšek, E.: LR algorithm with Laguerre shift for symmetric tridiagonal
matrices. Comput. J. 15(3), 268–270 (1972)
16. Hansen, E., Patrick, M.: A family of root finding methods. Numer. Math. 27, 257–269(1977)
17. Hansen, E., Patrick, M., Rusnak, J.: Some modifications of Laguerre’s method. BIT
17(4), 409–417 (1977)
18. Higham, N.J.: Accuracy and Stability of Numerical Algorithms. SIAM, Philadelphia,
PA (2002)
19. Jenkins, M.A., Traub, J.F.: Principles for testing polynomial zerofinding programs.
ACM Trans. Math. Software 1(1), 26–34 (1975)
20. Kahan, W.: Laguerre’s method and a circle which contains at least one zero of a polynomial. SIAM J. Numer. Anal. 4(3) (1967)
21. Laguerre, E.N.: Oeuvres de Laguerre. Gauthier-Villars and Fils, Paris (1898)
22. Lancaster, P.: Lambda-matrices and Vibrating Systems, International Series of Monographs on Pure and Applied Mathematics, vol. 94. Pergamon, Oxford, United Kingdom(1966)
23. Lenard, C.T.: Laguerre’s iteration and the method of traces for eigenproblems. Appl.
Math. Lett. 3(3), 73–74 (1990)
24. Leuze, M.R.: A hybrid Laguerre method. BIT 23(1), 132–138 (1983)
25. Maehly, V.H.: Zur iterativen auflösung algebraisher gleichungen. Z. Angew. Math. Phys.
5, 260–263 (1954)
26. Mekwi, W.: Iterative Methods for Roots of Polynomials. Master’s thesis, University of
Oxford, Oxford, England (2001)
27. MPFUN: Multiprecision Software. http://www.netlib.org/mpfun (2005). Accessed 18
June 2018
28. NAG: Nag Library, Mark 26. https://www.nag.com/numeric/fl/nagdoc_latest/html/
c02/c02aff.html (2017). Accessed 27 May 201829. Pan, V.Y.: Solving a polynomial equation: some history and recent progress. SIAM
Rev. 39(2), 187–220 (1997)
30. Parlett, B.: Laguerre’s method applied to the matrix eigenvalue problem. Math. Comp.
18(87), 464–485 (1964)
31. Pellet, A.E.: Sur un mode de séparation des racines des équations et la formule delagrange. Bull. Sci. Math. 5(2), 393–395 (1881)
32. Petković, M.S., Ilić, S., Rančić, L.: The convergence of a family of parallel zero-finding
methods. Comp. Math. Appl. 48, 455–467 (2004)
33. Petković, M.S., Ilić, S., Tričković, S.: A family of simultaneous zero finding methods.
Comput. Math. Appl. 34(10), 49–59 (1997)
34. Petković, M.S., Petković, L., Živković, D.: Laguerre-like methods for the simultaneous
approximation of polynomial zeros. In: Topics in Numerical Analysis, pp. 189–209.
Springer Vienna, Vienna (2001)
35. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.: Numerical Recipes: The
Art of Scientific Computing, 3rd edn. Cambridge University Press, New York, NewYork (2007)
36. Redish, K.A.: On Laguerre’s method. Int. J. of Math. Educ. in Sci. and Tech. 5(1),91–102 (1974)
37. Sitton, G.A., Burrus, C.S., Fox, J.W., , Treitel, S.: Factoring very-high-degree polynomials. IEEE Sig. Proc. Mag. 20, 27–42 (2003)
38. Smith, B.T.: A Zero Finding Algorithm using Laguerre’s Method. Master’s thesis,
University of Toronto, Toronto, Canada (1967)
39. Walsh, J.L.: On Pellet’s theorem concerning the roots of a polynomial. Ann. of Math.
26, 59–64 (1924)
40. Wilkinson, J.H.: Rounding Errors in Algebraic Processes. Prentice-Hall Inc. (1963)