# Validated Root Enclosures for Interval Polynomials with Multiplicities

Florent Bréhard, Adrien Poteaux, Léo Soudant

► To cite this version:

Florent Bréhard, Adrien Poteaux, Léo Soudant. Validated Root Enclosures for Interval Polynomials with Multiplicities. ISSAC'23: International Symposium on Symbolic and Algebraic Computation 2023, Jul 2023, Tromso, Norway. 10.1145/3597066.3597122. hal-04138791

HAL Id: hal-04138791

https://hal.science/hal-04138791

Submitted on 23 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Validated Root Enclosures for Interval Polynomials with Multiplicities

Florent Bréhard

Univ. Lille, CNRS, Centrale Lille, UMR
9189 CRISTAL
Lille, France

florent.brehard@univ-lille.fr

Adrien Poteaux

Univ. Lille, CNRS, Centrale Lille, UMR
9189 CRISTAL
Lille, France

adrien.poteaux@univ-lille.fr

Léo Soudant

École normale supérieure Paris-Saclay
Gif-sur-Yvette, France

soudant.leo@gmail.com

## ABSTRACT

Twenty years ago, Zeng [28, 30] proposed floating-point algorithms to compute multiple roots of univariate polynomials with real or complex coefficients beyond the so-called “attainable accuracy barrier”. Based on these foundations, we propose a validated numeric point of view on this problem. Our first contribution is an improvement of Zeng’s multiplicity detection algorithm using a simple trick that allows us to recover much higher multiplicities. As our main contribution, we propose two floating-point validated algorithms to compute rigorous enclosures for multiple roots. They consist in carefully combining the ideas underlying Zeng’s numerical algorithms with Newton-like fixed-point validation techniques. We also provide a prototype Julia implementation of these algorithms.

## CCS CONCEPTS

* Mathematics of computing → Interval arithmetic; Nonlinear equations; Solvers.

## KEYWORDS

multiple roots of univariate polynomials, fixed-point validation

# 1 MOTIVATION AND RELATED WORKS

This article is concerned with the numerical computation and validation of multiple roots of univariate polynomials.

*Numerical multiple root-finding.* Polynomial root-finding over real or complex numbers is one of the most common tasks in scientific computing and a major topic in numerical analysis and complexity theory, including nearly-optimal algorithms proposed in the 90s [14, 17]. The complexity of this problem [16] is however subject to the so-called “attainable accuracy barrier” in presence of multiple roots: generic root-finding algorithms executed on an input polynomial with coefficients truncated to *O*(n) digits cannot retrieve more than *O*(n/d) digits of accuracy for a root of multiplicity *d*. This could require multi-precision floating-point arithmetic, possibly far beyond the standard IEEE 53 bit double precision.

However, early remarks by Kahan [11] and Zeng’s algorithmic contribution [28, 30] showed how to turn this *ill-conditioned* problem into a *well-conditioned* one, provided the polynomial’s multiplicity structure $\ell = (\ell_1, \dots, \ell_k) \in \mathbb{N}_{>0}^k$ is given or can be inferred.

**PROBLEM 1.** Given $p = p_n x^n + p_{n-1} x^{n-1} + \dots + p_0 \in \mathbb{C}_n[x]$ and some maximum tolerated error on the input coefficients, recover the highest possible multiplicity structure $\ell \in \mathbb{N}_{>0}^k$ and compute root approximations $\tilde{z} \in \mathbb{C}^k$ such that $p \approx p_n \prod_{i=1}^k (x - \tilde{z}_i)^{\ell_i}$.

Existing approaches for this problem include Schröder’s modified Newton operator (see [7] and references therein) with quadratic convergence to a multiple root/a cluster of roots; or approximate GCDs (AGCD) of numerical polynomials [4], a notion closely related to multiple root-finding [1, 3, 5, 6, 23–25, 31]. Zeng’s approach [28, 30] that inspires this article combines AGCD computations with Gauss-Newton iterations on overdetermined polynomial systems encoding the roots to overcome this attainable accuracy barrier.

*Validated numerics.* Starting with Moore’s interval arithmetic [13] in the 60s, the field of validated numerics [15, 21, 27] aims at replacing heuristic or asymptotic error bounds with rigorous, explicitly computed enclosures. It relies on set-valued representations of mathematical objects, together with correctness-preserving algorithmic operations. For example, a real number $x$ is represented by an interval $[x] = [x, \bar{x}]$ with floating-point endpoints, i.e., $x \in [x]$. Arithmetic operations are defined on intervals to preserve this relation: if $x \in [x]$ and $y \in [y]$, then $[x] \star [y]$ for $\star \in \{+, -, \times, \div\}$ is built so that $x \star y \in [x] \star [y]$ holds. In particular, rounding errors are taken into account. Similar representations exist for complex numbers (which we will still refer to as intervals throughout this article), e.g., a pair of real intervals for the real and imaginary parts.

Most validated numerics libraries (e.g., INTLAB [19] or Arb [10]) can compute tight enclosures for simple roots of univariate polynomials. For multiple roots, a precise statement is needed to determine with respect to what the root enclosures are validated.

**PROBLEM 2.** Given $[p] = [p_n]x^n + \dots + [p_0]$ with interval coefficients, $0 \notin [p_n]$ and $\ell \in \mathbb{N}_{>0}^k$, compute intervals $[z] = ([z_1], \dots, [z_k])$ such that for any polynomial $p \in [p]$ having multiplicity structure $\ell$, there exist unique $z_i \in [z_i]$ such that $p = p_n \prod_{i=1}^k (x - z_i)^{\ell_i}$.

The rigorous root enclosure methods of [20] are subject to the attainable accuracy barrier, since they are correct w.r.t. all polynomials $p \in [p]$, including those having clusters of roots.

Methods based on deflation [12, 22] consist in making an overdetermined system invertible by adding extra variables for perturbations on the coefficients of $p$. For instance, for a root $z$ of multiplicity $k$, consider $\hat{p}(x) = p(x) + \varepsilon_{k-2}x^{k-2} + \dots + \varepsilon_0$ and solve the system $\{\hat{p}^{(j)}(z) = 0, 0 \le j < k\}$ for the variables $z, \varepsilon_0, \dots, \varepsilon_{k-2}$. The computed enclosure $[z]$ is guaranteed to contain a root of multiplicity $k$ of an $\varepsilon$-perturbation of $p$. Such a property, however, is not per se a solution to Problem 2: the correcting perturbation $\hat{p} - p$ is not guaranteed to be the same for all distinct roots, the resulting $\hat{p}$ does not necessarily belong to $[p]$ and its coefficients could be close to those of $p \in [p]$ while having completely different roots.
**Our contributions.** They are twofold. First, after recalling the main lines of Zeng's two numerical algorithms GCDROOT and PEJROOT in the beginning of Section 2, we improve the former (for multiplicity detection) using a simple trick. This gives Algorithm GCDROOT* in Section 2.3 which can correctly recover much higher multiplicities in our numerical experiments.

As our main contribution in Section 3, we propose an a posteriori validation method to simultaneously compute rigorous enclosures for the multiple roots, by carefully combining the ideas underlying Zeng's numerical algorithms with Newton-like fixed-point validation techniques. In particular, the combined nonlinear and overdetermined nature of Problem 2 is an obstacle for classical fixed-point validation. To overcome this, we design an original two-phase strategy based on two complementary validated algorithms VALGCDROOT and VALPEJROOT. Our focus is on validating multiple roots beyond the attainable accuracy barrier: all our examples were executed in machine double precision to highlight this point.

Our implementation ValidatedMultipleRoots.jl<fnref n="1" /> is a prototype open-source package in Julia [2]. It is based on Verzani's implementation of Zeng's algorithms in the Polynomials.jl package<fnref n="2" />. Although Zeng's original implementation in Matlab [29] is still available<fnref n="3" />, we found more convenient to rely on this recent Julia code, notably for multi-precision and interval computations.

**A note about complexity.** We do not focus on the complexity of our algorithm: a complete complexity analysis would require a priori bounds for the precision needed to solve our problem, which would go beyond the scope and length of this paper. Nevertheless, if we stick to finite precision computations and only use classical matrix multiplication, it is easy to check that our algorithms run in complexity $O(k^2 n)$, where $n$ is the input degree and $k$ the number of distinct roots. For instance, computations of the example shown in Table 1 are almost linear in $n$ ($k = 4$ therein is constant).

## 2 ZENG'S ALGORITHM

Zeng's method to compute multiple roots of inexact polynomials [28, 30] mainly consists in numerically solving overdetermined systems of linear and nonlinear equations, in some least-squares sense. The necessary background is recalled in Section 2.1. Then we summarise the main ideas of the two algorithms GCDROOT and PEJROOT composing Zeng's method in Section 2.2. After that, in Section 2.3, we present our first contribution, Algorithm GCDROOT*: a modification of GCDROOT that overcomes some of its limitations by reducing the impact of errors in the input.

### 2.1 Preliminaries.

The following paragraphs summarise notions such as pejorative manifolds, Sylvester matrices and approximate GCD. They make it possible to talk about multiple roots in a numerical context with precise mathematical definitions.

**Pejorative manifolds.** Given a multiplicity structure $\ell \in \mathbb{N}_{>0}^k$ with $\ell_1 + \cdots + \ell_k = n$, the function $G_\ell$ maps the roots $z \in \mathbb{C}^k$ to the

coefficients of the polynomial $p = \prod_{i=1}^k (x - z_i)^{\ell_i}$, i.e.,

$$G_{\ell} : z \in \mathbb{C}^{k} \mapsto \begin{pmatrix} g_{0}(z) \\ \vdots \\ g_{n-1}(z) \end{pmatrix} \in \mathbb{C}^{n},$$

$$\text{s.t.} \quad p = x^n + g_{n-1}(z)x^{n-1} + \cdots + g_0(z).$$

**Definition 2.1.** The pejorative manifold $\Pi_\ell$ is the image of the map $G_\ell : \mathbb{C}^k \to \mathbb{C}^n$. It is a manifold (with singularities) of the coefficient space $\mathbb{C}^n$ parameterised by the roots $(z_1, \dots, z_k) \in \mathbb{C}^k$.

The tangent space of $\Pi_\ell$ at $z \in \mathbb{C}^k$ is given by the Jacobian:

$$J_{G_{\ell}}(z) = \left( \begin{array}{c|c} q_1 & \cdots \\ \end{array} \right) \in \mathbb{C}^{n \times k}, \\ \text{where } q_i := q_{i,n-1}x^{n-1} + \cdots + q_{i,0} := -\frac{\ell_i}{x-z_i}p. \quad (1)$$

The singularities of $\Pi_\ell$ are the points where $J_{G_\ell}(z)$ is rank deficient. This happens when at least two roots $z_i$ and $z_j$ coincide: the corresponding point belongs to the pejorative submanifold $\Pi_{\ell'} \subseteq \Pi_\ell$ associated to a higher multiplicity structure $\ell' \in \mathbb{N}_{>0}^{k'}$ with $k' < k$.

As it will be seen in Section 2.2, Algorithm PEJROOT of [28] relies on Gauss-Newton iterations to converge to a solution of $G_\ell(z) = p$ on the pejorative manifold $\Pi_\ell$.

**Convolution and Sylvester matrices.** Notions like divisibility and GCD in $\mathbb{C}[x]$ may seem ill-posed in a numerical context with errors or uncertain data. Yet, using convolution and Sylvester matrices, these questions are reduced to classical problems in linear algebra. They can be solved using efficient algorithms in numerical linear algebra, such as the singular value decomposition (SVD) [8, §2.5].

Let $p \in \mathbb{C}_n[x]$. For a given $k \ge 0$, the multiplication by $p$ induces a morphism $\mathbb{C}_k[x] \to \mathbb{C}_{n+k}[x]$, the representation of which in the canonical bases is given by the $k$-th convolution matrix $C_k(p)$:

$$C_k(p) = \begin{pmatrix} p_0 \\ p_1 \\ \vdots \\ p_n \end{pmatrix} \begin{pmatrix} p_0 \\ p_1 \\ \vdots \\ p_n \end{pmatrix} \in \mathbb{C}^{(n+k+1) \times (k+1)}.$$

For two polynomials $p \in \mathbb{C}_n[x]$ and $q \in \mathbb{C}_m[x]$ of respective coefficients $p \in \mathbb{C}^{n+1}$ and $q \in \mathbb{C}^{m+1}$, the coefficients of $pq$ are given by the vector $C_m(p)q = C_n(q)p$.

**Definition 2.2.** If $p \in \mathbb{C}_n[x]$, $q \in \mathbb{C}_m[x]$, and $0 \le d \le \min(n, m)$, then the $d$-th Sylvester matrix of $(p, q)$ is

$$S_d(p, q) = \left( C_{m-d}(p) \mid C_{n-d}(q) \right) \in \mathbb{C}^{(n+m-d+1) \times (n+m-2d+2)}.$$

$S_d(p, q)$ represents the map

$$(w, v) \in \mathbb{C}_{m-d}[x] \times \mathbb{C}_{n-d}[x] \mapsto wp + vq \in \mathbb{C}_{n+m-d}[x].$$

These matrices are connected with the GCD of $p$ and $q$.

**LEMMA 2.3 ([31, PROP. 4]).** Let $p \in \mathbb{C}_n[x]$, $q \in \mathbb{C}_m[x]$, and $u$ be their GCD. Denote by $d$ the degree of $u$. Then $S_d(p, q)$ has full (column) rank for $d < i \le \min(n, m)$, whereas $S_d(p, q)$ has a kernel of dimension 1 spanned by the vector $(w, -v)$ corresponding to cofactors $w \in \mathbb{C}_{m-d}[x]$ and $v \in \mathbb{C}_{n-d}[x]$ such that $p = uv$ and $q = uw$.

<footnote n="1">¹https://gitlab.univ-lille.fr/florent.brehard/ValidatedMultipleRoots.jl/ - v. 0.0.1</footnote>

<footnote n="2">²https://juliamath.github.io/Polynomials.jl/</footnote>

<footnote n="3">³see "Supplemental Material" at https://dl.acm.org/doi/10.1145/992200.992209</footnote>
The numerical treatment of Sylvester matrices underlies the key concept of approximate GCD.

Approximate GCD (AGCD). Computing the exact GCD of two polynomials p and q is an ill-posed problem, since under generic perturbations of the input, the result is always 1. The notion of approximate GCD [5] [4, §1] offers a robust definition by asking for the exact GCD of perturbed polynomials p + δp and q + δq.

Definition 2.4. Let $p \in \mathbb{C}_n[x]$, $q \in \mathbb{C}_m[x]$ and $\varepsilon \ge 0$. We say that the polynomial $u$ of degree $d \le \min(n, m)$ is an $\varepsilon$-GCD for $(p, q)$ if:

(i) $u$ is the exact GCD of a pair $(p+\delta p, q+\delta q)$ with $\delta p \in \mathbb{C}_n[x]$, $\delta q \in \mathbb{C}_m[x]$, and $\|\delta p\|_2^2 + \|\delta q\|_2^2 \le \varepsilon^2$;

(ii) No polynomial of degree larger than $d$ satisfies condition (i);

(iii) $u$ minimises $\|\delta p\|_2^2 + \|\delta q\|_2^2$ among degree-$d$ polynomials satisfying (ii).

The method proposed in [5] relies on the observation that the GCD of $p + \delta p$ and $q + \delta q$ has degree at least $d$ if and only if $S_d(p+\delta p, q+\delta q)$ is rank-deficient. Calling $\sigma_{-1}$ the smallest singular value of $S_d(p, q)$ computed by SVD, they moreover prove that such an $\varepsilon$-perturbation $(\delta p, \delta q)$ exists only if $\sigma_{-1} \le \varepsilon\sqrt{\max(n, m) - d + 1}$. Combining this necessary condition with a quadratic optimisation strategy to refine the AGCD $u$ together with the cofactors $v$ and $w$, the authors of [31] proposed a fully automated algorithm AGCD that computes the $\varepsilon$-GCD of $p$ and $q$. This routine is a key ingredient of Algorithm GCDROOT of [28] summarised in Section 2.2 below.

## 2.2 Zeng's algorithm

Zeng's method is based on two algorithms called GCDROOT and PEJROOT. GCDROOT infers the multiplicity structure of the input polynomial up to some user-defined tolerance, and computes root approximations with moderate accuracy. These roots serve as an initialiser for PEJROOT which refines them to high accuracy.

### 2.2.1 Algorithm GCDROOT
It infers the multiplicity structure of $p \in \mathbb{C}_n[x]$ (together with root approximations) by inspecting the AGCD, defined in the previous section, of $p$ and $p'$.

We now specialise Lemma 2.3 to the case $(p, q) = (p, p')$. For $1 \le k \le n$ the $k$-th discriminant matrix $\Delta_k(p)$ is:

$$ \Delta_k(p) := S_{n-k}(p, p') = \left( C_{k-1}(p) \mid C_k(p') \right) \in \mathbb{C}^{(n+k) \times (2k+1)}. $$

LEMMA 2.5. $p \in \mathbb{C}_n[x]$ has $k$ distinct roots in $\mathbb{C}$ if and only if all $\Delta_i(p)$ for $1 \le i < k$ have full column rank and $\Delta_k(p)$ has a kernel of dimension 1. This kernel is spanned by the vector $(w, -v)$ corresponding to the coefficients of the cofactors

$$ v = \prod_{i=1}^{k} (x - z_i) \in \mathbb{C}_k[x], \quad w = \sum_{i=1}^{k} \ell_i \prod_{j \ne i} (x - z_j) \in \mathbb{C}_{k-1}[x], \quad (2) $$

where $z_i$ are the roots of $p$ with multiplicity $\ell_i$. The polynomial

$$ u = \prod_{i=1}^{k} (x - z_i)^{\ell_i - 1} \in \mathbb{C}_{n-k}[x], $$

is the GCD of $p$ and $p'$, and we have $uv = p, uw = p'$.

This result is used repeatedly (following somehow the symbolic square-free decomposition algorithm) to define algorithm GCDROOT of [30], that we briefly recall:

**Algorithm 1 GCDROOT(p, ε)**

Input: $p \in \mathbb{C}_n[x]$ and tolerance $\varepsilon$ on the input error

Output: roots $\tilde{z} \in \mathbb{C}^k$ together with multiplicities $\ell \in \mathbb{N}_{>0}^k$

1: Compute the AGCD $u$ of $p$ and $p'$ with cofactors $v, w$

2: Compute the (simple) roots $\tilde{z}$ of $v$

3: Set $u^{(1)} = u$ and for $t = 2, \dots, \min(n, m)$ compute $u^{(t)}$ as the AGCD of $u^{(t-1)}$ and $u^{(t-1)'}$ with cofactors $v^{(t)}, w^{(t)}$, until $\deg(u^{(t)}) = 0$

4: Compute the (simple) roots of all the $v^{(t)}$'s and match them to roots in $\tilde{z}$.

5: Set $\ell_i := \min\{t \text{ s.t. } \tilde{z}_i \text{ is not a root of } v^{(t+1)}\}$

6: Return $\tilde{z}, \ell$

In [30], more details are given on how to propagate the error tolerances in the successive AGCD computations (step 3).

### 2.2.2 Algorithm PEJROOT
The second algorithm proposed by Zeng in [28, 30], PEJROOT, takes as input initial root approximations $\tilde{z}_0$ and the corresponding multiplicity structure $\ell$. Such data is computed by GCDROOT, as previously seen.

Assume $p \in \mathbb{C}_n[x]$ is monic (otherwise rescale it by the inverse of its leading coefficient): $p = x^n + p_{n-1}x^{n-1} + \dots + p_0$ with $p \in \mathbb{C}^n$. $p$ is supposed to lie on the pejorative manifold $\Pi_\ell$, or nearby due to errors in the input coefficients of $p$. Moreover, the initial roots $\tilde{z}_0$ define an initial point $p_0 = G_\ell(\tilde{z}_0)$ on $\Pi_\ell$.

Algorithm PEJROOT consists in applying Gauss-Newton iterations to bring $p_0$ closer and closer to $p$ (or its projection onto $\Pi_\ell$). The underlying system of equations,

$$ G_{\ell}(z) = p \iff \left\{ \begin{array}{ll} g_0(z_1, \dots, z_k) = p_0, & \\ \vdots & \\ g_{n-1}(z_1, \dots, z_k) = p_{n-1}, & \end{array} \right. \quad (3) $$

consists of $n$ polynomial equations in $k$ variables $z_i$. Except in the case where all roots of $p$ are simple, we have $k < n$ and (3) is overdetermined. Its Jacobian $J_{G_\ell}(z)$, Eq. (1), is rectangular with full column rank. Instead of classical Newton iterations, which require a square, invertible Jacobian, PEJROOT uses Gauss-Newton iterations:

$$ \tilde{z}_{t+1} = \tilde{z}_t - J_{G_\ell}(\tilde{z}_t)^+ (G_\ell(\tilde{z}_t) - p), $$

where $J_{G_\ell}(\tilde{z}_t)^+$ is the Moore-Penrose inverse (a.k.a. pseudoinverse) of $J_{G_\ell}(\tilde{z}_t)$. In practice, one merely computes the least-squares solution $h_t \in \mathbb{C}^k$ of the overdetermined linear system:

$$ J_{G_\ell}(\tilde{z}_t) h_t = G_\ell(\tilde{z}_t) - p. $$

Proper convergence results may be found in [30, §3].

## 2.3 First contribution: improving GCDROOT

Larger errors in the coefficients of $p$ of multiplicity structure $\ell \in \mathbb{N}_{>0}^k$ increase the risk of the $\varepsilon$-GCD of $p$ and $p'$ having degree greater
**Algorithm 2 PEJROOT(p, ℓ, ˜z₀)**

**Input:** $p \in \mathbb{C}_n[x]$, $\ell \in \mathbb{N}_{>0}^k$ and initial root approx $\tilde{z}_0 \in \mathbb{C}^k$

**Output:** refined root approximations $\tilde{z} \in \mathbb{C}^k$

For $t = 0, 1, \dots$:

1: Compute the least-squares solution $\mathbf{h}_t \in \mathbb{C}^k$ minimising the residue $d_t = \|J_{G_t}(\tilde{z}_t)\mathbf{h}_t - G_t(\tilde{z}_t) + \mathbf{p}\|_2$.

2: Update the approximation: $\tilde{z}_{t+1} = \tilde{z}_t - \mathbf{h}_t$.

3: If $d_t$ is below some tolerance, return $\tilde{z}_{t+1}$, otherwise continue.

than $n-k$. If so, GCDROOT infers a wrong multiplicity structure $\ell'$, corresponding to a pejorative submanifold $\Pi_{\ell'} \subset \Pi_{\ell}$.

However, the repeated use of AGCD extraction in step 3 of GCD-ROOT is another significant source of errors that can be avoided if one can infer $\ell$ directly, once the $z_i$ have been computed. This is our first contribution: the procedure GCDROOT* based on the following simple observation to deduce the $\ell_i$ from the cofactors $v, w$.

**LEMMA 2.6.** If $z = (z_1, \dots, z_k)$ are the distinct multiple roots of $p \in \mathbb{C}_n[x]$, then the corresponding multiplicities $\ell_i$ are

$$ \ell_i = \frac{w(z_i)}{v'(z_i)}, $$

where $v$ and $w$ are the cofactors as in Lemma 2.5.

The first-order sensitivity of this formula w.r.t. perturbations $\delta v$ and $\delta w$ on the coefficients of $v$ and $w$ is given by:

$$ \delta \ell_i = \frac{v'(z_i)w'(z_i) - v''(z_i)w(z_i)}{v'(z_i)^3} (\delta v)(z_i) - \frac{w(z_i)}{v'(z_i)^2} (\delta w)'(z_i) + \frac{1}{v'(z_i)} (\delta w)(z_i) + O(\|\delta v\|^2 + \|\delta w\|^2). $$

**PROOF.** The formula for $\ell_i$ directly follows from (2) in Lemma 2.5: we have $w(z_i) = \ell_i \prod_{j \neq i} (z_i - z_j)$ and $v'(z_i) = \prod_{j \neq i} (z_i - z_j)$.

Now suppose $v$ and $w$ are perturbed by $\delta v$ and $\delta w$. The resulting perturbations $\delta z_i$ of the roots $z_i$ follow from $(v + \delta v)(z_i + \delta z_i) = 0$:

$$ \delta z_i = \frac{(\delta v)(z_i)}{v'(z_i)} + O(\|\delta v\|^2). $$

By differentiating the relation for $\ell_i$ w.r.t. $v, w, z_i$ and replacing $\delta z_i$, we obtain the desired estimate involving $\delta v$ and $\delta w$. $\square$

In particular, the sensitivity analysis shows that the technique used by Algorithm GCDROOT* to recover the multiplicity $\ell_i$ of the root $z_i$ (line 3) is robust w.r.t. the numerical errors of $v$ and $w$ computed in line 1, except, as expected, if this root is too close to other roots of $p$ (then $v'(z_i)$ tends to 0).

**Numerical example.** We compare our modified routine GCDROOT* to Zeng's original GCDROOT and show how this improves the total numerical multiple-root finding process solving Problem 1.

**Example 1.** Consider, as in [30, §4.6], the polynomial

$$ p_m = (x - 1)^{4m} (x - 2)^{3m} (x - 3)^{2m} (x - 4)^m, \quad (4) $$

for $m = 1, 2, \dots$, with coefficients truncated to standard double precision (53 bits). Table 1 summarises the multiplicity structures recovered by Zeng's GCDROOT and our GCDROOT*. GCDROOT fails as soon as $m = 5$, and the obtained $\ell_i$ rapidly diverge too much

**Algorithm 3 GCDROOT*(p, ε)**

**Input:** $p \in \mathbb{C}_n[x]$ and tolerance $\varepsilon$ on the input error

**Output:** roots $\tilde{z} \in \mathbb{C}^k$ together with multiplicities $\ell \in \mathbb{N}_{>0}^k$

1: Compute the AGCD $u$ of $p$ and $p'$ with cofactors $v, w$

2: Compute the (simple) roots $\tilde{z}$ of $v$

3: Compute multiplicities $\ell_i$ as in Lemma 2.6, and round them to the nearest integer

4: Return $\tilde{z}, \ell$

<table><caption>**Table 1:** Multiplicity structures of $p_m$ (4) recovered by GCD-ROOT and GCDROOT*. Boldface indicates incorrect results. The right part gives the first erroneous cases of GCDROOT*.</caption><thead><tr><th>m</th><th>GCDROOT</th><th>GCDROOT*</th><th>m</th><th>GCDROOT*</th></tr></thead><tbody><tr><td>4</td><td>[16, 12, 8, 4]</td><td>[16, 12, 8, 4]</td><td>49</td><td>[196, 147, 98, 50]</td></tr><tr><td>5</td><td>[22, 15, 9, 4]</td><td>[20, 15, 10, 5]</td><td>56</td><td>[224, 168, 112, 57]</td></tr><tr><td>6</td><td>[37, 11, 7, 5]</td><td>[24, 18, 12, 6]</td><td>59</td><td>[236, 178, 119, 58]</td></tr><tr><td>7</td><td>[47, 8, 10, 5]</td><td>[28, 21, 14, 7]</td><td>63</td><td>[252, 189, 127, 62]</td></tr><tr><td>8</td><td>[56, 8, 11, 5]</td><td>[32, 24, 16, 8]</td><td>64</td><td>[256, 193, 130, 62]</td></tr><tr><td>9</td><td>[68, 7, 10, 5]</td><td>[36, 27, 18, 9]</td><td>65</td><td>[260, 195, 129, 66]</td></tr></tbody></table>

from the actual values. On the other hand, using the same tolerance parameters, GCDROOT* succeeds up to $m = 48$, and the first incorrectly recovered multiplicities (right part of the table) only differ from the actual ones by a few units.

**Remark 2.7.** In Zeng's original Matlab implementation [29], GCD-ROOT on the same example succeeds for $m$ up to 7 and root approximations are more accurate. This is probably due to the fact that the use of rescaling weights in the AGCD procedure was not re-implemented in the more recent Julia package Polynomials. j1 we rely on. It would be interesting to see the effect of weights on both GCDROOT and GCDROOT*.

**About conditioning and numerical stability.** The conditioning of Problem 1 for a given $\ell$, i.e., the perturbation on the roots $z$ induced by a multiplicity-preserving perturbation of the coefficients of $p$ in system (3), is given by the norm of $J_{G_\ell}(z)^+$, which is equal to the smallest singular value of $J_{G_\ell}(z)$ [30, §3.4]. In a floating-point arithmetic context, we may prefer the relative structure-preserving condition number (see, e.g., [26, §12] or [9, §1.6]):

$$ \kappa_{\ell} := \frac{\|J_{G_{\ell}}(z)^+\|}{\|z\|/\|\mathbf{p}\|}. $$

Hence, knowing $\ell$, the problem of computing multiple roots is not subject to infinite sensitivity: only *finitely* many bits are lost. Moreover, experiments in [30, §3.4] show that $\kappa_\ell$ is moderate, even with high multiplicities, if the distinct roots are well separated.

However, rounding errors during execution also affect the computed roots. For PEJROOT, each iteration solves a linear least-squares problem using the QR factorisation, which is backward stable [9, Thm. 20.3] but does not exploit the specific definition of $J_{G_\ell}(z)$. Therefore, rounding errors grow with the matrix condition number $\kappa(J_{G_\ell}(z)) = \|J_{G_\ell}(z)\|\|J_{G_\ell}(z)^+\|$. Similarly, the SVD used by the
AGCD routine in `GCDROOT` and `GCDROOT*` to compute the cofactors $v$ and $w$ is backward stable [26, Thm. 19.4]. Therefore, the rounding errors in $v$ (from which the (simple) roots are extracted afterwards) grow with $\kappa(\Delta_k(p)) = \|\Delta_k(p)\|\|\Delta_k(p)^+\|$.

Although it is claimed in [30] that `PEJROOT` is more accurate than `GCDROOT`, our experiments show it is not a general fact (see Table 3). A deeper numerical analysis of both methods would be necessary to understand in which cases one is better than the other.

## 3 VALIDATING MULTIPLE ROOTS

To solve Problem 2, we propose validation algorithms that compute rigorous error bounds on the numerical roots obtained previously. Algorithm `VALPEJROOT` computes accurate bounds by translating Problem 2 into the straightforward, well-conditioned system (3) already used by `PEJROOT`. Unfortunately, this system combines two difficulties: nonlinearity and overdetermination (see Section 3.3). The consequence is that `VALPEJROOT` must rely on a priori root enclosures to compute tighter rigorous bounds. Such initial enclosures are computed by Algorithm `VALGCDROOT`, detailed in Section 3.2. It relies on a different system to encode the roots, inspired by Zeng's `GCDROOT`, whose structure makes it possible to depend on no a priori root enclosures. Before describing these two algorithms, which together solve Problem 2 (numerical examples are given in Section 3.4), we first provide preliminaries about interval arithmetic and a posteriori validation based on Newton-like fixed-point operators.

### 3.1 Preliminaries for the validation strategy.

Interval arithmetic. We call $\mathbb{R}$ and $\mathbb{C}$ the sets of real or complex intervals representable using floating-point arithmetic. Intervals are denoted using brackets: $[a]$, $[x]$, $[y]$, etc. Then $a \in [a]$ simply means that $a$, as an exact real or complex number, is contained in the mathematical set represented by $[a]$. These notations extend entrywise to vectors: $a \in [a]$ with $a \in \mathbb{R}^n$ (or $\mathbb{C}^n$) and $[a] \in \mathbb{R}^n$ (or $\mathbb{C}^n$) if $a_i \in [a_i]$ for all $i$. Analogous notations are used for polynomials (given by their coefficients) and matrices (given by their entries). Whenever $a$ is already defined but not $[a]$, the latter denotes the singleton $\{a\}$ by default. Moreover, given $a \in \mathbb{R}$ or $\mathbb{C}$ and $r \in \mathbb{R}_{\ge 0}$, $[a \pm r]$ is the interval of radius $r$ centred at $a$. $[a \pm r] \in \mathbb{R}^n$ or $\mathbb{C}^n$ is defined analogously for $a \in \mathbb{R}^n$ or $\mathbb{C}^n$ and $r \in \mathbb{R}_{\ge 0}$. The magnitude of $[a]$ is defined as $\max\{|a|, a \in [a]\} \in \mathbb{R}_{\ge 0}$, and it is extended entrywise to vectors and matrices: $\text{mag}[a] \in \mathbb{R}^n$, $\text{mag}[A] \in \mathbb{R}^{n \times k}_{\ge 0}$.

A function $[f] : \mathbb{C}^k \to \mathbb{C}$ is said to be an interval extension of $f : \mathbb{C}^k \to \mathbb{C}$ if for all $a_1 \in [a_1], \dots, a_k \in [a_k]$, $f(a_1, \dots, a_k) \in [f]([a_1], \dots, [a_k])$. Most interval arithmetic libraries provide interval extensions for basic functions such as arithmetic operations $+$, $-$, $\times$, $\div$ and some elementary functions. They will be the building blocks of our validation algorithms.

Validation using Newton-like fixed-point operators. Let $H: \mathbb{C}^k \to \mathbb{C}^k$ be a system of equations, and let $x^* \in \mathbb{C}^k$ denote an exact solution of $H(x) = 0$. In many situations, $x^*$ is not expressed as a finite composition of arithmetic operations and elementary functions, so that basic interval arithmetic is not sufficient to get an interval enclosure of it.

<figure><img src="image_2.png" /><figcaption>Figure 1: Theorem 3.2 in dimension $k = 2$, with the image (green solid shape) of the box $[\tilde{x} \pm r]$ contained in the dashed box $[N_H(\tilde{x}) \pm \Lambda r]$, thus satisfying inequality (6).</figcaption></figure>

A posteriori validation consists in expressing $x^*$ as a fixed point of a contracting operator, then deducing rigorous and tight error bounds for any input approximation $\tilde{x}$ of $x^*$. A popular option for locally invertible $H$ is to use a Newton-like operator:

$$N_H(x) := x - A H(x), \quad (5)$$

where $A \in \mathbb{C}^{k \times k}$ is an approximate inverse of the Jacobian $J_H(\tilde{x})$ of $H$ at $\tilde{x}$. Then the zeros of $H$ are exactly the fixed point of $N_H$ (provided $A$ is injective, which will be a byproduct of the validation). We will rigorously check that $N_H$ is a contraction around $\tilde{x}$.

Definition 3.1. $N_H$ is said to be $\Lambda(r)$-Lipschitz over $[\tilde{x} \pm r] \in \mathbb{C}^k$ with $r \in \mathbb{R}_{\ge 0}^k$ if $\Lambda(r) \in \mathbb{R}_{\ge 0}^{k \times k}$ satisfies:

$$|N_H(x) - N_H(x')| \le \Lambda(r) |x - x'| \quad \text{for all } x, x' \in [\tilde{x} \pm r].$$

Such a Lipschitz matrix $\Lambda(r)$ for $N_H$ can be rigorously computed by bounding the Jacobian of $N_H$ over $x \in [\tilde{x} \pm r]$:

$$J_{N_H}(x) = I_k - A J_H(x).$$

We also need to compute a rigorous bound $\delta \in \mathbb{R}_{\ge 0}^k$ on the defect:

$$\delta \ge |N_H(\tilde{x}) - \tilde{x}| = |AH(\tilde{x})|.$$

Then the following generalisation of Banach's fixed-point principle (see e.g., [18, Thm. 1], or [21, Thm. 10.6] for an analogous statement due to Krawczyk) asserts the existence and uniqueness of a fixed point $x^*$ in $[\tilde{x} \pm r]$.

Theorem 3.2. If the following inequality holds entrywise:

$$\Lambda(r) r + \delta < r, \quad (6)$$

then $N_H$ is a contraction over $[\tilde{x} \pm r]:$ the spectral radius (i.e., the modulus of the largest eigenvalue) $\lambda := \rho(\Lambda(r))$ satisfies $\lambda < 1$, and $N_H$ admits a unique fixed point $x^*$ inside $[\tilde{x} \pm r]$.

For the sake of completeness, we provide a routine `VALIDATEDBOUNDS` which automates the process of finding a tight $r$ satisfying (6) above. It is similar to Rohn's routine [18, §3] or Rump's $\epsilon$-inflation method [21, Algo. 10.7], except that we here assume $\Lambda(r)$ to depend on $r$, since we are concerned with nonlinear systems. Moreover $\Lambda$ is assumed to be nondecreasing: if $r \le r'$ (entrywise), then $\Lambda(r) \le \Lambda(r')$ (entrywise).
**Algorithm 4 VALIDATEDBOUNDS(Λ, δ)**

**Input:** procedure $\mathbf{r} \in \mathbb{R}_{\ge 0}^k \mapsto \Lambda(\mathbf{r}) \in \mathbb{R}_{\ge 0}^{k \times k}$, and $\delta \in \mathbb{R}_{\ge 0}^k$

**Output:** $\mathbf{r} \in \mathbb{R}_{\ge 0}^k$ satisfying $\Lambda(\mathbf{r})\mathbf{r} + \delta < \mathbf{r}$, or failure message

**Parameters:** $\varepsilon = 10^{-10}$, $t_{\max} = 20$ (default values)

1. $\eta_i \leftarrow \begin{cases} \varepsilon\delta_i & \text{if } \delta_i > 0 \\ \varepsilon\delta_0 & \text{if } \delta_i = 0 \end{cases}$ for $i = 1 \dots k$, where $\delta_0 \leftarrow \min_{\delta_i > 0} \delta_i$

2. $\delta^+ \leftarrow \delta + \eta$

3. $\mathbf{r}_0 \leftarrow 0 \in \mathbb{R}^k$

4. **for** $t = 0 \dots t_{\max} - 1$ **do**

5. &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{r}_{t+1} \leftarrow \Lambda(\mathbf{r}_t) \mathbf{r}_t + \delta^+$

6. &nbsp;&nbsp;&nbsp;&nbsp; **if** $\mathbf{r}_{t+1} - \mathbf{r}_t < \eta$ **then** return $\mathbf{r}_t$ **end if**

7. **end for**

8. **return** FAIL

**PROPOSITION 3.3.** If $\Lambda : \mathbb{R}_{\ge 0}^k \to \mathbb{R}_{\ge 0}^{k \times k}$ is continuous nondecreasing and $\delta \ne 0 \in \mathbb{R}_{\ge 0}^k$, then VALIDATEDBOUNDS($\Lambda, \delta$) satisfies:

(i) If it succeeds in a finite number of iterations, then the computed vector $\mathbf{r}$ satisfies inequality (6).

(ii) If inequality (6) with $\delta$ replaced by the $\varepsilon$-inflated vector $\delta^+$ (see lines 1–2) admits a solution, then the algorithm computes a solution $\mathbf{r}$ satisfying (6) in a finite number of iterations.

**PROOF.** The proof is similar to the linear case in [18, Thm. 2]. □

*The overdetermined case.* Particular care is needed when extending this strategy to the overdetermined case $H : \mathbb{C}^k \to \mathbb{C}^n$ with $k < n$. Linear overdetermined systems have been investigated in [18]. Algorithms VALGCDROOT and VALPEJROOT in the following rely on overdetermined nonlinear systems (7) and (3).

First, the Jacobian $J_H(\tilde{\mathbf{x}}) \in \mathbb{C}^{n \times k}$ is still assumed to have full column rank, but it is no longer square and hence not invertible. Therefore, to define $N_H(5)$, one takes for $A$ a numerical approximation of its pseudoinverse:

$$A = J^+ = (J^* J)^{-1} J^* \in \mathbb{C}^{k \times n}, \quad \text{where } J = J_H(\tilde{\mathbf{x}}) \text{ and } J^* = \bar{J}^T.$$

Since $J^+J = I_k$, $N_H$ is still contracting locally around $\tilde{\mathbf{x}}$, and VALIDATEDBOUNDS provides rigorous bounds $\mathbf{r}$. However, $A$ is no longer injective, implying that the fixed point $\mathbf{x}^*$ of $N_H$ in $[\tilde{\mathbf{x}} \pm \mathbf{r}]$ (by Theorem 3.2) may no longer be a solution of $H(\mathbf{x}) = 0$. For instance, overdetermined systems have no solution in the generic case. This explains the careful phrasing of Propositions 3.5 and 3.6 below, and the need for a priori root enclosures for VALPEJROOT.

## 3.2 Validated initial enclosures

Here we assume that $p \in \mathbb{C}_n[x]$ of multiplicity structure $\ell \in \mathbb{N}_{\ge 0}^k$ is represented by an interval polynomial $[p] \in \mathbb{C}_n[x]$, and that approximations $\tilde{z} \in \mathbb{C}^k$ of its roots have already been computed. Algorithm VALGCDROOT([p], $\ell$, $\tilde{z}$) constructs enclosures $[z] \in \mathbb{C}^k$ of its exact roots, based on an alternative system to Eq. (3). By Lemma 2.5, the roots of $p$ are the simple roots of the square-free cofactor $v$:

$$wp - vp' = 0, \quad \text{and} \quad \prod_{i=1}^{k} (x - z_i) = v.$$

Let $v \in \mathbb{C}^k$ and $w \in \mathbb{C}^{k-1}$ denote the coefficients of $v$ and $w$ except the (fixed) leading ones: $v = x^k + v_{k-1}x^{k-1} + \dots + v_0$, and $w = nx^{k-1} + w_{k-2}x^{k-2} + \dots + w_0$. Form the vector $\mathbf{x} = (\mathbf{w}, -\mathbf{v}, \mathbf{z}) \in \mathbb{C}^{3k-1}$, and let $\mathbf{b} \in \mathbb{C}^{n+k-1}$ denote the coefficients of $x^k p' - nx^{k-1} p$. Then, with 1 denoting $(1, \dots, 1) \in \mathbb{N}_{\ge 0}^k$, the system to solve is:

$$F_k(\mathbf{x}) = F_k\left(\begin{pmatrix} \mathbf{w} \\ -\mathbf{v} \\ \mathbf{z} \end{pmatrix}\right) = \begin{pmatrix} \mathbf{b} \\ 0 \end{pmatrix} \iff \left\{ \begin{array}{l} \Delta_{k-1}(p) \begin{pmatrix} \mathbf{w} \\ -\mathbf{v} \end{pmatrix} = \mathbf{b}, \\ G_1(\mathbf{z}) - \mathbf{v} = 0. \end{array} \right. \quad (7)$$

This system has $n+2k-1$ equations in $3k-1$ variables. It is hence overdetermined too. However, the overdetermination entirely lies in the upper part of the system, which is linear. This is the key point making VALGCDROOT not relying on a priori root enclosures.

The Jacobian of $F_k$ is the $(n+2k-1) \times (3k-1)$ matrix:

$$J_{F_k} \begin{pmatrix} \mathbf{w} \\ -\mathbf{v} \\ \mathbf{z} \end{pmatrix} = \begin{pmatrix} \Delta_{k-1}(p) & O_{(n+k-1) \times k} \\ (O_{k \times (k-1)} & I_k & J_{G_1}(\mathbf{z}) \end{pmatrix}, \quad (8)$$

which is necessarily non-surjective when $k < n$. One builds the Newton-like validation operator:

$$N_{F_k}(\mathbf{x}) = \mathbf{x} - A(F_k(\mathbf{x}) - (\mathbf{b}, 0)), \quad (9)$$

where $A$ is computed numerically as the pseudoinverse of $J_{F_k}(\mathbf{x})$ according to the following Lemma.

**LEMMA 3.4.** The Jacobian $J_{F_k}(\mathbf{x})$ is injective whenever $p \in \mathbb{C}_n[x]$ has $k$ distinct roots $z_i$. Therefore, the pseudoinverse $J_{F_k}(\mathbf{x})^+$ satisfies $J_{F_k}(\mathbf{x})^+ J_{F_k}(\mathbf{x}) = I_{3k-1}$, and it is equal to:

$$J_{F_k}(\mathbf{x})^+ = \begin{pmatrix} A_1 & O_{(2k-1) \times k} \\ A_3 & A_2 \end{pmatrix}, \\ \text{where } A_1 = \Delta_{k-1}(p)^+, \quad A_2 = J_{G_1}(\mathbf{z})^{-1},$$

$$A_3 = -A_2 (O_{k \times (k-1)} \quad I_k) A_1 = -A_2 (A_{1,ij})_{\substack{k \le i \le 2k-1 \\ 1 \le j \le n+k-1}}.$$

**PROOF.** The injectivity of $J_{F_k}(\mathbf{x})$ follows from the one of $\Delta_{k-1}(p)$ by Lemma 2.5, together with the invertibility of $J_{G_1}(\mathbf{z})$. It automatically implies that $J_{F_k}(\mathbf{x})^+$ is a left inverse for $J_{F_k}(\mathbf{x})$. Finally, for any $c \in \mathbb{C}^{n+2k-1}$, $\mathbf{y} = J_{F_k}(\mathbf{x})^+ \mathbf{c}$ is by definition the least-squares solution of $J_{F_k}(\mathbf{x})\mathbf{y} = \mathbf{c}$. Since $J_{F_k}(\mathbf{x})$ is lower triangular by block and $J_{G_1}(\mathbf{z})$ is invertible, $\mathbf{y}$ is built by using $A_1 := S_{k-1}(p)^+ \text{ to define its } (\mathbf{w}, -\mathbf{v})$ component (thus minimising the errors on the first } n+k-1 \text{ equations), and by adjusting the } z \text{ component of } \mathbf{y} \text{ so as to cancel the remaining } k \text{ equations. This leads to the given block matrix expression for } J_{F_k}(\mathbf{x})^+.$$ □

Once $N_{F_k}$ is built, Algorithm VALGCDROOT proceeds by rigorously satisfying the hypotheses of Theorem 3.2 using interval arithmetic to deduce bounds $\mathbf{r} \in \mathbb{R}_{\ge 0}^{3k-1}$ for an approximation $\tilde{\mathbf{x}} = (\tilde{w}, -\tilde{v}, \mathbf{z})$ of system (7).

**PROPOSITION 3.5.** Let $[p] \in \mathbb{C}_n[x]$, $\ell \in \mathbb{N}_{\ge 0}^k$ and $\tilde{z} \in \mathbb{C}^k$. If VALGCDROOT([p], $\ell$, $\tilde{z}$) does not fail, then it computes enclosures $[z] \in \mathbb{C}^k$ such that for all $p \in [p]$ with multiplicity structure $\ell$, there exists a unique vector of roots $z \in [z]$ such that $p = p_n \prod_{i=1}^k (x - z_i)^{\ell_i}$.

Moreover, no $p \in [p]$ can have less than $k$ distinct roots, and those having exactly $k$ distinct ones have multiplicity structure $\ell$.
**Algorithm 5 VALGCDROOT([p], ℓ, ˜z)**

**Input:** interval polynomial $[p] \in \mathbb{C}_n[x]$, multiplicity structure $\ell \in \mathbb{N}_{>0}^k$, root approximations $\tilde{z} \in \mathbb{C}^k$

**Output:** root enclosures $[z] \in \mathbb{C}^k$ or failure message

1: compute $\tilde{v}, \tilde{w}$ as in (2) to form $\tilde{x} = (\tilde{w}, -\tilde{v}, \tilde{z})$

2: compute $A \approx J_{F_k}(\tilde{x})^+$ numerically, as in Lemma 3.4
  ▷ From now on, use interval arithmetic

3: $\delta \leftarrow \text{mag}([A]([F_k](\tilde{x}) - ([b], 0)))$

4: $\Lambda : \mathbf{r} \in \mathbb{R}^{3k-1}_{>0} \mapsto \text{mag}(I_{3k-1} - [A][J_{F_k}]([\tilde{x} \pm \mathbf{r}]))$

5: $\mathbf{r} \leftarrow \text{VALIDATEDBOUNDS}(\Lambda, \delta)$

6: $[\mathbf{x}] = ([\mathbf{w}], -[\mathbf{v}], [z]) \leftarrow [\tilde{x} \pm \mathbf{r}]$

7: **for** $i = 1 \dots k$ **do**

8: &nbsp;&nbsp;&nbsp;&nbsp;**if** $\frac{[w]}{[v]'}([z_i]) \cap \mathbb{N} \neq \{\ell_i\}$ **then** **return** **FAIL** **end if**

9: **end for**

10: **return** $[z]$

**PROOF.** Given the candidate solution $\tilde{x} = (\tilde{w}, -\tilde{v}, \tilde{z})$ computed in line 1 and the Newton-like operator $N_{F_k}$ (9) defined from the pseudoinverse $A$ computed in line 2, the use of interval arithmetic beginning on line 3 guarantees that for any $p \in [p]$:

(1) $\delta$ is a rigorous entrywise upper bound for the defect $N_{F_k}(\tilde{x}) - \tilde{x}$;
(2) $\Lambda(\mathbf{r})$ with arbitrary $\mathbf{r} \in \mathbb{R}_{\ge 0}$ is a rigorous entrywise bound for the Jacobian of $N_{F_k}$ taken at any point inside the interval vector $[\tilde{x} \pm \mathbf{r}]$, hence a valid Lipschitz matrix for $N_{F_k}$ over this domain.

According to Theorem 3.2 and Proposition 3.3, the bounds $\mathbf{r}$ computed in line 5 (if VALIDATEDBOUNDS does not fail), ensure that the spectral radius of $\Lambda(\mathbf{r})$ is smaller than 1, that the Newton-like operator $N_{F_k}$ (for any $p \in [p]$) is a contraction over the $\mathbf{r}$-inflated interval vector $[\mathbf{x}] \in \mathbb{C}^{3k-1}$, and therefore admits a unique fixed point $\mathbf{x}^* = (\mathbf{w}^*, -\mathbf{v}^*, \mathbf{z}^*)$ in it. Moreover, since $\Lambda(\mathbf{r})$ has the form $\begin{pmatrix} \Lambda_1 & O_{(2k-1)\times k} \\ \Lambda_3 & \Lambda_2(\mathbf{r}) \end{pmatrix}$ (with $\Lambda_1, \Lambda_3$ independent of $\mathbf{r}$), $\rho(\Lambda(\mathbf{r})) < 1$ implies that $\rho(\Lambda_1) < 1$ and $\rho(\Lambda_2(\mathbf{r})) < 1$, so that $A_1 \Delta_{k-1}(p)$ and $A_2 J_{G_1}(z)$ are invertible (hence also $A_2$ since it is a square matrix).

Let us now assume that $p \in [p]$ has the desired multiplicity structure $\ell$. Then by Lemma 2.5, there are unique polynomials $v$ and $w$ of the form $v = x^k + v_{k-1}x^{k-1} + \dots + v_0$ and $w = nx^{k-1} + w_{k-2}x^{k-2} + \dots + w_0$ satisfying $wp - vp' = 0$. The corresponding coefficients therefore satisfy $\Delta_{k-1}(p) \begin{pmatrix} w \\ -v \end{pmatrix} = \mathbf{b}$, and since $A_1 \Delta_{k-1}(p)$ is invertible, they are equal to $(\mathbf{w}^*, -\mathbf{v}^*)$. Now from the last $k$ equations of $A (F_k(\mathbf{x}^*) - (\mathbf{b}, 0)) = 0$,

$$A_3 \underbrace{\left( \Delta_{k-1}(p) \begin{pmatrix} \mathbf{w}^* \\ \mathbf{v}^* \end{pmatrix} - \mathbf{b} \right)}_{=0} + A_2 (G_1(\mathbf{z}^*) - \mathbf{v}^*) = 0,$$

we deduce that $G_1(\mathbf{z}^*) = \mathbf{v}^*$, by invertibility of $A_2$. Hence, $\mathbf{z}^*$ are the (simple) roots of $v$, which are also the (multiple) roots of $p$.

Finally, the tests in lines 7–9 carried out in interval arithmetic guarantee that the obtained roots are matched with the correct multiplicities, following Lemma 2.6.

Now for the second claim, if $p \in [p]$ has less than $k$ distinct roots, then $\Delta_{k-1}(p)$ is rank deficient, and so is the Jacobian $J_{F_k}$. Therefore $\rho(\Lambda(\mathbf{r})) \ge 1$ and line 5 will necessarily raise an error. $\square$

## 3.3 Refined validated enclosures

Given a polynomial $p \in \mathbb{C}_n[x]$ of multiplicity structure $\ell \in \mathbb{N}_{>0}^k$, represented by $[p] \in \mathbb{C}_n[x]$, Algorithm VALPEJROOT aims at validating approximate roots $\tilde{z} \in \mathbb{C}^k$, seen as an approximate solution of system (3). The Newton-like validation operator,

$$N_{G_{\ell}}(z) := z - A (G_{\ell}(z) - \mathbf{p}), \quad (10)$$

is built from a numerical approximation $A$ of the pseudoinverse $J_{G_{\ell}}(\tilde{z})^+ \in \mathbb{C}^{k \times n}$.

Using interval arithmetic and Theorem 3.2, VALPEJROOT constructs enclosures $[z] \in \mathbb{C}^k$ around $\tilde{z}$ containing a unique fixed point $\mathbf{z}^*$ of $N_{G_{\ell}}$. However, since $A$ is not injective, one cannot deduce that $\mathbf{z}^*$ is a solution of system (3). Moreover, the argument used in the proof of Proposition 3.5 for VALGCDROOT was specific to the structure of system (7), whose overdetermined part was linear. This is why VALPEJROOT([p], $\ell$, $\tilde{z}$, [z₀]) relies on an extra argument: an a priori, possibly not very tight, enclosure $[z₀] \in \mathbb{C}^k$ of the roots. This may be instantiated as the output of VALGCDROOT([p], $\ell$, $\tilde{z}$).

**Algorithm 6 VALPEJROOT([p], ℓ, ˜z, [z₀])**

**Input:** $[p] \in \mathbb{C}_n[x], \ell \in \mathbb{N}_{>0}^k$, root approximations $\tilde{z} \in \mathbb{C}^k$, and initial root enclosures $[z_0] \in \mathbb{C}^k$

**Output:** refined root enclosures $[z] \in \mathbb{C}^k$, or failure message

1: compute $A \approx J_{G_{\ell}}(\tilde{z})^+$ numerically

▷ From now on, use interval arithmetic

2: $[\mathbf{p}] \leftarrow (\frac{p_0}{p_n}, \dots, \frac{p_{n-1}}{p_n})$

3: $\delta \leftarrow \text{mag}([A]([G_{\ell}](\tilde{z}) - [\mathbf{p}]))$

4: $\Lambda_0 \leftarrow \text{mag}(I_k - [A][J_{G_{\ell}}][z_0])$

5: VALIDATEDBOUNDS($\Lambda_0, \delta$) ▷ raises an error if $\rho(\Lambda_0) \ge 1$

6: $\Lambda : \mathbf{r} \in \mathbb{R}_{>0}^k \mapsto \text{mag}(I_k - [A][J_{G_{\ell}}]([\tilde{z} \pm \mathbf{r}]))$

7: $\mathbf{r} \leftarrow \text{VALIDATEDBOUNDS}(\Lambda, \delta)$

8: $[z] \leftarrow [\tilde{z} \pm \mathbf{r}]$

9: **if** $[z] \not\subseteq [z_0]$ **then** **return** **FAIL** **end if**

10: **return** $[z]$

**PROPOSITION 3.6.** Let $[p] \in \mathbb{C}_n[x], \ell \in \mathbb{N}_{>0}^k, \tilde{z} \in \mathbb{C}^k$ and $[z_0] \in \mathbb{C}^k$. Assume that for all $p \in [p]$ with multiplicity structure $\ell$, there exists a unique vector of roots $z \in [z_0]$ such that $p = p_n \prod_{i=1}^k (x - z_i)^{\ell_i}$. Then if VALPEJROOT([p], $\ell$, $\tilde{z}$, [z₀]) does not fail, it computes refined root enclosures $[z] \subseteq [z_0]$ such that $z \in [z]$.

**PROOF.** The use of interval arithmetic after line 2 ensures that $\delta$ is a rigorous upper bound for the defect $N_{G_{\ell}}(\tilde{z}) - \tilde{z}$, and $\Lambda_0$ (resp. $\Lambda(\mathbf{r})$) a rigorous Lipschitz matrix for $N_{G_{\ell}}$ over $[z_0]$ (resp. $[\tilde{z} \pm \mathbf{r}]$).

The first call to VALIDATEDBOUNDS with the constant Lipschitz matrix $\mathbf{r} \mapsto \Lambda_0$ (the second argument, $\delta$, is not relevant) ensures that $\rho(\Lambda_0) < 1$ by Proposition 3.3 and Theorem 3.2 (otherwise this routine raises an error). As a result, the initial root enclosure $[z_0]$ contains at most one fixed point of $N_{G_{\ell}}$.

The second call to VALIDATEDBOUNDS, if it does not fail, returns an error vector $\mathbf{r} \in \mathbb{R}_{>0}^k$ such that $N_{G_{\ell}}$ contains a unique fixed point $\mathbf{z}^*$ inside $[z] := [\tilde{z} \pm \mathbf{r}]$. Also, if the inclusion $[z] \subseteq [z_0]$ in line 9 is satisfied, then it proves that this fixed point is also the unique fixed point of $N_{G_{\ell}}$ in $[z_0]$.
<table><caption>Table 2: Maximum relative errors on the root approximations and enclosures for $p_m$ - Eq. (4). Results in parentheses for VALPEJROOT indicate not fully rigorous enclosures due to the failure of the contraction test in line 5.</caption><thead><tr><th>m</th><th>GCDROOT*</th><th>PEJROOT</th><th>VALGCDROOT</th><th>VALPEJROOT</th></tr></thead><tbody><tr><td>4</td><td>1.25e-07</td><td>4.59e-15</td><td>5.62e-07</td><td>4.28e-13</td></tr><tr><td>5</td><td>4.99e-07</td><td>1.33e-14</td><td>1.05e-06</td><td>4.44e-13</td></tr><tr><td>6</td><td>4.01e-07</td><td>9.92e-15</td><td>2.08e-06</td><td>4.98e-13</td></tr><tr><td>7</td><td>4.35e-07</td><td>7.55e-15</td><td>2.89e-06</td><td>4.78e-13</td></tr><tr><td>8</td><td>6.47e-07</td><td>4.59e-15</td><td>4.80e-06</td><td>4.99e-13</td></tr><tr><td>9</td><td>1.09e-06</td><td>3.11e-15</td><td>6.74e-06</td><td>4.99e-13</td></tr><tr><td>10</td><td>1.13e-06</td><td>6.77e-15</td><td>8.92e-06</td><td>5.11e-13</td></tr><tr><td>11</td><td>1.83e-06</td><td>8.88e-16</td><td>1.22e-05</td><td>4.93e-13</td></tr><tr><td>12</td><td>9.01e-07</td><td>5.92e-16</td><td>1.61e-05</td><td>5.13e-13</td></tr><tr><td>13</td><td>3.29e-06</td><td>7.40e-15</td><td>2.00e-05</td><td>5.18e-13</td></tr><tr><td>14</td><td>2.52e-06</td><td>3.40e-15</td><td>2.62e-05</td><td>(5.29e-13)</td></tr><tr><td>15</td><td>6.33e-06</td><td>5.03e-15</td><td>3.14e-05</td><td>(5.32e-13)</td></tr><tr><td>16</td><td>1.24e-05</td><td>3.55e-15</td><td>—</td><td>(5.26e-13)</td></tr><tr><td>17</td><td>1.58e-05</td><td>7.85e-15</td><td>—</td><td>(5.34e-13)</td></tr><tr><td>18</td><td>1.36e-05</td><td>9.33e-15</td><td>—</td><td>(5.31e-13)</td></tr><tr><td>19</td><td>2.99e-05</td><td>1.33e-15</td><td>—</td><td>(5.44e-13)</td></tr><tr><td>20</td><td>4.91e-05</td><td>8.88e-16</td><td>—</td><td>(5.37e-13)</td></tr></tbody></table>

Now let $p \in [p]$ have multiplicity structure $\ell$. By the assumption on $[z_0]$, there exists a unique $z \in [z_0]$ such that $p = p_n \prod (x - z_i)^{\ell_i}$. Clearly, $z$ is a fixed point of $N_{G_\ell}$. Hence, $z = z^* \in [z]$. □

**Remark 3.7.** Given $[p] \in \mathbb{C}_n[x]$, Algorithms VALGCDROOT and VALPEJROOT compute enclosures $[z] \in \mathbb{C}^k$ for the roots of any $p \in [p]$ of multiplicity structure $\ell$. They do however not prove that $[p]$ contains at least one such $p$. This can be done by checking that $p_n \prod_{i=1}^k (x - \tilde{z}_i) \in [p]$, with $\tilde{z}$ the computed numerical roots.

## 3.4 Validated Numerical Examples

**Example 1 (validated version).** We computed tight coefficientwise enclosures $[p_m]$ (up to machine double precision) for polynomials $p_m$ of Example 1, and applied the full validation chain: from the root approximations $\tilde{z}$ computed by our routine GCDROOT* followed by Zeng's PEJROOT, an initial enclosure $[z_0]$ is computed by VALGCDROOT and refined by VALPEJROOT. The maximum relative errors obtained on the four roots are summarised in Table 2. We observe that the relative radius of $[z_0]$ grows with $m$ roughly as fast as the relative error of the initial approximations computed by GCDROOT*. This follows from the growth with $m$ of the condition number of Sylvester matrices, and hence of system (7). On the other hand, the growth with $m$ of the condition number associated to system (3) is moderate, which explains why only two to three digits are lost in the intervals computed by VALPEJROOT. However, for $m \ge 14$, VALPEJROOT cannot fully certify its enclosures due to the failure of the contraction test in line 5: VALGCDROOT either fails to compute initial enclosures $[z_0]$, or they are too large. Nevertheless, these results are encouraging and show that the full validation chain solving Problem 2 on this example works for degrees up to 130 and multiplicities up to 52 using double precision only.

<figure><img src="image_3.png" /><figcaption>Figure 2: Accuracy of the computed root enclosures in function of accuracy of interval coefficients of $[p_m]$ - Eq. (4).</figcaption></figure>

**Example 1 (continued: the effect of uncertainties).** For $m$ up to 5, Figure 2 shows the relative accuracy of the root enclosures $[z]$ in function of the relative accuracy of the interval coefficients of $[p_m]$. First, we observe that the number of lost bits in the result is equal to the number of lost bits in the input plus a constant corresponding to the conditioning. This corroborates our claim that this validation method has finite conditioning: it is not subject to the “attainable accuracy barrier”.

Unfortunately, the input accuracy threshold over which our method fails decreases with $m$, although the computed bounds themselves do not increase excessively. The first reason is the condition number of Sylvester matrices used by VALGCDROOT. Another source of failure we observed is the potentially very large overestimations in the interval evaluations of the maps $G_\ell$ and $J_{G_\ell}$ over non-thin intervals $[z]$, something well-known as the wrapping effect in interval analysis [21, §9.2].

**Example 2 (nearby multiple roots).** Figure 3 compares obtained accuracies of the approximations and validated enclosures for

$$ (x - 1 + \varepsilon)^{m_1} (x - 1)^{m_2} (x + 0.5)^{m_3}, \quad (11) $$

already used as example in [28, §5.2]. The accuracy deteriorates when the gap $\varepsilon$ between the two nearby multiple roots tends to 0. This is a consequence of the growth of the structure-preserving condition number when we get closer to the pejorative submanifold $\Pi_{m_1+m_2,m_3}$. On the other hand, the multiplicities do not seem to impact the accuracy significantly. The overestimation factor of the validated enclosures compared to the approximations is rather mild. However, higher multiplicities lower the threshold on $\varepsilon$ over which the combination of VALGCDROOT and VALPEJROOT fails.

**Example 3 (complex roots).** The following polynomial:

$$ (x^2 - x + 1)^{m_1} (x^2 + 4x + 7)^{m_2} (x^2 - x - 1)^{m_3} (x^2 + 2x + 2)^{m_4} \quad (12) $$

has 8 complex roots depicted in Figure 4. The relative accuracies of the computed root approximations and enclosures for different values of the $m_i$ are summarised in Table 3. It turns out that, in
<figure><img src="image_4.png" /><figcaption>Figure 3: Relative accuracy of root approximations (dashed lines) and root enclosures (solid lines) for the polynomial (11) in function of $\epsilon$ and multiplicities ($m_1, m_2, m_3$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Complex roots of polynomial (12).</figcaption></figure>

this example, the condition number of $J_{G\epsilon}(z)$ grows significantly faster than the condition number of the Sylvester matrices $\Delta_k(p)$. As a result, the approximations produced by $GCDROOT^*$ are better than those of $PEJROOT$ (we thus kept the former for the validation), and $VALGCDROOT$ performs better than $VALPEJROOT$. Such cases are not discussed in [28, 30], and further research is needed to predict when they occur.

## 4 CONCLUDING REMARKS

In this article, we improved Zeng’s numerical method [28, 30] to solve Problem 1 and proposed a validation method based on two algorithms $VALGCDROOT$ and $VALPEJROOT$ to solve Problem 2.

A first future work is to refine our implementations of the validation routines so that their performances get closer to the purely numerical routines in terms of high multiplicities and obtained

<table><caption>Table 3: Maximum relative errors on the root approximations and enclosures for the polynomials (12).</caption><thead><tr><th>m<sub>1</sub>, m<sub>2</sub>, m<sub>3</sub>, m<sub>4</sub></th><th>GCDROOT*</th><th>PEJROOT</th><th>VALGCDROOT</th><th>VALPEJROOT</th></tr></thead><tbody><tr><td>2, 2, 1, 1</td><td>2.47e-15</td><td>3.98e-15</td><td>3.17e-13</td><td>6.75e-13</td></tr><tr><td>3, 1, 1, 3</td><td>4.87e-15</td><td>3.95e-14</td><td>6.01e-12</td><td>8.67e-11</td></tr><tr><td>5, 3, 3, 1</td><td>4.39e-15</td><td>6.83e-14</td><td>4.56e-13</td><td>5.12e-10</td></tr><tr><td>6, 3, 3, 1</td><td>2.02e-14</td><td>4.22e-14</td><td>7.97e-13</td><td>2.30e-09</td></tr><tr><td>2, 5, 5, 4</td><td>5.02e-13</td><td>1.16e-12</td><td>5.52e-11</td><td>—</td></tr><tr><td>8, 5, 2, 1</td><td>6.95e-15</td><td>1.59e-12</td><td>6.01e-13</td><td>—</td></tr><tr><td>8, 10, 3, 7</td><td>2.91e-10</td><td>8.12e-10</td><td>6.77e-09</td><td>—</td></tr><tr><td>20, 11, 7, 5</td><td>1.25e-13</td><td>1.26e-04</td><td>3.13e-12</td><td>—</td></tr><tr><td>25, 30, 3, 5</td><td>6.42e-10</td><td>1.83e-02</td><td>3.04e-08</td><td>—</td></tr><tr><td>25, 30, 3, 10</td><td>1.15e-07</td><td>8.04e-03</td><td>—</td><td>—</td></tr></tbody></table>

accuracy without using multi-precision. This, in particular, concerns Algorithm $VALPEJROOT$ where refined evaluation methods could lower the overapproximations on the defect and Lipschitz matrices, e.g., Discrete Fourier Transform (DFT) based evaluation methods to reduce the wrapping effect. Also, in cases where the iterative condition number is large, one can consider combining the “one-root-at-a-time” deflation method of [22] with our two-stage validation approach to make it compliant with Problem 2.

On the theoretical side, a thorough complexity analysis of algorithms to solve Problem 2 will require a deeper analysis of the growth of the various condition numbers involved.

**Acknowledgments.** The authors are grateful to Bruno Salvy for useful discussions on detecting and converging to clusters of roots; to John Verzani for his kind help in using the Polynomials.jl package; and to the anonymous referees for relevant suggestions.

## REFERENCES

[1] Bernhard Beckermann and George Labahn. 1998. When are Two Numerical Polynomials Relatively Prime? *Journal of Symbolic Computation* 26, 6 (1998), 677–689. https://doi.org/10.1006/jsco.1998.0234

[2] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. 2017. Julia: A fresh approach to numerical computing. *SIAM review* 59, 1 (2017), 65–98. https://doi.org/10.1137/141000671

[3] Dario A. Bini and Paola Boito. 2007. Structured Matrix-Based Methods for Polynomial $\epsilon$-Gcd: Analysis and Comparisons. In *Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation (ISSAC '07)*. Association for Computing Machinery, New York, NY, USA, 9–16. https://doi.org/10.1145/1277548.1277551

[4] Paola Boito. 2012. Structured matrix based methods for approximate polynomial GCD. Vol. 15. Springer Science & Business Media. https://doi.org/10.1007/978-88-7642-381-9

[5] Robert M. Corless, Patrizia M. Gianni, Barry M. Trager, and Stephen M. Watt. 1995. The Singular Value Decomposition for Polynomial Systems. In *Proceedings of the 1995 International Symposium on Symbolic and Algebraic Computation (ISSAC '95)*. ACM, 195–207. https://doi.org/10.1145/220346.220371

[6] Ioannis Z. Emiris, André Galligo, and Henri Lombardi. 1997. Certified approximate univariate GCDs. *Journal of Pure and Applied Algebra* 117-118 (1997), 229–251. https://doi.org/10.1016/S0022-4049(97)00013-3

[7] Marc Giusti, Grégoire Lecerf, Bruno Salvy, and Jean-Claude Yakoubsohn. 2005. On location and approximation of clusters of zeros of analytic functions. *Foundations of Computational Mathematics* 5 (2005), 257–311. https://doi.org/10.1007/s10208-004-0144-z

[8] Gene H. Golub and Charles F. Van Loan. 1996. *Matrix computations* (3rd ed.). Johns Hopkins University Press. https://doi.org/10.56021/9781421407944

[9] Nicholas J. Higham. 2002. *Accuracy and Stability of Numerical Algorithms* (second ed.). Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898718027

[10] Fredrik Johansson. 2017. Arb: Efficient Arbitrary-Precision Midpoint-Radius Interval Arithmetic. *IEEE Trans. Comput.* 66, 8 (2017), 1281–1292. https://doi.org/10.1109/TT.2017.66.8.8
org/10.1109/TC.2017.2690633

[11] William Kahan. 1972. *Conserving confluence curbs ill-condition*. Technical Report. https://apps.dtic.mil/sti/citations/AD0766916

[12] Angelos Mantzaflaris and Bernard Mourrain. 2011. Deflation and Certified Isolation of Singular Zeros of Polynomial Systems. In *Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation* (San Jose, California, USA) (ISSAC '11). ACM, 249–256. https://doi.org/10.1145/1993886.1993925

[13] Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud. 2009. *Introduction to Interval Analysis*. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898717716

[14] C. Andrew Neff and John H. Reif. 1996. An Efficient Algorithm for the Complex Roots Problem. *Journal of Complexity* 12, 2 (1996), 81–115. https://doi.org/10.1006/jcom.1996.0008

[15] Arnold Neumaier. 1991. *Interval Methods for Systems of Equations*. Cambridge University Press. https://doi.org/10.1017/CBO9780511526473

[16] Victor Y. Pan. 1997. Solving a Polynomial Equation: Some History and Recent Progress. *SIAM Rev.* 39, 2 (1997), 187–220. https://doi.org/10.1137/S0036144595288554

[17] Victor Y. Pan. 2002. Univariate Polynomials: Nearly Optimal Algorithms for Numerical Factorization and Root-finding. *Journal of Symbolic Computation* 33, 5 (2002), 701–733. https://doi.org/10.1006/jsco.2002.0531

[18] Jiří Rohn. 1996. Enclosing solutions of overdetermined systems of linear interval equations. *Reliable Computing* 2 (1996), 167–171. https://doi.org/10.1007/BF02425920

[19] Siegfried M. Rump. 1999. *INTLAB—INTerval LABoratory*. Springer. https://doi.org/10.1007/978-94-017-1247-7_7

[20] Siegfried M. Rump. 2003. Ten methods to bound multiple roots of polynomials. *J. Comput. Appl. Math.* 156, 2 (2003), 403–432. https://doi.org/10.1016/S0377-0427(03)00381-9

[21] Siegfried M. Rump. 2010. Verification methods: Rigorous results using floating-point arithmetic. *Acta Numerica* 19 (2010), 287–449. https://doi.org/10.1017/S096249291000005X

[22] Siegfried M. Rump and Stef Graillat. 2010. Verified error bounds for multiple roots of systems of nonlinear equations. *Numerical Algorithms* 54 (2010), 359–377. https://doi.org/10.1007/s11075-009-9339-3

[23] Tateaki Sasaki and Matu-Tarow Noda. 1989. Approximate square-free decomposition and root-finding of ill-conditioned algebraic equations. *J. Inform. Process* 12, 2 (1989), 159–168.

[24] Arnold Schönhage. 1985. Quasi-GCD computations. *Journal of Complexity* 1, 1 (1985), 118–137. https://doi.org/10.1016/0885-064X(85)90024-X

[25] Hans J Stetter. 2004. *Numerical polynomial algebra*. SIAM.

[26] Lloyd N Trefethen and David Bau. 1997. *Numerical Linear Algebra*. Vol. 181. Society for Industrial and Applied Mathematics.

[27] Warwick Tucker. 2011. Validated Numerics: A Short Introduction to Rigorous Computations. Princeton University Press. https://doi.org/10.2307/j.ctvcm4g18

[28] Zhonggang Zeng. 2003. A Method Computing Multiple Roots of Inexact Polynomials. In *Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation (ISSAC '03)*. ACM, 266–272. https://doi.org/10.1145/860854.860907

[29] Zhonggang Zeng. 2004. Algorithm 835: MultRoot—a Matlab Package for Computing Polynomial Roots and Multiplicities. *ACM Trans. Math. Softw.* 30, 2 (jun 2004), 218–236. https://doi.org/10.1145/992200.992209

[30] Zhonggang Zeng. 2005. Computing multiple roots of inexact polynomials. *Math. Comp.* 74, 250 (2005), 869–903. https://doi.org/10.1090/S0025-5718-04-01692-8

[31] Zhonggang Zeng and Barry H. Dayton. 2004. The Approximate GCD of Inexact Polynomials. In *Proceedings of the 2004 International Symposium on Symbolic and Algebraic Computation (ISSAC '04)*. ACM, 320–327. https://doi.org/10.1145/1005285.1005331