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

HAL
open science

# Validated Numerical Newton-Puiseux Algorithm

Florent Bréhard, Fabien Corinaldesi, Adrien Poteaux

► To cite this version:

Florent Bréhard, Fabien Corinaldesi, Adrien Poteaux. Validated Numerical Newton-Puiseux Algorithm. 2026. hal-05510315

HAL Id: hal-05510315

https://hal.science/hal-05510315v1

Preprint submitted on 13 Feb 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Validated Numerical Newton-Puiseux Algorithm

Florent Bréhard

florent.brehard@univ-lille.fr

Univ. Lille, CNRS, Centrale Lille, UMR

9189 CRISTAL

Lille, France

Fabien Corinaldesi

fabien.corinaldesi@univ-lille.fr

Univ. Lille, CNRS, Centrale Lille, UMR

9189 CRISTAL

Lille, France

Adrien Poteaux

adrien.poteaux@univ-lille.fr

Univ. Lille, CNRS, Centrale Lille, UMR

9189 CRISTAL

Lille, France

## Abstract

We present the first fully validated numerical adaptation of the Newton-Puiseux algorithm. Adapting the modular-numeric strategy of [19] to the ARNP algorithm [22], we develop filters to correctly match the recursive calls of the numeric algorithm against the branches of the trees computed modulo *p*, ensuring the preservation of Newton polygons and multiplicities at every step. Our algorithm is proven to be correct (the output interval Puiseux series contain the exact ones) and convergent (output intervals can be made as tight as desired), together with an arithmetic complexity analysis. We finally illustrate these theoretical results with numerical examples using our prototype implementation in Julia.

## Keywords

Puiseux series, validated numerics, symbolic-numeric algorithm

## 1 Introduction and motivation

This paper provides a validated modular-numeric algorithm that computes Puiseux series of a bivariate polynomial with validated structure and error bounds for its coefficients.

Let $\mathbb{K}$ be a number field, $x$ and $y$ two indeterminates over $\mathbb{K}$ and $F \in \mathbb{K}[x, y]$ a bivariate polynomial primitive, separable and monic<fnref n="1" /> in $y$. We denote $D$ the total degree of $F$, $d = \deg_y(F)$ and $\tau$ the size of its coefficients. Let $\overline{\mathbb{K}}$ be the algebraic closure of $\mathbb{K}$ and $\delta$ the $x$-valuation of the discriminant of $F$ in $y$. The Puiseux theorem states that for any $x_0 \in \overline{\mathbb{K}}$, the roots of $F$ (viewed as a univariate polynomial in $y$) may be expressed as fractional Laurent power series in $(x - x_0)$ with coefficients in $\overline{\mathbb{K}}$. These are the (classical) Puiseux series<fnref n="2" /> of $F$ above $x_0$, fundamental objects of the theory of algebraic curves [4, 30]. Many applications are given in [21, 22].

The main difficulty is usually to compute the singular parts of Puiseux series (that contain the relevant information about the singularities of the associated curve; remaining terms can be computed up to an arbitrary precision in quasi-linear time via Newton iterations). The best arithmetic complexity is $O(D^3)$ [23]. However, the binary complexity remains an important challenge, explaining how Puiseux series mainly remain a theoretical tool not so used in practice. To our knowledge, the only binary complexity results are due to Chistov [6] ("the complexity is polynomial") and Walsh [31, 32] (complexity bound in $O(D^{36} \tau^2)$). And practical behaviour of the symbolic execution of this algorithm is indeed not promising. On another hand, floating point computation of Puiseux series is not straightforward. The slightest approximation causes Newton-Puiseux algorithm to miss essential information and numerical

instability. To overcome these difficulties, a modular-numeric approach has been introduced in [18, 19] (see Section 3.2). Up to now, the numerical computations were made via floating-point computations, without any certification of the result. This papers aims so validate these computations (the modular part of the algorithm has been fully accomplished in [20, 21]; see Section 3.2. We will only add a minor adaptation to fit the ARNP algorithm; see Section 3.2). Combining these different results, we get:

**THEOREM 1.** Let $F \in \mathbb{Z}[x, y]$ of total degree $D$ and coefficient size $\tau$. Given a numerical precision<fnref n="3" /> $\eta$, there exists a deterministic algorithm that either fails from insufficient precision or computes Puiseux all series of $F$ above $x = 0$ with certified structured and validated coefficients in less that $O(D^4 (\tau + \eta))$ binary operations.

**THEOREM 2.** Let $F \in \mathbb{Z}[x, y]$ of total degree $D$ and coefficient size $\tau$. Given a numerical precision<fnref n="3" /> $\eta$, there exists a Las Vegas algorithm that either fails from insufficient precision or computes Puiseux all series above all critical points $F$ with certified structured and validated coefficients in less that $O(D^5 (\tau + \eta + D))$ binary operations.

The paper is organised as follows. We start with a short introduction on validated numerics in Section 2, then briefly remind the Newton-Puiseux algorithm and the modular-numeric strategy in Section 3, before the presentation of our two stage filters in Section 4 (a first filter according to Newton polygons, then a second one according to the multiplicity structures of the characteristic polynomials). Finally, we deal with numerical validation, focusing on the Weierstrass preparation theorem in Section 5, then the Newton-Puiseux algorithm itself in Section 6. We conclude by a few words on our prototype Julia implementation and numerical examples.

Let us conclude this introduction by summarising our contributions and precise what was already done before this paper. The modular-numeric strategy comes from [18, 19], and is kept unchanged. The only novelty here is the required adaptations on polygon trees to fit the ARNP algorithm of [22], who was published after the papers on modular reduction [20, 21]. The main idea of a two-stage filters was given (in french only) in [19]. The Newton polygon filter is here adapted to the new definition of polygon trees, while the filter based on multiplicity structures has been greatly improved thanks to our previous work on validated numerical roots of univariate polynomials [3] (results that are slightly adapted to our context too). The validated Weierstrass preparation theorem mainly comes from our work on Hensel lifting [2]. Finally, the validation of the main algorithm is entirely new. Note that we do not address the divide and conquer algorithm of [23] here. A better (for practical performance) strategy would actually be to factorise $F$ via the OM-algorithm [24], and then to run the ARNP Newton-Puiseux algorithm of [22] on each irreducible factor. We plan to

<footnote n="1">¹the monicity assumption can be easily lifted via Hensel lifting and the use of reciprocal polynomials for roots above infinity. See [23, §4.5] for details.</footnote>

<footnote n="2">²terms written in italics in this introduction are defined in Section 3.1.</footnote>

<footnote n="3">³Analysing the required precision to achieve a target accuracy is left to future work.</footnote>
work on a validated modular numeric version of the OM algorithm
in the future, but the validation part of the current paper is needed
anyway for dealing with the irreducible case.

## 2 A concise introduction to validated numerics

We assume a simplified floating-point (FP) arithmetic model with unbounded exponents and a precision *t* adaptable on demand, i.e., a set **F** of FP numbers *x* = $\frac{m}{2^{t-1}} \cdot 2^e$ with *m*, *e* ∈ **Z** s.t. *m* = 0 or $2^{t-1} \le |m| < 2^t$. Following the IEEE 754 standard [16], the arithmetic operations (+, -, ×, /, √) are correctly rounded w.r.t. the rounding modes *round-up* RU, *round-down* RD and *round-to-nearest* RN, and each rounding induces a relative error $|\tilde{x} - x|/|x| \le u_{\text{prec}}$ for $x \ne 0$ and $\tilde{x} = \text{RN}(x)$, or $2u_{\text{prec}}$ for $\tilde{x} = \text{RU}(x)$ or $\text{RD}(x)$, with $u_{\text{prec}} = 2^{-t}$ the rounding unit.

Ball arithmetic. Ball arithmetic [12, 28] is a computationally effi-
cient variant of Moore's interval arithmetic [15]. Intervals used
to rigorously approximate real numbers are represented using
midpoint-radius balls $\boldsymbol{x} = [m \pm r] = [m - r, m + r]$, with mid($\boldsymbol{x}$) =
$m \in \mathbf{F}_f$ for a high, adaptive precision $t$, and rad($\boldsymbol{x}$) = $r \in \mathbf{F}_{t'}$
with a lower, usually fixed precision $t'$ ($u'_{\text{prec}} = 2^{-t'}$). Hence, most
operations on balls, starting with arithmetic ones, involve only
one high-precision operation in $\mathbf{F}_t$ plus additional low-precision
operations in $\mathbf{F}_{t'}$.

The cornerstone of interval/ball arithmetic, the inclusion prin-
ciple, states that, e.g., $x + y \in x + y$ must hold for any $x \in x$ and
$y \in y$, where boldface letters indicate balls. Arithmetic operations
are notably implemented with this property. We denote by $\square\mathbb{R}$ the
set of real balls, and by $\circ\mathbb{R}$ the set of $0$-balls, i.e., $[\pm r] = [0 \pm r]$ with
$r \in \mathbf{F}_{t'}$. Call $|-|_{\circ\mathbb{R}}$ the operation that maps any $\boldsymbol{x} = [m \pm r] \in \square\mathbb{R}$ to
$|\boldsymbol{x}|_{\circ} = [\pm|m| + r] \in \circ\mathbb{R}$. For any finite-dimensional $\mathbb{R}$-vector space
$V$ with a canonical basis (e.g., $\mathbb{C}$, $\mathbb{C}^n$, $\mathbb{C}[x, y]$, etc.), $\square V$ and $\circ V$ is
the set of balls $v = (v_1, \dots, v_n)$ with real balls $v_i$ for each coordi-
nate, and $v \in v$ means component-wise membership. Also, $\|v\|_{\infty}$
denotes $\max_i |v_i|$. A ball extension of $f: X \approx \mathbb{K}^n \to Y \approx \mathbb{K}^m$ is a
map $\square f: \square X \to \square Y$ s.t. for all $\boldsymbol{x} \in \square X$ and $\boldsymbol{x} \in \boldsymbol{x}$, $f(\boldsymbol{x}) \in \square f(\boldsymbol{x})$.
Moreover, $\square f$ is isotonic if for all $\boldsymbol{x} \subseteq \boldsymbol{x}' \in \square X$, $\square f(\boldsymbol{x}) \subseteq \square f(\boldsymbol{x}')$.
Also, $\square f$ is Lipschitz-continuous at $\boldsymbol{x} \in X$ if there exists $c \in \mathbb{R}^{n \ge 0}$
and $d, d' \ge 0$ such that for all $\boldsymbol{x} \ni \boldsymbol{x}$ and sufficiently small rad($\boldsymbol{x}$),
$u_{\text{prec}}$, $u'_{\text{prec}}$, rad($\square f_i(\boldsymbol{x})$) $\le (c^T \cdot \text{rad}(\boldsymbol{x}) + d) u_{\text{prec}}$ $(1+d' u'_{\text{prec}})$. The
notion extends to a compact $K \subseteq X$ if the constants $c, d, d'$ are
independent of $x$. If $f: X \to Y$ is a function given by a straight-
line program (SLP) (i.e., a directed acyclic graph whose nodes are
arithmetic operations and leaves are input variables or constants),
then its natural ball extension $\square f$ obtained by applying the corre-
sponding ball operations, is isotonic and Lipschitz-continuous over
any compact of its domain of definition (see, e.g., [11, Prop. 3.1]).

For a continuously differentiable function $f: X \to Y$, we also
define a 0-ball extension of its derivative $df$ at some $x$ to be a map
$\circ df(x) : \circ X \to \circ Y$ s.t. for all $e \in \circ X$ and $e, e' \in e$, $df(x+e) \cdot e' \in \circ df(x) \cdot e$. We may simply write $\circ df$ when $x$ is clear from the
context. Such maps $\circ df$ are traditionally implemented as $e \mapsto M_e \cdot e$
with $M_e \in \square \mathbb{R}^{m \times n}$ an enclosure of $df(x+e)$ for any $e \in e$, but we
will also use SLPs to achieve better complexity bounds.

Krawczyk's method. When the above self-validating operations
fail at computing tight enclosures of the solution $x^*$ to an equation

$f(x) = 0$, a posteriori validation using Krawczyk's method [15, §8.2]
is often a good alternative. It reconstructs a rigorous $\boldsymbol{x} = \tilde{\boldsymbol{x}} + \boldsymbol{e} \in \square X$
($e \in \circ X$) from a purely numerical approximation $\tilde{\boldsymbol{x}}$. Let us assume
that $f: X \approx \mathbb{R}^n \to Y \approx \mathbb{R}^m$ is continuously differentiable and that
$df(x)$ is injective in a neighbourhood of $x^*$ (which notably implies
$m \ge n$). The method first computes an effective linear operator
$A: Y \to X$ (e.g., represented by a matrix or a SLP) approximating
the inverse of $df(\tilde{\boldsymbol{x}})$ if $m = n$, or its Moore-Penrose inverse when
$m > n$ [25]. This defines the Krawczyk operator:

$$
\kappa : X \to X, \quad x \mapsto x - A \cdot f(x), \tag{1}
$$

of which $x^*$ is a fixed point. We expect $\kappa$ to be contracting around
$x^*$ since $d\kappa(x) = \text{id}_X - A \cdot df(x)$ should be small by definition of A.

THEOREM 2.1 (BANACH FIXED POINT, see e.g. [11, Thm. 2.1]). If

$$
\kappa(\tilde{x}) + d\kappa(\tilde{x} + e) \cdot e \subseteq \dot{e}, \quad (2)
$$

then $\tilde{x} + e$ contains a unique fixed point of $\kappa$, and $\text{id}_X - d\kappa(\tilde{x} + e) = A \cdot df(x + e)$ is invertible for all $e \in e$.

If $n = m$, both $A$ and $df(x+e)$ are invertible, and the obtained fixed point is a zero of $f$. When $m > n$, we must a priori postulate the existence of a solution to the overdetermined system $f(x) = 0$.

Krawczyk's method, and more precisely Rump's $\epsilon$-inflation strat-
egy [26, §10.5], turns Theorem 2.1 into an effective validation proce-
dure, detailed as Algorithm KRAWCZYK. The following proposition
(with proof postponed to App. A), rigorously establishes its correct-
ness and termination in finite-precision ball arithmetic, without the
restriction that $m = n$ nor that $\circ d\kappa$ must compute with matrices.

**Algorithm 1: KRAWCZYK(d, ◦dκ)**

In: defect $d \in \circ X$, 0-ball enclosure $\circ d\kappa$ of $d\kappa$

Out: $e \in \circ X$ s.t. $\tilde{x} + e$ contains a unique fixed point of $\kappa$

1 $\eta \leftarrow u_{\text{prec}} \|d\|_{\infty}$; $e'_1 \leftarrow d$; $t \leftarrow 0$

2 repeat

3 $t \leftarrow t + 1$; $e_t \leftarrow e'_t + \eta$

4 **if** $\|e_t\|_{\infty} \ge 2\|e_1\|_{\infty}$ **then throw** Not Enough Precision

5 **else** $e'_{t+1} \leftarrow d + \circ d\kappa \cdot e_t$

6 **until** $e'_{t+1} \subseteq e_t$

7 return $e \leftarrow e_t$

PROPOSITION 2.2.

(i) KRAWCZYK(d, ◦dκ) either returns e or an error.

(ii) If $\kappa(\tilde{x}) \in \tilde{x} + d$ (the defect), $\circ d\kappa$ is a 0-ball extension of $d\kappa$ at $\tilde{x}$, and KRAWCZYK(d, ◦dκ) returns e, then e satisfies (2) and $\tilde{x} + e$ contains a unique fixed-point of $\kappa$.

(iii) If, as $u_{\text{prec}} = 2^{-t} \to 0$, $\text{rad}(d) = O(u_{\text{prec}})$, and $\|\circ d\kappa \cdot d^+\|_{\infty} \le \frac{1}{4}\|d^+\|_{\infty}$ where $d^+ = \{x \in X \mid \|x\|_{\infty} \le \frac{3}{2}\|d\|_{\infty}\}$, then KRAWCZYK(d, ◦dκ) eventually never raises an error and returns e with $\text{rad}(e) = O(u_{\text{prec}})$.

Proposition 2.2 does not bound the number of iterations of
KRAWCZYK, but as they only involve operations in OX at coarse
precision *u*<sub>prec</sub> independent of the target accuracy, we choose to dis-
card them in our complexity bounds. Propositions accompanying
the validated algorithms in the next sections follow this pattern:
(i) **Correctness:** if no error occurs, the inclusion property is true w.r.t. any point $x \in \mathbf{x}$ (the input) satisfying a given property;

(ii) **Lipschitz-Continuity:** For $x \in X$ satisfying a given property, as $u_{\text{prec}} \to 0$ and $\mathbf{x} \not\to \mathbf{x}$ with rad($\mathbf{x}$) = O($u_{\text{prec}}$), the algorithm eventually never fails and returns $\mathbf{y}$ s.t. rad($\mathbf{y}$) = O($u_{\text{prec}}$) (the hidden constant may depend on $\mathbf{x}$ but not on $u_{\text{prec}}$);

(iii) **Arithmetic Complexity:** an asymptotic bound on the number of FP/ball operations at precision $u_{\text{prec}}$.

# 3 Newton-Puiseux modular-numeric strategy

## 3.1 Definitions and algorithm ARNP of [22]

Let $\mathbb{K}$ be a finite extension of $\mathbb{Q}$ or $\mathbb{F}_p$, and $H(x, y) \in \mathbb{K}[x, y]$ be of degree $d$ in $y$. Up to a change of variable $x \leftarrow x + x_0$, we only consider solving $H(x, y) = 0$ for $y(x)$ above $x_0 = 0$. Under the assumption $p = 0$ (see, e.g., [4, 9, 30]) or $p > d$ [5, §IV.6], the well-known Puiseux theorem asserts:

**THEOREM 3.1 (PUISEUX, 1850).** Let $H \in \mathbb{K}[x][y]$ of degree $d$ be separable and monic. Then it has $d$ distinct roots in $\mathbb{K}((x))$:

$$S_{ij}(x) = \sum_{k=0}^{\infty} \alpha_{i,k} \zeta_{e_i}^{jk} x^{\frac{k}{e_i}}$$

with $e_1, \dots, e_s \in \mathbb{N}_{>0}$ s.t. $\sum_{i=1}^s e_i = d$, $0 \le j \le e_i - 1$, $1 \le i \le s$ and $\zeta_{e_i}$ primitive root of unity of order $e_i$. Moreover, the $\{\alpha_{i,k}\}$ belong to a finite extension of $\mathbb{K}$.

The $S_{ij}$ are the *Puiseux series* of $H$ and $e_i$ its ramification indices. One can group the Puiseux series conjugated by the ramification or the coefficient field they belong to, leading to Duval's *rational Puiseux expansions* (RPE) [8], i.e. parametrisations of the curve of the form $(\lambda T^e, S(T))$. Given $S = \sum_{\alpha \in \mathbb{Q}} c_\alpha x^\alpha$ a Puiseux series and $n \in \mathbb{Q}$, we define the truncation $[S]^n = \sum_{\alpha \le n} c_\alpha x^\alpha$, and extend it to any bivariate $H(x, y)$ coefficient-wise. The *regularity index r* of a Puiseux series $S$ of $H$ is the minimum number of terms needed to separate it from any other Puiseux series $S'$ of $H$, i.e. the least integer $N \in \mathbb{N}$ s.t. $[S]^\frac{N}{e} = [S']^\frac{N}{e}$ implies $S = S'$. $[S]^\frac{r}{e}$ is the *singular part of S* in $H$.

Algorithm ARNP [22] (Algorithm 2 below) is an improvement of Duval's rational Newton-Puiseux algorithm [8] that computes all RPEs of $H$. In a nutshell (see [23, §2.2] for more details):

* The Newton polygon $\mathcal{N}(H)$ of $H$ is the lower convex hull of its support $\text{Supp}(H) = \{(i, j) \in \mathbb{N}^2 \mid \text{coeff}(H, x^j, y^i) \neq 0\}$. Each edge $\Delta \in \mathcal{N}(H)$ is supported by a line $mi + qj = l$, $m, q, l \in \mathbb{N}$ with $q$ minimal. Its slope$^4 \lambda = \frac{m}{q} \in \mathbb{Q}_{\ge 0}$ indicates that at least one Puiseux series $S$ of $H$ has valuation $\lambda$.

* The roots of the associated characteristic polynomial $\phi_\Delta(z) = \sum_{(i,j) \in \Delta} \text{coeff}(H, x^j, y^i) z^{\frac{i-l_0}{q}}$ (with $l_0$ the leftmost point of $\Delta$) are all possible values $\xi$ s.t. $S(x) = x^{\frac{m}{q}} (\xi^{\frac{1}{q}} + o(1))$. For each $\xi$, Duval's rational Puiseux substitution $(x, y) \leftarrow (\xi^t x^q, x^m(y + \xi^s))$ (with $s, q \in \mathbb{Z}$ s.t. $sq - tm = 1$) is applied to focus on the $o(x^{\frac{m}{q}})$ tail of $S$.

* Over a finite extension $\mathbb{K}$ of $\mathbb{Q}$ or $\mathbb{F}_p$, a routine FACTOR($\mathbb{K}, \phi_\Delta$) factorises $\phi_\Delta$ into $\mathbb{K}$-irreducible factors $\phi$ of multiplicity $\mu$.

* WPT($H_2, n_1$) computes the Weierstrass factor of $H$ to precision $n_1$ whose Puiseux series are those starting with $\xi^{\frac{1}{q}}x^{\frac{m}{q}}$ (see Sect. 5).

**Algorithm 2: ARNP(H, n, π = (x, y))**

In: • $H = \sum_{i=0}^{d} A_i(x)y^i \in \mathbb{K}[x, y]$ separable and monic in $y$

• truncation order $n \in \mathbb{N}$;

• current truncated parametrisation $\pi(x, y)$.

Out: • all RPEs $R_i$ of $H$ at precision $n/e_i \ge r_i/e_i$;

• polygon tree $\mathcal{T}(H)$.

1. $B \leftarrow A_{d-1}/d$; $\pi_1 \leftarrow [\pi(x, y - B)]^n$

2. **if** $d = 1$ **then** return $(\{\pi_1(T, 0)\}, \{\})$

3. $\mathcal{T} \leftarrow \{\}$; $\mathcal{R} \leftarrow \{\}$; $H_1 \leftarrow [H(x, y - B)]^n$

4. **foreach** $\Delta$ **in** $\mathcal{N}(H_1)$ **do**

$$ \begin{array}{l} \text{\# } \phi(\xi) = 0, \quad \text{sq} - tm = 1, \quad n_1 = qn - l \\ \pi_2 \leftarrow \pi_1(\xi^l x^q, x^m(y + \xi^s)) \\ H_2 \leftarrow [H_1(\xi^l x^q, x^m(y + \xi^s))/x^l]^{n_1} \\ H_3 \leftarrow \text{WPT}(H_2, n_1) \\ (\mathcal{R}_0, \mathcal{T}_0) \leftarrow \text{ARNP}(H_3, n_1, \pi_2) \\ \mathcal{R} \leftarrow \mathcal{R} \cup \mathcal{R}_0; \text{Add } (\frac{m}{q}, \mu, \mathcal{T}_0) \text{ deg}(\phi) \text{ times to } \mathcal{T} \\ \text{end} \end{array} $$

12. end

13. return $(\mathcal{R}, \mathcal{T})$

* the “Abhyankar’s trick” (shift by $B$) extracts the part common to all Puiseux series to ensure each recursive call discovers a separation between at least two groups of Puiseux series.

* In addition, ARNP returns the polygon tree $\mathcal{T}(H)$ of $H$, defined in the next section and crucial to our modular-numeric approach.

## 3.2 Modular-numeric strategy and polygon trees

We follow the modular-numeric strategy sketched in [18] and detailed (in french) in [19]:

1. Compute the Puiseux series modulo a *well-chosen* prime number $p$

2. Follow the structure of the computed series to conduct numerical computation of the coefficients.

The criteria guaranteeing that $p$ preserves the structure of the Puiseux series (greater than $d$, not dividing the denominators of the coefficients of $H$ nor the trailing coefficient of the discriminant) are proved in [20, 21]. The structure is stored in a *polygon tree* that keeps track of the successive slopes $\lambda$ and multiplicities $\mu$, which could barely be recovered purely numerically. We use a slightly different, yet completely equivalent, definition to that of [21, §4.3].

**Definition 3.2.** A *polygon tree* $\mathcal{T}$ is a (possibly empty) list of $(\lambda_i, \mu_i, \mathcal{T}_i)$ where $\lambda_i \in \mathbb{Q}_{\ge 0}$ is a slope in the Newton polygon, $\mu_i \in \mathbb{N}_{>0}$ a multiplicity, and $\mathcal{T}_i$ a child (sub-tree) of $\mathcal{T}$. We abusively say $\lambda \in \mathcal{T}$ (resp. $(\lambda, \mu) \in \mathcal{T}$) if there is $i$ s.t. $\lambda = \lambda_i$ (resp. $(\lambda, \mu) = (\lambda_i, \mu_i)$).

In particular, the Newton polygon $\mathcal{N}(\mathcal{T})$ representing the “root” of $\mathcal{T}$ is recovered as follows. For each $\lambda = \frac{m}{q} \in \mathcal{T}$, build an edge $\Delta$ of slope $\lambda$ and width$^5 w(\Delta) = q \sum (\lambda_i, \mu_i, \mathcal{T}_i) \in \mathcal{T}$ $\mu_i$. Then $\mathcal{N}(\mathcal{T})$ is the Minkowski sum of these edges, and the multiplicity of each edge $\Delta$ is the list of $\mu_i$ s.t. $(\lambda, \mu_i, \mathcal{T}_i) \in \mathcal{T}$. Finally, given $\lambda = \frac{m}{q} \in \mathcal{T}$, we call $\ell(\mathcal{T}, \lambda)$ the integer $\ell$ s.t. $\Delta$ is supported by the line $mi + qj = l$.

<footnote>$^4$This definition actually corresponds to the opposite of the geometric slope.</footnote>

<footnote>$^5$The width of an edge $\Delta = [(i_r, j_r), (i_r, j_r)]$ is $w(\Delta) = i_r - i_l$.</footnote>
Having this polygon tree as an oracle giving the exact Newton polygon and the multiplicity structures on each edge, we will use validated numerics to compute rigorous approximations of the Puiseux series. However, inferring the reduction modulo $p$ of a complex number (in $\mathbb{Q}$) by using only converging numerical approximations is an undecidable problem. Hence, two roots $\xi_i, \xi_i' \in \mathbb{C}$ of a characteristic polynomial $\phi_\Delta$ with same multiplicity $\mu$ cannot be matched to the corresponding roots in $\overline{\mathbb{F}}_p$ computed in the first step, and therefore to the corresponding sub-trees $\mathcal{T}_i, \mathcal{T}_i'$. Consequently, we must keep together all the polynomials obtained from a Puiseux substitution with a root of $\phi_\Delta$ of multiplicity $\mu$ on one side, and all trees $\mathcal{T}_i$ s.t. $(\lambda, \mu, \mathcal{T}_i) \in \mathcal{T}$ on another side. At the next recursive call, we might have several Newton polygon candidates for the numerical polynomials: they will be sorted thanks to the filter POLFILTER in Sect. 4.1. Then for each slope $\lambda$, we might have several possible multiplicity structures for the numerical characteristic polynomials. This is solved in Sect. 4.2 using a second filter, MULTFILTER, which also computes rigorous enclosures of the multiple roots.

**Complete polygon trees.** A last difficulty to be mentioned here is that the modular criterion does not prevent the possible cancellation of Puiseux series coefficients of non-essential terms. They correspond to integer slopes for which the characteristic polynomial has a root equal to 0 modulo $p^6$. In that case, ARNP executed on $\overline{H}$ (the reduction of $H$ mod $p$) does not “see” this edge and returns a tree $\mathcal{T}(\overline{H})$ different from $\mathcal{T}(H)$ [21, §5.2]. This issue was solved in [21] by the use of generic Newton polygons. As we are running a faster algorithm, this strategy would not lead to a good complexity, and we now solve it thanks to complete polygon trees.

**Definition 3.3.** A polygon tree $\mathcal{T} = \{(\lambda_i, \mu_i, \mathcal{T}_i)\}_i$ is said complete if $\max_i \lambda_i \le \lfloor \min_i \lambda_i \rfloor \in \mathbb{N}$, and if, recursively, all $\mathcal{T}_i$ are complete.

The completion of a polygon tree $\mathcal{T} = \{(\lambda_i, \mu_i, \mathcal{T}_i)\}_i$ with $\lambda_i = \frac{m_i}{q_i}$ into a complete polygon tree $\mathcal{T}^*$ is done via the following recursive procedure. Let $\lambda = \lceil \lambda_{\min} \rceil$ and start with an empty tree $\mathcal{T}^*$. Denoting $\mathcal{T}_i^*$ the (recursively) completed tree of $\mathcal{T}_i$, do:

* For any $(\lambda_i, \mu_i, \mathcal{T}_i)$ s.t. $\lambda_i \le \lambda$, add $(\lambda_i, \mu_i, \mathcal{T}_i^*)$ to $\mathcal{T}^*$.

* If there is exactly one $(\lambda_i, \mu_i, \mathcal{T}_i)$ with $\lambda_i > \lambda$ and if additionally $\lambda_i \in \mathbb{N}$, add $(\lambda, \mu_i, \mathcal{T}_0)$ to $\mathcal{T}^*$ with $\mathcal{T}_0$ defined as the tree $\{(\lambda'_j + \lambda_i - \lambda, \mu'_j, \mathcal{T}_j^{**})\}_{j \in \mathcal{T}_i}$. If $\mathcal{T}_i = \{(\lambda'_j, \mu'_j, \mathcal{T}_j')\}_{j \in \mathcal{T}_i}$.

* Otherwise, let $\mathcal{I} = \{i \mid \lambda_i > \lambda\}$. Define $\mu = \sum_{i \in \mathcal{I}} q_i \mu_i$ and $\mathcal{T}_0 = \{(\lambda_i - \lambda, \mu_i, \mathcal{T}_i^*)\}_{i \in \mathcal{I}}$. Finally add $(\lambda, \mu, \mathcal{T}_0^*)$ to $\mathcal{T}^*$.

**THEOREM 3.4.** *Under the modular criterion, $\mathcal{T}(H)^* = \mathcal{T}(\overline{H})^*$.*

**PROOF.** (sketch) This is not difficult but slightly technical. We omit details only because of space constraints. The main idea is that [21, Thm 48] deals with everything except possibly with integer slopes. However, the Abhyankar trick of ARNP<fnref n="7" /> “removes” all integer slopes having a single root of multiplicity equal to the degree of the input polynomial. The only possible integer slopes we could miss are therefore bigger than the rightmost slope. The above strategy takes care of all possible cases. $\square$

$\mathcal{T}(\overline{H})^*$ computed by ARNP over $\overline{\mathbb{F}}_p$ is therefore a “correct” data representing the structure of the Puiseux series.

# 4 Polygon and multiplicity filters

We detail the filters mentioned in Section 3.2. For a collection $\mathcal{H} = (H_\alpha)_{\alpha \in A}$ of polynomials in $\mathbb{C}[x, y]$, and $\mathcal{L} = (\mathcal{T}_\beta)_{\beta \in B}$ a collection of polygon trees, the relation $\mathcal{H} \sim \mathcal{L}$ means that there exists a bijection $\sigma: A \to B$ s.t. $\mathcal{T}(H_\alpha)^* = \mathcal{T}_{\sigma(\alpha)}$. Each call to VALIDATEDPUISEUX is indeed given a collection of polynomials with the collection of the corresponding completed polygon trees, but $\sigma$ is unknown. Here, “collection” refers to any data structure implementing indexed sets $\mathcal{X} = \{X_\alpha\}_{\alpha \in A}$: each $X_\alpha$ is uniquely identified by an index $\alpha \in A$. This allows us to retrieve the provenance of each $X_\alpha$ regardless of any set reordering, splitting, or merging.

By inspecting the polygons and multiplicities in all $\mathcal{T} \in \mathcal{L}$, filters refine $\mathcal{H} \sim \mathcal{L}$ into a sound partition, namely a set $Q$ of pairs $(\mathcal{H}_i, \mathcal{L}_i) = (\{H_\alpha\}_{\alpha \in A_i}, \{\mathcal{T}_\beta\}_{\beta \in B_i})$ s.t. the index sets $A_i$ (resp. $B_i$) form a partition of the index set $A$ (resp. $B$), and $\mathcal{H}_i \sim \mathcal{L}_i$ for all $i$. Since VALIDATEDPUISEUX is given a collection $\mathcal{H}$ of ball polynomials $H_\alpha \in \square\mathbb{C}[x, y]$, we write $\mathcal{H} \in \mathcal{H}$ if $H_\alpha \in \mathcal{H}_\alpha$ for all $\alpha$, and we say that $Q = \{(\mathcal{H}_i, \mathcal{L}_i)\}_i$ is a sound partition of $\mathcal{H} \sim \mathcal{L}$ if for all $\mathcal{H} \in \mathcal{H}$ s.t. $\mathcal{H} \sim \mathcal{L}$, the induced partition $Q = \{(\mathcal{H}_i, \mathcal{L}_i)\}_i$ is sound.

## 4.1 Filtering by Newton polygons

POLFILTER, produces a polygon-partition of $\mathcal{H} \sim \mathcal{L}$, i.e., a sound partition where for each $i$, all $\mathcal{T}_\beta \in \mathcal{L}_i$ have the same polygon $N(\mathcal{T}_\beta)$. It is important to note that we cannot build the Newton polygon of $H_\alpha$ by checking the ball coefficients of $H_\alpha$ alone. Indeed, for $c \in c \in \square\mathbb{C}$, $0 \notin c$ is a proof of $c \neq 0$, but conversely $0 \in c$ does not imply $c = 0$ without additional information. Therefore, POLFILTER works by process of elimination, using this “one-directional information” combined with the possible polygons given by $\mathcal{L}$.

Assuming the $N_\beta = N(\mathcal{T}_\beta)$ are not all equal (otherwise return $\{(\mathcal{H}, \mathcal{L})\}$), consider the set of edges that appear in some of the $\mathcal{T}_\beta$, but not all. Among it, pick the discriminating edge, i.e., $\Delta = [(i_l, j_l), (i_r, j_r)]$ with the smallest slope $\lambda$, and maximum width in case of a tie<fnref n="8" />. Then there is a sequence<fnref n="9" /> $\Delta_1, \dots, \Delta_k$ of adjacent edges joining $(i_r, j_r)$ to $(d, 0)$ in all $N_\beta$. By checking the ball coefficients of each $H_\alpha$ on $\Delta$, POLFILTER either certifies that $\Delta \in N(H_\alpha)^* := N(\mathcal{T}(H_\alpha)^*)$, or supposes that $\Delta \in N(H_\alpha)^*$ (this guess may be wrong if the balls are too coarse). However, by reading in $\mathcal{L}$ the number of $\alpha$ s.t. $\Delta \in N(H_\alpha)^*$, it decides whether we obtained a partition according to $\Delta$ (then it recursively partitions the two subsets), or whether more precision is necessary.

**LEMMA 4.1.** Let $\mathcal{H} \sim \mathcal{L}$ and $H \in \mathcal{H}$.

(i) If $\lambda \in \mathbb{Q} \setminus \mathbb{N}$, $\Delta \in N(H)^* \Leftrightarrow (i_l, j_l) \in \text{Supp}(H)$.

(ii) If $\lambda \in \mathbb{N}$, $\Delta \in N(H)^* \Leftrightarrow \text{Supp}(H) \cap (\Delta \setminus \{(i_r, j_r)\}) \neq \emptyset$.

**PROOF.** (i) If $\Delta \in N(H)^*$, then $\Delta \in N(H)$ since $\lambda \notin \mathbb{N}$. Hence $(i_l, j_l) \in \text{Supp}(H)$. Suppose conversely that $(i_l, j_l) \in \text{Supp}(H)$. Since $N(H)^* = N_\beta$ for some $\beta$ and $(i_r, j_r)$ is already a vertex of $N_\beta$, by minimality of $\lambda$, $(i_l, j_l)$ cannot be below $\Delta$. Hence $N(H)^* = N_\beta$ has an edge $\Delta_0$ of slope $\lambda$, and $\Delta_0 = \Delta$ by maximality of $w(\Delta)$.

(ii) If $\Delta \in N(H)^*$, either $\Delta$ is equal to or extends an edge $\Delta_0$ of $N(H)$, and then $\text{Supp}(H) \cap (\Delta \setminus \{(i_r, j_r)\})$ contains at least the left endpoint of $\Delta_0$. Or $\lambda$ is a new slope introduced by the completion of $N(H)$, i.e., $\lambda = \lceil \lambda_k \rceil > \lambda_1 > \dots > \lambda_k$. Since all $N_\beta$ are complete and

<footnote n="6">⁶From the modular criterion, we cannot have two such roots on the same edge.</footnote>

<footnote n="7">⁷Note that this algorithm has been written in [22] several years after [21].</footnote>

<footnote n="8">⁸Such a tie can only happen for non-integer $\lambda$ by definition of complete trees.</footnote>

<footnote n="9">⁹The sequence is empty if $(i_r, j_r) = (d, 0)$.</footnote>
contain $\Delta_1, \dots, \Delta_k$, they must contain $\Delta$, which is a contradiction. If conversely $(i, j) \in \text{Supp}(H) \cap (\Delta \setminus \{(i_r, j_r)\})$, by minimality of $\lambda$, $(i, j)$ cannot be below $\Delta$. Hence it is connected to $(i_r, j_r)$ in $\mathcal{N}(H)^*$ by an edge $\Delta_0$ of slope $\lambda$. Since $\lambda \in \mathbb{N}$ and all polygons are complete, $\Delta_0 = \Delta$ is the edge of slope $\lambda$ joining the vertical axis to $(i_r, j_r)$. $\square$

**Algorithm 3: POLFILTER($\mathcal{H}, \mathcal{L}$)**

In: • a collection $\mathcal{H}$ of polynomials $H_\alpha \in \square\mathbb{C}[x, y]$;

• a collection $\mathcal{L}$ of polygon trees $\mathcal{T}_\beta$

Out: a polygon-partition of $\mathcal{H} \sim \mathcal{L}$.

1 if all $\mathcal{N}(\mathcal{T}_\beta)$ are equal then return $\{(\mathcal{H}, \mathcal{L})\}$

$$ \# \Delta = [(i_l, j_l), (i_r, i_r)] \text{ the discriminating edge, of slope } \lambda $$

2 $\mathcal{L}_0 \leftarrow \{\mathcal{T}_\beta \mid \Delta \in \mathcal{N}(\mathcal{T}_\beta)\}$

3 if $\lambda \in \mathbb{N}$ then

4 $\mathcal{H}_0 \leftarrow \{H_\alpha \mid \forall (i, j) \in \Delta \setminus (i_r, j_r), 0 \notin \text{coeff}(H_\alpha, x^j, y^i)\}$

5 else $\mathcal{H}_0 \leftarrow \{H_\alpha \mid 0 \notin \text{coeff}(H_\alpha, x^{j'}, y^{i'})\}$

6 if $\#\mathcal{H}_0 \neq \#\mathcal{L}_0$ then throw Not Enough Precision

7 return POLFILTER($\mathcal{H}_0, \mathcal{L}_0$) $\cup$ POLFILTER($\mathcal{H} \setminus \mathcal{H}_0, \mathcal{L} \setminus \mathcal{L}_0$)

**PROPOSITION 4.2.** Assume all $\mathcal{N}_\beta$ have same width $d$.

(i) If POLFILTER($\mathcal{H}, \mathcal{L}$) does not fail, then it returns a polygon-partition of $\mathcal{H} \sim \mathcal{L}$.

(ii) POLFILTER($-$, $\mathcal{L}$) is Lipschitz-continuous at any $\mathcal{H} \sim \mathcal{L}$.

(iii) POLFILTER($\mathcal{H}, \mathcal{L}, \lambda$) requires $O(\text{snd})$ FP/ball operations at precision $\text{u}_{\text{prec}}$, where $s = \#\mathcal{H}$ and all $H \in \square\mathbb{C}[x]_{\le n}[y]_{\le d}$.

**PROOF.** (i) Let $\mathcal{H} \in \mathcal{H}$ s.t. $\mathcal{H} \sim \mathcal{L}$. The claim is trivial if all $\mathcal{N}_\beta$ are equal. Otherwise, direction $\Leftarrow$ of Lemma 4.1(i) and (ii) guarantees that for all $H_\alpha$ s.t. $H_\alpha \in \mathcal{H}_0$, $\Delta \in \mathcal{N}(H_\alpha)^*$ holds. Line 6 ensures that all such $H_\alpha$ have been collected into $\mathcal{H}_0$, so that $\{(\mathcal{H}_0, \mathcal{L}_0), (\mathcal{H} \setminus \mathcal{H}_0, \mathcal{L} \setminus \mathcal{L}_0)\}$ is a sound partition of $\mathcal{H} \sim \mathcal{L}$ based on $\Delta$, and POLFILTER recursively refines it until the base case.

(ii) Direction $\Rightarrow$ of Lemma 4.1(i) and (ii) ensures that all $H_\alpha$ s.t. $\Delta \in \mathcal{N}(H_\alpha)^*$ will be collected into $\mathcal{H}_0$, provided that $\mathcal{H}$ is tight enough to ensure that all relevant nonzero coefficients are detected.

(iii) The next discriminating edge $\Delta'$ in the recursive call on $\mathcal{H}_0 \sim \mathcal{L}_0$ necessary lies at the left of $\Delta$ (since $\Delta$ is common to all $\mathcal{T} \in \mathcal{L}_0$). That of $(\mathcal{H} \setminus \mathcal{H}_0) \sim (\mathcal{L} \setminus \mathcal{L}_0)$ either lies above $(\lambda' > \lambda)$ or is strictly included in $\Delta$, i.e., $\lambda' = \lambda$ and $w(\Delta') < w(\Delta)$. The latter case necessarily implies $\lambda \notin \mathbb{N}$, since all polygons are complete. Therefore, each point $(i, j) \in [0, d] \times [0, n]$ is visited at most once in total, making a total of $O(\text{snd})$ 0-tests for the coefficients of $\mathcal{H}$. $\square$

## 4.2 Filtering by multiplicities

At this stage, we may assume that all trees $\mathcal{T}_\beta \in \mathcal{L}$ have the same Newton polygon $\mathcal{N}(\mathcal{T}_\beta) = \mathcal{N}$. Let $\Delta \in \mathcal{N}$ of slope $\lambda \in \frac{m}{q}$ and width $w(\Delta)$. All characteristic polynomials $\phi_{\Delta, \alpha}$ have degree $d_\Delta = \frac{w(\Delta)}{q}$. However, any $\phi \in \phi_{\Delta, \alpha}$ picked at random has $d_\Delta$ simple roots with probability 1. Therefore, any strategy to recover both roots and correct multiplicities for $\phi_{\Delta, \alpha}$ must exploit the multiplicity structures $M_\lambda(\mathcal{T}_\beta)$ along $\Delta$, defined as the multiset $\{\mu_i^{(\nu_i)}\}_i$ of multiplicities $\mu_i$ s.t. $\mathcal{T}_\beta$ has exactly $\nu_i > 0$ children labelled by $(\lambda, \mu_i)$.

In Sect. 4.2.1, we adapt the results of [3] to provide a procedure □MULTROOTS that computes the roots and multiplicities of a ball univariate polynomial based solely on the given number $k$ of distinct roots. Then in Sect. 4.2.2 we present MULTFILTER that matches the $\phi_{\Delta, \alpha}$ against the candidate multiplicity structures of the $\mathcal{T}_\beta$.

### 4.2.1 Univariate multiple root-finding.
As explained in, e.g., [3, Prop. 3.5] or [13, Thm. 4] for a (tight enough) $\phi$, knowing the number $k$ of distinct roots by some oracle completely determines the multiplicity structure of any $\phi \in \phi$ having $k$ distinct roots, and makes it possible to compute tight root enclosures.

Given $\phi \in \mathbb{C}[z]$ of degree $d$ and $k \le d$, let $S_k(\phi) \in \mathbb{C}^{(d+k) \times (2k+1)}$ denote the Sylvester-like matrix representing the linear map $(U, V) \in \mathbb{C}[z]_{\le k} \times \mathbb{C}[z]_{\le k-1} \mapsto V\phi - U\phi' \in \mathbb{C}[z]_{d+k-1}$.

**LEMMA 4.3** ([3, Lem. 2.5]). Let $\phi(z)$ have $k$ distinct roots $\xi_i$ in $\mathbb{C}$.

(i) $S_\ell(\phi)$ is injective for $0 \le \ell < k$.

(ii) The kernel of $S_k(\phi)$ has rank one. For any $(U, V) \ne 0$ in it, the roots of $U(z)$ are all simple and equal to the $\xi_i$. Furthermore, $V(\xi_i)/U'(\xi_i)$ is equal to the multiplicity $\mu_i$ of $\xi_i$ in $\phi$.

In particular, by letting $U(z) = z^k + u_{k-1}z^{k-1} + \cdots + z_0$ and $V(z) = dz^{k-1} + v_{k-2}z^{k-2} + \cdots + v_0$, $U(z)$ and $V(z)$ can be recovered from the injective but overdetermined linear system $S_{k-1} \cdot (u, v) = b$, where $b \in \mathbb{C}^{d+k-1}$ is the coefficients of $z^k\phi' - dz^{k-1}\phi$. We turn this lemma into a validated procedure. A first routine □COUNTROOTS($\phi, k$) validates a lower bound on the actual number of distinct roots of any $\phi \in \phi$ by solving the above linear system for a candidate $k$:

(1) Compute ball matrix $S_{k-1}$ and vector $b$ for $S_{k-1}$ and $b$;

(2) Compute a numerical Moore-Penrose inverse $A$ of $S_{k-1}$ and apply it onto $b$ to obtain numerical coefficients $(\tilde{u}, \tilde{v})$;

(3) Compute $\delta = |A \cdot (S_{k-1} \cdot (\tilde{u}, \tilde{v}) - b)|_\circ$ (the defect) and $\Lambda = |\text{id}_{2k-1} - A \cdot S_{k-1}|_\circ$, for the (constant) Jacobian of the Krawczyk operator $\kappa(u, v) = (u, v) - A \cdot (S_{k-1} \cdot (u, v) - b)$;

(4) Use KRAWCZYK($\delta$, ◦$\text{dk}$) with ◦$\text{dk}$ : $e \mapsto \Lambda \cdot e$, and, if it does not fail, return the obtained enclosures $U, V$ for $U, V$.

A second routine □COUNTROOTS($\phi, k$) eventually computes a $k$-validated multiple-root enclosure for $\phi$, i.e., a set $\Xi$ of pairs $(\xi_i, \mu_i) \in \square\mathbb{C} \times \mathbb{N}_{>0}$ ($1 \le i \le k$) s.t. for all $\phi \in \phi$ having exactly $k$ distinct roots$^{10}$, there exist $c \in \mathbb{C}$ and $\xi_i \in \xi_i$ s.t. $\phi(z) = c \prod_{i=1}^k (z - \xi_i)^{\mu_i}$:

(1) Compute $U, V$ using □COUNTROOTS($\phi, k$);

(2) Use a standard validated routine □ROOTS($U$) to compute enclosures $\xi_1, \dots, \xi_k$ of the (necessarily simple) roots of $U$;

(3) Check for each $i$ that the ball $V(\xi_i)/U'(\xi_i)$ contains a unique integer $\mu_i$ (fail otherwise), and return $\Xi = \{(\xi_i, \mu_i)\}_{i=1}^k$.

**PROPOSITION 4.4.**

(i) □COUNTROOTS($\phi, k$) not failing guarantees that any $\phi \in \phi$ has at least $k$ distinct roots. If □MULTROOTS($\phi, k$) succeeds, it computes a $k$-validated multiple-root enclosure for $\phi$.

(ii) □COUNTROOTS($-$, $k$) is Lipschitz-continuous at any $\phi \in \mathbb{C}[z]$ with at least $k$ distinct roots, and so is □MULTROOTS($-$, $k$) at any $\phi \in \mathbb{C}[z]$ having exactly $k$ distinct roots,

(iii) Both require $O(\text{dk}^2)$ FP/ball operations at precision $\text{u}_{\text{prec}}$.

**PROOF.** (i) If KRAWCZYK in □COUNTROOTS succeeds, then $S_{k-1}$ associated to any $\phi \in \phi$ is injective by Prop. 2.2(ii) and Thm. 2.1,

<footnote>$^{10}$It is perfectly possible that $\phi$ contain no polynomial $\phi$ with exactly $k$ distinct roots.</footnote>
hence Lemma 4.3(i) concludes. Moreover, $(u, v)$ contains a unique fixed-point of $\kappa$. But $\kappa$ being affine, its fixed-point is globally unique. Now if $\phi \in \Phi$ has exactly $k$ distinct roots, then $U, V$ in Lemma 4.3(ii) yield a fixed-point $(u, v)$ of $\kappa$. Hence $(u, v) \in (u, v)$, and □MULTROOTS, if it does not fail, computes the desired multiple-root enclosure.

(ii) For □COUNTROOTS, let $\phi$ has at least $k$ roots. Since $S_{k-1}(\phi)$ is injective by Lemma 4.3(i), its Moore-Penrose inverse is a proper left inverse. We prove the hypotheses of Proposition 2.2(ii). As $u_{\text{prec}} \to 0$, since rad$(S_{k-1})$, rad$(b) = O(u_{\text{prec}})$, that $\tilde{c}$ converges to the exact solution of the linear system and $A$ to the exact Moore-Penrose inverse, we have rad$(\delta)$, rad$(\Lambda) = O(u_{\text{prec}})$. It is clear that eventually $\| \circ d\kappa(\delta^+) \|_{\infty} = \|\Lambda \cdot \delta^+\|_{\infty} \le \frac{1}{4}\|\delta^+\|_{\infty}$ where $\delta^+ = \{e \in \mathbb{R}^{2k-1} \mid \|e\|_{\infty} \le \frac{3}{2}\|\delta\|_{\infty}\}$. The Lipschitz-continuity of □MULTROOTS essentially boils down to that of □ROOTS.

(iii) The arithmetic complexity is dominated by the computation of $A$ and $\Lambda$, in $O(dk^2)$ (KRAWCZYK only involves ball operations in coarse precision $u'_{\text{prec}}$). Finally □ROOTS($u$) runs in arithmetic complexity $O^*(k)$ (see, e.g., [17, Thm. 2.1.1]) □

### 4.2.2 The filtering-by-multiplicity algorithm.
Exploiting routines □COUNTROOTS and □MULTROOTS above, MULTFILTER($\mathcal{H}$, $\mathcal{L}$, $\lambda$) returns a multiplicity-partition of $\mathcal{H} \sim \mathcal{L}$ along $\Delta \in \mathcal{N}$ of slope $\lambda$, i.e., a sound partition s.t. all $\mathcal{T}_{\beta} \in \mathcal{L}$ have the same $M_{\lambda}(\mathcal{T}_{\beta})$. It also returns the multiple-root enclosures $\Xi_{\alpha}$ associated to each $\phi_{\Delta, \alpha}$. The filtering strategy essentially relies on the possible values of $k$ (number of distinct roots) for the $\phi_{\Delta, \alpha}$. Indeed, for $M_{\lambda}(\mathcal{T}_{\beta}) = \{\mu_i^{(v_i)}\}_i$, let $k_{\beta} = \sum_i v_i$ and consider the sorted multiset $\mathcal{K}$ of the $k_{\beta}$,

$$ \mathcal{K} = \{k_1^{(\ell_1)} < \dots < k_s^{(\ell_s)}\}, \quad (\ell_1 + \dots + \ell_s = \#\mathcal{L}). \quad (3) $$

For each $i$ and $\alpha$, recall that if □COUNTROOTS($\phi_{\Delta, \alpha}$, $k_i$) terminates without error, then $k_i \le k_{\alpha}$ where $k_{\alpha}$ is the actual number of distinct roots. Hence the number of $\alpha$ that pass this test must be exactly $\sum_{j \ge i} \ell_j$. If the $k_{\alpha}$ are all correctly recovered, then □MULTROOTS eventually computes roots and (necessarily correct) multiplicities, which yields the desired multiplicity-partition.

<caption>**Algorithm 4: MULTFILTER($\mathcal{H}$, $\mathcal{L}$, $\lambda$)**</caption>

In:
* a collection $\mathcal{H}$ of polynomials $H_\alpha \in \square\mathbb{C}[x, y]$;
* a collection $\mathcal{L}$ of trees $\mathcal{T}_\beta$ with common polygon $\mathcal{N}$
* $\lambda$ the slope of an edge $\Delta \in \mathcal{N}$.
Out:
* a multiplicity-partition of $\mathcal{H} \sim \mathcal{L}$;
* multiple-root enclosures $\Xi_\alpha$ for each $\phi_{\Delta, \alpha}$.

\# $\mathcal{K} = \{k_1^{(\ell_1)} < \dots < k_s^{(\ell_s)}\}$ as in Equation (3)

1 for $i = 1 \dots s$ do
2 &nbsp;&nbsp; for $H_\alpha \in \mathcal{H}$ do
3 &nbsp;&nbsp;&nbsp;&nbsp; try: $\square\text{COUNTROOTS}(\phi_{\Delta, \alpha}, k_i)$; $k_\alpha \leftarrow k_i$
4 &nbsp;&nbsp; end
5 &nbsp;&nbsp; if $\#\{\alpha | k_\alpha = k_i\} \ne \sum_{j \ge i} \ell_j$ then
6 &nbsp;&nbsp;&nbsp;&nbsp; throw Not Enough Precision
7 end
8 for $H_\alpha \in \mathcal{H}$ do $\Xi_\alpha \leftarrow \square\text{MULTROOTS}(\phi_{\Delta, \alpha}, k_\alpha)$
9 \# $M_\alpha$ = multiplicity structure given by $\Xi_\alpha$
10 return (partition of $\mathcal{H} \sim \mathcal{L}$ by $M_\alpha, M_\beta, \{\Xi_\alpha\}_\alpha$)

**PROPOSITION 4.5.** Let $\mathcal{L}$ be s.t. all $\mathcal{T}_{\beta} \in \mathcal{L}$ have the same Newton polygon $\mathcal{N}(\mathcal{T}_{\beta}) = \mathcal{N}$, and let $\Delta \in \mathcal{N}$ be of slope $\lambda = \frac{m}{q}$.

(i) If MULTFILTER($\mathcal{H}$, $\mathcal{L}$, $\lambda$) does not fail, then it returns a multiplicity-partition of $\mathcal{H} \sim \mathcal{L}$ w.r.t. $\Delta$ and validated-multiple root enclosures for the characteristic polynomials $\phi_{\Delta, \alpha}$.

(ii) MULTFILTER($-$, $\mathcal{L}$, $\lambda$) is Lipschitz-continuous at any $\mathcal{H} \sim \mathcal{L}$.

(iii) Algorithm MULTFILTER($\mathcal{H}$, $\mathcal{L}$, $\lambda$) requires $O(\text{sd}_{\Delta}^3)$ FP/ball operations at precision $u_{\text{prec}}$, where $s = \#\mathcal{H}$ and $d_{\Delta} = w(\Delta)/q$.

**PROOF.** (i) Let $\mathcal{H} \in \mathcal{H}$ with $\mathcal{H} \sim \mathcal{L}$. By Proposition 4.4(i), for each $k_i$, □COUNTROOTS($\phi_{\Delta, \alpha}$, $k_i$) necessarily fails if $\phi_{\Delta, \alpha}$ has strictly less than $k_i$ distinct roots. Hence at most $\sum_{j \ge i} \ell_j$ polynomials $\phi_{\Delta, \alpha}$ pass the test, and having equality in line 5 means that they have all been correctly identified. At the end, if no error occurs, the correct $k_{\alpha}$ has been found for each $\phi_{\Delta, \alpha}$, and the rest follows by the correctness of □MULTROOTS.

(ii) is a straightforward consequence of Prop. 4.4(ii).

(iii) To reach the claimed complexity bound, □COUNTROOTS must be implemented with the following improvement. For a given $\phi \in \square\mathbb{C}[z]$ of degree $d_{\Delta}$, as noted in [33, §4.2.1], the QR factorisation<fnref n="11" /> of $S_{k_i-1}(\phi)$ requires only $O(d_{\phi}k_i(k_i-k_{i-1}))$ additional operations when reusing that of $S_{k_{i-1}-1}(\phi)$, and similarly for $\Lambda$. Hence, the overall cost of all calls to □COUNTROOTS drops to $O(\text{sd}_{\Delta}^3)$. □

## 5 Validated Weierstrass Preparation Theorem

After each Puiseux change of variable w.r.t. a root $\xi$ in the Newton-Puiseux algorithm, we get a new polynomial $H(x, y) \in \square[x]_{\le n}[y]_{\le d}$ that satisfies $v_y(H(0, y)) = \mu > 0$, where $\mu$ is the multiplicity of $\xi$ and $v_y$ the valuation in $y$. Extracting the part of $H$ that contains the Puiseux series above $\xi$ is exactly the purpose of the Weierstrass preparation theorem: it computes the unique decomposition $H = [WP]^n$ s.t. $W$ (the Weierstrass polynomial) in monic of degree $\mu$ in $y$ with $W(0, y) = y^{\mu}$, and $P(0, 0) \ne 0$ (the discarded part containing the other Puiseux series).

Calling $X_1 = x\square[x]<n[y]<\mu, \tilde{X}_1 = y^{\mu}+X_1, X_2 = \square[x]_{\le n}[y]_{\le d-\mu}$, and $Y = y^{\mu}X_2 + X_1$, we have $W \in \tilde{X}_1, P \in X_2$ and $H \in Y$. We seek the unique zero of $f: \tilde{X}_1 \times X_2 \to Y$, $(W, P) \mapsto [WP]^n - H$ with $W(0, y) = y^{\mu}$. This is a particular case of Hensel lifting, and we here adapt the numerically validated algorithm of [2] to this setting<fnref n="12" />.

We first assume the existence of a numerical Hensel lifting routine HENSEL to compute approximate numerical factors $\tilde{W}, \tilde{P}$ as well as approximations $\tilde{S}, \tilde{T}$ of the Bézout cofactors $S, T$ at precision $n$ in $x$, i.e., $[SW + TP]^n = 1$ (see [2, §4.1]). Then a Krawczyk operator $\kappa: \tilde{X}_1 \times X_2 \to \tilde{X}_1 \times X_2$ is built by constructing an explicit approximate inverse $A$ of $df(\tilde{W}, \tilde{P})$, using $\tilde{W}, \tilde{P}, \tilde{S}, \tilde{T}$ instead of $W, P, S, T$ in the following lemma (see [2, Lem. 7]).

**LEMMA 5.1.** The differential $df(W, P) : (E_1, E_2) \in X_1 \times X_2 \mapsto [PE_1 + WE_2]^n \in Y$ where $P(0, 0) \ne 0$ is invertible with inverse given by

$$ df(W, P)^{-1}(E) = \begin{pmatrix} [TE] \% W]^n \\ [SE + ([TE] \div W) P]^n \% y^{d-\mu+1} \end{pmatrix}, $$

where $A \div B$ and $A \% B$ denote the quotient and remainder of the Euclidean division of $A$ by $B$.

<footnote n="11">¹¹needed to compute the Moore-Penrose inverse $A$ in □COUNTROOTS.</footnote>

<footnote n="12">¹²The norm-wise validation in [2] (as opposed of the component-wise validation using Krawczyk's method here) produced coarser enclosures.</footnote>
Algorithm □WPT essentially applies Krawczyk's method on $κ$. It requires a validated routine □DivRem($A, B, n$) computing $\lceil A \div B \rceil^n$ and $\lceil A \% B \rceil^n$, using Krawczyk's method again (see [2, §3]).

<caption>**Algorithm 5:** □WPT($H, μ, n$)</caption>

**In:**
* a polynomial $H \in \mathbb{C}[x]_{\le n}[y]_{\le d}$;
* the $y$-valuation $μ \in \mathbb{N}_{>0}$;
* the $x$-truncation order $n \in \mathbb{N}$.

**Out:**
$W \in \square \tilde{X}_1$ the Weierstrass polynomial of $H$.
1. $(\tilde{W}, \tilde{P}, \tilde{S}, \tilde{T}) \leftarrow \text{HENSEL}(\tilde{H}, y^μ, [\tilde{H} \% x] \div y^μ, n)$
2. $E \leftarrow [\tilde{W}\tilde{P} - H]^n$, $E' \leftarrow [\tilde{S}\tilde{W} + \tilde{T}\tilde{P} - 1]^n$
3. $(Q, R) \leftarrow \square \text{DivRem}(\tilde{T}E, \tilde{W}, n)$
4. $(D_1, D_2) \leftarrow (|R|_\circ, |\lceil \tilde{S}E + \tilde{Q}\tilde{P} \rceil^n \% y^{d-μ+1}|_\circ)$
5. $\circ dK \leftarrow \text{proc}(E_1, E_2)$:
6. $\quad (Q', R') \leftarrow \square \text{DivRem}(E'E_1 + 2\tilde{T}E_1E_2, \tilde{W}, n)$
7. $\quad \text{return } (R', \lceil E'E_2 + 2\tilde{S}E_1E_2 + \tilde{G}Q' \rceil^n \% y^{d-μ+1})$
8. end
9. $(E_1, E_2) \leftarrow \text{KRAWCZYK}((D_1, D_2), \circ dK)$
10. $\text{return } W \leftarrow \tilde{W} + E_1$

**PROPOSITION 5.2.**

(i) If □WPT($H, μ, n$) does not fail, then it computes $W \in \square \tilde{X}_1$ containing the Weierstrass polynomial $W$ of any $H \in H$.

(ii) □WPT(-, μ, n) is Lipschitz-continuous at $H$ if $v_y(H(0, y)) = μ$.

(iii) Algorithm □WPT($H, μ, n$) with $H \in \mathbb{C}[x]_{\le n}[y]_{\le d}$ requires $\mathcal{O}(nd)$ FP/ball operations at precision $u_{\text{prec}}$.

**PROOF.** (i) follows from the correctness of KRAWCZYK (Prop. 2.2(ii)). By Lemma 5.1, $(D_1, D_2)$ is an enclosure of the defect $A \cdot f(\tilde{W}, \tilde{P})$ and $\circ dK$ is a 0-ball enclosure of $dk$ at $(\tilde{W}, \tilde{P})$ by [2, Lem. 8].

(ii) We prove the hypotheses of Proposition 2.2(ii). Since, for fixed $H$ s.t. $v_y(H(0, y)) = μ$, the call to HENSEL is essentially a SLP, $\tilde{W}, \tilde{P}, \tilde{S}, \tilde{T}$ converge to $W, P, S, T$ in $\mathcal{O}(u_{\text{prec}})$. Therefore, rad($E$), rad($E'$) = $\mathcal{O}(u_{\text{prec}})$, and rad($D_1$), rad($D_2$) = $\mathcal{O}(u_{\text{prec}})$. Also, when $u_{\text{prec}} \to 0$, $\|\circ dK(D^+)\|_{\infty} \le \frac{1}{4}\|D^+\|_{\infty}$, where $D^+ = \{(E_1, E_2) \in X_1 \times X_2 | \frac{3}{2}\|(E_1, E_2)\|_{\infty} \le \|(D_1, D_2)\|_{\infty}\}$, holds because in $\mathcal{O}(K, E_1)$ and $E_2$ are multiplied<fnref n="13" /> together or with $E'$, which are all $\mathcal{O}(u_{\text{prec}})$.

(iii) Numerical HENSEL costs $\mathcal{O}^{-}(nd)$ FP operations at precision $u_{\text{prec}}$. Lines 2-4 involve additions, multiplications and Euclidean divisions of polynomials in $\mathbb{C}[x]_{\le n}[y]_{\le O(d)}$, hence a cost of $\mathcal{O}^{-}(nd)$ ball operations at precision $u_{\text{prec}}$. Finally, $\circ dK$ and KRAWCZYK only involve operations at coarser precision $u_{\text{prec}}$. □

# 6 The validated Newton-Puiseux algorithm

Given $H \in \mathbb{C}[x, y]$ with $\mathcal{T} = \mathcal{T}(H)^*$, Algorithm VALIDATEDPUISEUX is initially called on $\mathcal{H} = \{H\}$ s.t. $H \in H$, $\Pi = \{(x, y)\}$, and $\mathcal{L} = \{\mathcal{T}\}$. If, however, for some edge $\Delta \in \mathcal{N}(\mathcal{T})$ of slope $\lambda$, the characteristic polynomial $\phi_\Delta$ has several roots $\xi_\alpha$ of same multiplicity $\mu$, then we cannot decide yet which polynomial $H_\alpha$ obtained from

$H$ by the Puiseux substitution<fnref n="14" /> with $\xi_\alpha$, is associated to which $\mathcal{T}_\beta$ s.t. $(\lambda, \mu, \mathcal{T}_\beta) \in \mathcal{T}$. Hence VALIDATEDPUISEUX is recursively called with collections $\mathcal{H} \sim \mathcal{L}$ of polynomials and trees.

Using POLFILTER, we obtain a first partition $Q_1$, s.t. in each $(\mathcal{H}_1, \mathcal{L}_1) \in Q_1$, all $\mathcal{T} \in \mathcal{L}_1$ have the same Newton polygon $N$. We write $\lambda \in \mathcal{L}_1$ if $N$ has an edge $\Delta$ of slope $\lambda$. For each such $\lambda$, MULTFILTER computes a second partition $Q_2$ s.t. in each $(\mathcal{H}_2, \mathcal{L}_2) \in Q_2$, the multiplicity structure $M_\lambda(\mathcal{T})$ is the same for all $\mathcal{T} \in \mathcal{L}_2$. Again, we say $(\lambda, \mu) \in \mathcal{L}_2$ if $\mu \in M_\lambda(\mathcal{T})$. This allows for the safe computation of the roots of the associated characteristic polynomials, and, after the Puiseux substitution and a call to □WPT, for a recursive call of the algorithm.

<caption>**Algorithm 6:** VALIDATEDPUISEUX($\mathcal{H}, \Pi, \mathcal{L}, n$)</caption>

**In:**
* a collection $\mathcal{H}$ of polynomials $H_\alpha \in \square \mathbb{C}[x, y]$;
* a collection $\Pi$ of associated parametrisations $\pi_\alpha$;
* a collection $\mathcal{L}$ of polygon trees $\mathcal{T}_\beta$.

**Out:**
the list of all RPEs associated to the $H_\alpha$.
1. $\mathcal{R} \leftarrow \{\}$
2. for $H_\alpha \in \mathcal{H}$ do
3. &nbsp;&nbsp;&nbsp;&nbsp;\# $B = \text{coeff}(H_\alpha, y^{d-1})/d \in \square \mathbb{C}[x]$ 
4. &nbsp;&nbsp;&nbsp;&nbsp;$\pi_\alpha \leftarrow \pi_\alpha(x, y - B)$; $H_\alpha \leftarrow [\mathcal{H}_\alpha(x, y - B)]^n$
5. end
6. for $\mathcal{H}_1, \mathcal{L}_1 \in \text{POLFILTER}(\mathcal{H}, \mathcal{L})$ and $\lambda \in \mathcal{L}_1$ do
7. &nbsp;&nbsp;&nbsp;&nbsp;($Q, \Xi) \leftarrow \text{MULTFILTER}(\mathcal{H}_1, \mathcal{L}_1, \lambda)$
8. &nbsp;&nbsp;&nbsp;&nbsp;for $(\mathcal{H}_2, \mathcal{L}_2) \in \mathcal{Q}$ and $\mu$ s.t. $(\lambda, \mu) \in \mathcal{L}_2$ do
9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;\# $\lambda = m/q$, $\ell = \ell(\mathcal{L}_2, \lambda)$, $n_1 = qn - \ell$
10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \mathcal{H}' \leftarrow \{[\mathcal{H}_\alpha(x^q, x^m(y + \xi^\frac{1}{q}))/x^\ell]^n \}_{n_1} $$, $\mathcal{H}_\alpha \in \mathcal{H}_2, (\xi, \mu) \in \Xi_\alpha \} $$
11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \Pi' \leftarrow \{\pi_\alpha(x^q, x^m(y + \xi^\frac{1}{q})), \pi_\alpha \text{ s.t. } \mathcal{H}_\alpha \in \mathcal{H}_2, (\xi, \mu) \in \Xi_\alpha \} $$
12. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Apply □WPT(-, $\mu$, $n_1$) to each $H \in \mathcal{H}'$
13. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \mathcal{L}' \leftarrow \{\mathcal{T} \mid (\lambda, \mu, \mathcal{T}) \in \mathcal{T}_\beta, \mathcal{T}_\beta \in \mathcal{L}_2\} $$
14. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \mathcal{R} \leftarrow \mathcal{R} \cup \text{VALIDATEDPUISEUX}(\mathcal{H}', \Pi', \mathcal{L}', n_1) $$
15. end
16. return $\mathcal{R}$

**PROPOSITION 6.1.**

(i) If VALIDATEDPUISEUX($\mathcal{H}, \Pi, \mathcal{L}, n$) does not fail, then it computes a collection $\mathcal{R}$ of ball RPES s.t. for all $\Pi = \{\pi_\alpha\}_{\alpha \in \Pi}$, all $\mathcal{H} = \{H_\alpha\}_{\alpha \in \mathcal{H}}$ with $\mathcal{H} \sim \mathcal{L}$, all $H_\alpha \in \mathcal{H}$ and all RPE $R$ of $H_\alpha$, $\lceil \pi_\alpha \circ R \rceil^n$ is contained in exactly one RPE of $\mathcal{R}$.

(ii) VALIDATEDPUISEUX(-, -, $\mathcal{L}, n$) is Lipschitz-continuous at any $\mathcal{H} \sim \mathcal{L}$ and $\Pi$.

(iii) It requires $\mathcal{O}^{-}(s(nd^2+d^3))$ FP/ball operations at precision $u_{\text{prec}}$, where $s = \#\mathcal{H}$ and $d$ the common $y$-degree of the $H_\alpha$.

**PROOF.** (i) and (ii) are direct consequences of the discussion above combined with the correctness and Lipschitz-convergence

<footnote n="13">¹³Line 6 contains a call to DivRem, but as shown in [2, §4.2.2], it can be avoided by precomputing (outside $\circ dK$) the quotient of $y^{d+\mu-2}$ by $\tilde{W}$, so that $\circ dK$ really involves only additions and multiplications of polynomials.</footnote>

<footnote n="14">¹⁴Note that we here simply use the classical Puiseux substitution instead of Duval's rational one, since the main advantage of the latter is to minimise field extension degrees when working in exact arithmetic.</footnote>
of the filters and □WPT stated in Propositions 4.2, 4.5 and 5.2. The remaining operations, i.e., the shifts in Abhyankar's trick and in the Puiseux substitutions, are simply SLPs executed in ball arithmetic, so that they are also correct and Lipschitz-continuous.

(iii) We prove the claimed complexity bound inductively.

* First, each shift of $y$ by $B$ in $\pi_{\alpha}$ and $H_{\alpha}$ can be done in $\mathcal{O}(nd)$ ball operations, which sums to $\mathcal{O}(snd)$ [1, Pb. 2.6].

* The same bound holds for POLFILTER by Prop. 4.2(iii).

From now on, we may assume w.l.o.g. that all $\mathcal{T} \in \mathcal{L}$ have the same polygon $N$. Let $\Delta \in \mathcal{N}$ be of slope $\lambda = \frac{m}{q}$ and width $w(\Delta)$. Call $d_{\Delta} = \frac{w(\Delta)}{q}$ the degree of the characteristic polynomials, and note that $\sum_{\Delta} qd_{\Delta} = \sum_{\Delta} w(\Delta) = d$.

* By Proposition 4.5(iii), MULTFILTER costs $\mathcal{O}(sd_{\Delta}^3)$ ball operations. Summing over $\Delta \in \mathcal{N}$ yields a bound of $\mathcal{O}(sd^3)$.

Assume again w.l.o.g. that all $\mathcal{T} \in \mathcal{L}$ have the same multiplicity structure: $M_{\lambda}(\mathcal{T}) = M = \{\mu_i^{(v_i)}\}$ with $\sum_{i} v_i \mu_i = d_{\Delta}$. Let $\mu^{(v)} \in M$.

* For each $H_{\alpha}$ and $(\xi, \mu) \in \Xi_{\alpha}$, the cost of the Puiseux substitution is $\mathcal{O}(n_1 d) \le \mathcal{O}(ndq)$ (essentially a shift of $y$ by $\xi$). Similarly, the cost of □WPT is $\mathcal{O}(ndq)$ by Proposition 5.2(iii). Repeating these operations $sv$ times (for each $H_{\alpha}$ and $\xi$) and summing over $\Delta, \mu$,

$$ \sum_{\Delta} \sum_{\mu} ndqsv = SND \sum_{\Delta} q \sum_{\mu} v \leq SND \sum_{\Delta} qd_{\Delta} = SND^2, $$

gives a total contribution of $\mathcal{O}(snd^2)$ to the complexity bound.

* Finally, the recursive call to VALIDATEDPUISEUX costs $\mathcal{O}(s'(qnd'^2+d'^3))$ by induction hypothesis, with $s' = \#\mathcal{L}' = sv$ and $d' = \mu$. Noting that $\sum_{\Delta} \sum_{\mu} qs'd' = s \sum_{\Delta} q \sum_{\mu} v\mu = sd$, we have

$$ \sigma\left(\sum_{\Delta} \sum_{\mu} s'(qnd'^2 + d'^3)\right) \subseteq \sigma(s(nd^2 + d^3)). \quad \square $$

**PROOF OF THEOREM 1.** From [21, Prop. 56], one can compute a prime $p$ of size $\mathcal{O}(d\tau)$ that has local good reduction. Compute the Puiseux series of $F$ modulo $p$ in $\mathcal{O}(D^4\tau)$ via e.g. [23] and use Sect. 3.2 to get the polygon tree. Conclude with Proposition 6.1, since $n$ can be the valuation of the discriminant [22, 23], i.e. $O(D^2)$. $\square$

**PROOF OF THEOREM 2.** Compute the discriminant in $\mathcal{O}(D^2\tau)$ [29, Cor 11.20], then its square-free factorisation in $\mathcal{O}(D^6 + D^5\tau)$ (Yun's algorithm, combined with the bit complexity of the gcd computation [29, Cor. 11.14] and the size $(D\tau, D^2)$ of the discriminant). From [21, Prop. 56], one can take $p$ of size $O(\log(D\tau))$ that will not change the multiplicities of the roots of the discriminant in an average of at most 2 iterations. For each squarefree factor, the Puiseux series above each of its root modulo $p$ to get a set of polygon trees via Section 3.2. It takes less than $\mathcal{O}(D^3)$ operations via [23]. Finally, isolate its numerical roots for a total cost of $\mathcal{O}(D^6 + D^5\tau)$ [7, Prop 2.23], then start VALIDATEDPUISEUX with this set of polygon trees and the set of bivariate polynomials one get by shifting the $x$ variable from each numerical root of the factor. This is less $\mathcal{O}(D^5\eta)$ from Proposition 6.1 (the sum of $s$ is bounded by the degree of the discriminant, thus $O(D^2)$, but so is the sum of the $n$ $s$ since $n$ can be bounded by the valuation of the discriminant [22, 23]). $\square$

# 7 Implementation and numerical examples

We provide implementation of VALIDATEDPUISEUX and sub-algorithms in an open source package NUMPUISEUX.JL<fnref n="15" />. It is based on NEMO.JL [10], itself based on FLINT/ARB [27]

An element of the ring $\square\mathbb{C}[x]$ is hence represented by the type ACBPOLYRINGELEM and an element of the ring $\square\mathbb{C}[x, y]$ by the type GENERIC.POLY{ACBPOLYRINGELEM}. Concerning implementation of WPT and VALGCDROOT, we would like to use complex floating point polynomial in order to do plain floating point computations, however, complex floating points in ABSTRACTALGEBRA.JL don't have the ring structure so we can't define this types of polynomial. As a result, plain floating points computations are replaced by ball computations and extractions of center. We now present a few examples computed with our Julia package that illustrates some potentially surprising difficulties, but also its efficiency (in comparison with symbolic computations) in most case. All of them are available in a folder examples on the repository.

**An apparently simple example.** The table 1 illustrate an example where we consider a polynomial $P$ with the following characteristics: $P$ has degree 5 in $y$ (36 in $x$) and its Puiseux series are Taylor series with regularity index at most 10. Coefficients of these Puiseux series are not particularly close, defined as rational numbers with less than 5 digits for their numerator and denominator. Despite looking like a really basic example, this polynomial required 100 bits of precision to compute Puiseux series. The reason why is associated to the first time we call algorithm □WPT: here, the Bézout cofactors needed for Hensel lighting have a pole at $\approx 0.15$. The main consequence is that when the degree in $x$ of the lifting increases, the size of the coefficients increase consequently: as the condition number of the Weierstrass preparation theorem is high if coefficients of Bézout cofactors are high, this example have a bad condition number. The situation was even worse while using the Hensel lifting of [2] (we then needed 150 bits of precision to compute the Puiseux series). This is why we switched to Krawczyk's method in this paper and implementation.

**Accuracy of □COUNTROOTS and □MULTROOTS.** We encapsulate lines 2 and 3 of MULTFILTER into a function in our implementation. We run MULTFILTER with two polynomials having the same Newton polygon but different multiplicity structure. The first one have a multiplicity structure [1, 1, 4] with two roots of distance less than $10^{-i}$ and the second one of multiplicity structure [1, 5]. Using 256 bits of precision, our validated Newton-Puiseux algorithm succeeds up to $i = 18$. Even if already satisfactory, we realised that the failure when $i > 18$ is due to □MULTROOTS, while the □COUNTROOTS part of the algorithm succeeds until $i = 38$. We believe that the two step validations of our current algorithm (first □COUNTROOTS, then □MULTROOTS) explains this “bad” behaviour of the whole filter in comparison to the sole □COUNTROOTS step. We intend to work on a a-posteriori validation of the two steps together to improve this.

**Random examples.** We present in table 2, a set of examples having same polygon tree and where coefficients are randomly generated. Each Puiseux series of each examples are Taylor series and 7 of them are of regularity index 3 and the last one is of regularity index

<footnote n="15">¹⁵The package is available here: https://gitlab.univ-lille.fr/adrien.poteaux/numpuiseux.jl</footnote>
2. Over 20 examples, we remark that 14 of them need only 53 bits to compute Puiseux series. Otherwise, example that cannot pass with 53 bits of precision have roots that are near.

**Computing Puiseux series over several critical points.** We finally illustrate Theorem 2 with one example having nearby singularities. This example is defined using so-called Mignotte's polynomials [14]: $M_{a,d}(x) = x^d - 2(a * x - 1)^2$ is proved to have two real roots distant by less than $2a^{-d/2}$. Here we consider $P = (y^3 - M_{6,10}(x))(y^3 - M_{3,10}(x)) + x^5 y^2$, and are interested at Puiseux series above the roots of its discriminant. Over $\mathbb{Q}$, the squarefree factorisation of the discriminant is $x^{15}R(x)$, with $R$ of degree 39 having an irreducible factor $P$ of degree 30. We focused on this factor $P$. Modular computation can be done with $p = 101$, which is proved to be a global good prime [21], so that this part is not a problem. Puiseux series are not complex (ramification 2, regularity index at most 1), but working with an algebraic number of degree 30 (that has a defining polynomial with coefficient-size bigger than $2^{46}$) generate a huge coefficient growth that make symbolic computation of a few terms difficult (more than 2 minutes to compute 7 terms on a laptop and the Maple implementation for instance; remember that the total degree is less than 10). However, numerical computations can be run really quickly, and work with a low precision: just computing the roots of $P$ with Nemo (thus Arb) requires 100 bits of precision. But then, above 29 of the 30 roots, one can compute the Puiseux series with only 53 bits of precision. For the last one, we had to increase the precision up to 75: this root is at distance less than $2^{-20}$ of another root of $R$ (that happens not to be a root of $P$).

## References

[1] Dario Bini and Victor Pan. *Polynomial and matrix computations. Vol. 1. Fundamental algorithms*. Birkhäuser, 1994.

[2] Florent Bréhard, Jasmin Krüger, Adrien Poteaux, and Arthur Vinciguerra. Validated numerical Hensel lifting. *Mathematics in Computer Science*, 2025. In press.

[3] Florent Bréhard, Adrien Poteaux, and Léo Soudant. Validated root enclosures for interval polynomials with multiplicities. In *Proceedings of the 2023 International Symposium on Symbolic and Algebraic Computation*, ISSAC '23, pages 90–99, New York, NY, USA, 2023. Association for Computing Machinery.

[4] E. Brieskorn and H. Knörrer. *Plane Algebraic Curves*. Birkhäuser, 1986.

[5] Claude Chevalley. *Introduction to the Theory of Algebraic Functions of One Variable*, volume 6 of *Mathematical Surveys*. American Mathematical Society, 1951.

[6] Alexander L. Chistov. Polynomial complexity of the Newton-Puiseux algorithm. In *Mathematical Foundations of Computer Science 1986*, pages 247–255, London, UK, 1986. Springer-Verlag.

[7] Daouda Niang Diatta, Sény Diatta, Fabrice Rouillier, Marie-Françoise Roy, and Michael Sragaloff. Bounds for polynomials on algebraic numbers and application to curve topology. *Discrete & Computational Geometry*, 67(3):631–697, 2022.

[8] Dominique Duval. Rational Puiseux expansions. *Compositio Math.*, 70(2):119–154, 1989.

[9] Martin Eichler. *Introduction to the Theory of Algebraic Numbers and Functions*. Academic Press, 1966.

[10] Claus Fieker, William Hart, Tommy Hofmann, and Fredrik Johansson. Nemo/hecke: Computer algebra and number theory packages for the julia programming language. In *Proceedings of the 2017 ACM on International Symposium on Symbolic and Algebraic Computation*, ISSAC '17, pages 157–164, New York, NY, USA, 2017. ACM.

[11] Alexandre Guillemot and Pierre Lairez. Validated numerics for algebraic path tracking. In *Proceedings of the 2024 International Symposium on Symbolic and Algebraic Computation*, ISSAC '24, pages 36–45, New York, NY, USA, 2024. Association for Computing Machinery.

[12] F. Johansson. Arb: efficient arbitrary-precision midpoint-radius interval arithmetic. *IEEE Transactions on Computers*, 66:1281–1292, 2017.

[13] Kurt Mehlhorn, Michael Sragaloff, and Pengming Wang. From approximate factorization to root isolation with application to cylindrical algebraic decomposition. *Journal of Symbolic Computation*, 66:34–69, 2015.

<table><caption>**Table 1:** Error bounds on coefficients of $P$ and an approximation of the conditioning number $κ$ depending on $n$ the index to which we lift. Coefficients $c_1, c_2, c_3$ and $c_4$ corresponds respectively to $\text{coeff}(P, x^{n/2}, y^0)$, $\text{coeff}(P, x^n, y^0)$, $\text{coeff}(P, x^{n/2}, y^3)$ and $\text{coeff}(P, x^n, y^3)$. $\rho \approx 0.15$ the pole of Bézout cofactors.</caption><thead><tr><th>n</th><th>κ</th><th>c<sub>1</sub>ρ<sup>n</sup></th><th>c<sub>2</sub>ρ<sup>n</sup></th><th>c<sub>3</sub>ρ<sup>n</sup></th><th>c<sub>4</sub>ρ<sup>n</sup></th><th>c<sub>1</sub>ρ<sup>n/2</sup></th><th>c<sub>2</sub>ρ<sup>n/2</sup></th><th>c<sub>3</sub>ρ<sup>n/2</sup></th><th>c<sub>4</sub>ρ<sup>n/2</sup></th></tr></thead><tbody><tr><td>1</td><td>2.7e1</td><td>0.0</td><td>0.0</td><td>1.3e-39</td><td>1.3e-29</td><td>0.0</td><td>0.0</td><td>8.8e-39</td><td>8.8e-29</td></tr><tr><td>10</td><td>1.1e11</td><td>8.6e-34</td><td>2.0e-27</td><td>6.0e-32</td><td>4.5e-26</td><td>1.1e-29</td><td>2.6e-23</td><td>7.9e-28</td><td>5.9e-22</td></tr><tr><td>20</td><td>1.4e20</td><td>1.1e-35</td><td>1.1e-25</td><td>2.6e-34</td><td>6.6e-25</td><td>1.9e-27</td><td>1.9e-17</td><td>4.5e-26</td><td>1.1e-16</td></tr><tr><td>30</td><td>8.5e28</td><td>1.2e-38</td><td>5.8e-25</td><td>1.0e-37</td><td>2.6e-24</td><td>2.8e-26</td><td>1.3e-12</td><td>2.3e-25</td><td>5.9e-12</td></tr><tr><td>40</td><td>3.3e37</td><td>3.7e-42</td><td>1.4e-24</td><td>2.2e-41</td><td>3.5e-23</td><td>1.1e-25</td><td>4.3e-8</td><td>6.6e-25</td><td>1.0e-6</td></tr><tr><td>50</td><td>1.0e46</td><td>7.3e-46</td><td>2.4e-22</td><td>3.6e-45</td><td>1.2e-21</td><td>2.8e-25</td><td>9.8e-2</td><td>1.4e-24</td><td>5.1e-1</td></tr><tr><td>60</td><td>3.1e54</td><td>1.1e-49</td><td>9.6e-21</td><td>5.0e-49</td><td>7.5e-20</td><td>5.8e-25</td><td>5.0e5</td><td>2.6e-24</td><td>3.9e6</td></tr></tbody></table>

[14] Maurice Mignotte. *Mathematics for computer algebra*. Springer-Verlag, New York, 1992.

[15] Ramon E. Moore, R. Baker Kearfott, and Michael J. Cloud. *Introduction to Interval Analysis*. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2009.

[16] Jean-Michel Muller, Nicolas Brunie, Florent de Dinechin, Claude-Pierre Jeannerod, Mioara Joldes, Vincent Lefevre, Guillaume Melquiond, Nathalie Revol, and Serge Torres. *Handbook of floating-point arithmetic*, volume 1. Springer, 2018.

[17] Victor Y. Pan. Univariate polynomials: Nearly optimal algorithms for numerical factorization and root-finding. *Journal of Symbolic Computation*, 33(5):701–733, 2002.

[18] Adrien Poteaux. Computing monodromy groups defined by plane algebraic curves. In *Proceedings of the 2007 International Workshop on Symbolic-numeric Computation*, pages 36–45, New-York, 2007. ACM.

[19] Adrien Poteaux. Calcul de développements de Puiseux et application au calcul de groupe de monodromie d'une courbe algébrique plane. PhD thesis, Université de Limoges, 2008.

[20] Adrien Poteaux and Marc Rybowicz. Good reduction of puiseux series and complexity of the newton-puiseux algorithm over finite fields. In *ISSAC '08: Proceedings of the twenty-first international symposium on Symbolic and algebraic computation*, pages 239–246, New York, NY, USA, 2008. ACM.

[21] Adrien Poteaux and Marc Rybowicz. Good reduction of Puiseux series and applications. *Journal of Symbolic Computation*, 47(1):32–63, 2012.

[22] Adrien Poteaux and Marc Rybowicz. Improving complexity bounds for the computation of Puiseux series over finite fields. In *Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation*, ISSAC '15, pages 299–306, New York, NY, USA, 2015. ACM.

[23] Adrien Poteaux and Martin Weimann. Computing Puiseux series: a fast divide and conquer algorithm. *Annales Henri Lebesgue*, 4:1061–1102, 2021.

[24] Adrien Poteaux and Martin Weimann. Local polynomial factorisation: Improving the Montes algorithm. In *Proceedings of the 2022 ACM on International Symposium on Symbolic and Algebraic Computation*, ISSAC '22, pages 149–158, New York, NY, USA, 2022. ACM.

[25] Jiří Rohn. Enclosing solutions of overdetermined systems of linear interval equations. *Reliable Computing*, 2:167–171, 1996.

[26] Siegfried M. Rump. Verification methods: Rigorous results using floating-point arithmetic. *Acta Numerica*, 19:287–449, 2010.

[27] The FLINT team. *FLINT: Fast Library for Number Theory*, 2025. Version 3.4.0, https://flintlib.org.

[28] J. van der Hoeven. Ball arithmetic. Technical report, HAL, 2009. http://hal.archives-ouvertes.fr/hal-00432152/.

[29] Joachim von zur Gathen and Jürgen Gerhard. *Modern Computer Algebra*. Cambridge University Press, New York, NY, USA, 3rd edition, 2013.

[30] R. J. Walker. *Algebraic Curves*. Springer-Verlag, 1950.

[31] P Walsh. A polynomial-time complexity bound for the computation of the singular part of a Puiseux expansion of an algebraic function. *Mathematics of Computation*, 69(231):1167–1182, 2000.

[32] P. G. Walsh. On the complexity of rational Puiseux expansions. *Pacific J. of Math.*, 188:369–387, 1999.

[33] Zhonggang Zeng. Computing multiple roots of inexact polynomials. *Mathematics of Computation*, 74(250):869–903, 2005.

## A Proof of Proposition 2.2

*PROOF.* Consider the map $\tau: e \in \mathbb{X} \mapsto (d + \mathop{\mathrm{dk}} \cdot e) + \eta \in \mathbb{X}$, so that $e_{t+1} = \tau(e_t)$ in Algorithm KRAWCZYK. Call $d = \text{rad}(d)$, $\eta = \text{rad}(\eta)$, $e_t = \text{rad}(e_t)$ and $e'_t = \text{rad}(e'_t)$.
<table><caption>**Table 2:** Random examples with number of mantissa's bits needed to compute them. Column denoted by *k* is the inverse of the product of differences of roots.</caption><thead><tr><th>prec</th><th>k</th><th>prec</th><th>k</th><th>prec</th><th>k</th><th>prec</th><th>k</th></tr></thead><tbody><tr><td>53</td><td>9.4e-30</td><td>53</td><td>3.6e-8</td><td>70</td><td>5.2e1</td><td>90</td><td>2.5e36</td></tr><tr><td>53</td><td>3.7e-5</td><td>53</td><td>3.3e-29</td><td>120</td><td>6.3e59</td><td>53</td><td>2.7e-21</td></tr><tr><td>53</td><td>7.9e-10</td><td>53</td><td>2.3e-12</td><td>53</td><td>1.1e-25</td><td>53</td><td>0.6</td></tr><tr><td>53</td><td>0</td><td>80</td><td>8.7e25</td><td>70</td><td>1.0e12</td><td>53</td><td>4.0e-22</td></tr><tr><td>53</td><td>2.8e2</td><td>53</td><td>0</td><td>53</td><td>3.1e-27</td><td>120</td><td>5.8e57</td></tr></tbody></table>

For (i), the isotonicity of $\odot d k$ implies that of $\tau$. Moreover, it is clear that $e_1 = d + \eta \subseteq e_2$. Hence the sequence ($e_t$) is increasing (for $\subseteq$), leading to the following alternative:

* Either all components of $e_t$ remain bounded, which make them converge to a limit. Therefore, since all components of $\eta$ are nonzero, we have $e_{t+1} \le e_t + \eta$ for some $t$. Note also that $e_{t+1} \ge e'_{t+1} + \eta$ by line 3. Hence $e'_{t+1} \le e_t$, so that the termination condition in line 6 is fulfilled and the algorithm returns $e_t$.

* Or at least one component of $e_t$ tends to infinity, and the termination-with-error condition in line 4 eventually becomes true.

For (ii), the termination condition in line 6 together with the construction of $e'_{t+1}$ in line 5 implies hypothesis (2) of Theorem 2.1. For (iii), observe that

$$ \begin{aligned} \|\tau(d^+)\|_{\infty} &\le (\|d\|_{\infty} + \frac{1}{4}\|d^+\|_{\infty} + \|\eta\|_{\infty})(1 + O(u'_{\text{prec}})) \\ &\le (\frac{11}{8}\|d\|_{\infty} + O(u'_{\text{prec}}))(1 + O(u'_{\text{prec}})) \\ &< \frac{3}{2}\|d\|_{\infty} = \|d^+\|_{\infty}. \end{aligned} $$

Hence $\tau(d^+) \subseteq d^+$. Since morever $e_1 \subseteq d^+$, we have by isotonicity of $\tau$ that $e_t \subseteq d^+$ for all $t$. Because $\|d^+\|_{\infty} = \frac{3}{2}\|d\|_{\infty} < 2\|e_1\|_{\infty}$, the test in line 4 never leads to an error, and by (i), KRAWCZYK($d$, $\odot d k$) eventually returns $e$ with $\text{rad}(e) = O(\|d^+\|_{\infty}) = O(u_{\text{prec}})$. $\square$

## B Global example

In this section, we present a numerical computation of RPEs of a polynomial $P$ who summarise main parts of VALIDATEDPUISEUX.

The polynomial $P$ is construct from his Puiseux series. These Puiseux series have rational coefficients of two digits at the numerator and the denominator. Some coefficients of the developed form of $P$ have rational coefficients on the order of $10^{57}$ at the numerator and $10^{49}$ at the denominator, for example in $yx^{59}$. The degree of $P$ on $y$ is 14 and the degree on $x$ is 87.

A local good prime for $P$ is $p = 29$. We execute ARNP on the reduction of $P$ modulo $p$, and apply the procedure of the section 3.2 on the resulting tree. We show this tree$^{16}$ in figure 1.

We focus on a few points (labeled (1) to (4) on the tree) that illustrate the different filters, but also the “completion” of the polygon tree of Section 3.2 and why we need to consider collections.

(1) We have two roots of multiplicity 7. Over $\mathbb{Q}$, their value are $r_1 = \frac{59}{52}$ and $r_2 = \frac{16}{15}$, i.e. $r_1 = 24 \mod p$ and $r_2 = 3 \mod p$. Numerically, we get $r_1 = [1.13 \pm 10^{-2}]$ and $r_2 = [1.06 \pm 10^{-2}]$.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Complete polygon tree which can be used to compute ball RPEs of $P$.</figcaption></figure>

As there is no way to connect the numerical approximations to the roots modulo $p$, we will treat simultaneously the two recursive polynomials.

(2) The Newton polygons at this stage are not the one we get from modular computations. Indeed, over $\mathbb{Q}$, we have a coefficient $\frac{29}{4} x^4$ for 4 Puiseux series. Because this coefficient does not lead to any root separation, its cancellation is not detected by our prime criterion. But the procedure of section 3.2 solves the issue, since here $\lambda_{\min} = 1$ is “prolongated” (adding the multiplicity 4). Without that, we would miss this coefficient, and fail later on.

(3) The two roots of multiplicity 2 are treated together. The children of the polygon tree lead to two different Newton polygons, therefore needing POLFILTER. This imply to consider the coefficient in $x^4$ of the two numerical polynomials here.

(4) Similarly, the two roots of multiplicity 4 leads to a collection of two polynomials. Here we have the same Newton polygon, but not the same structure multiplicity. This is therefore MULTFILTER that will enable us to separate the two polynomials and compute the associated roots.

<footnote>$^{16}$This tree is built from Definition 3.2 by applying the Newton polygon and multiplicity structures computations described just below this definition.</footnote>