Joris Van Der Hoeven
email: vdhoeven@lix.polytechnique.fr

Reliable homotopy continuation * Third preliminary version

Keywords: Homotopy continuation, polynomial system solving, ball arithmetic, interval arithmetic, reliable computing 65H20, 65H04, 65G20, 30C15, 13P15

In this paper, we present several algorithms for certified homotopy continuation. One typical application is to compute the roots of a zero dimensional system of polynomial equations. We both present algorithms for the certification of single and multiple roots. We also present several ideas for improving the underlying numerical path tracking algorithms, especially in the case of multiple roots.

Introduction

Besides Gröbner basis computations, homotopy methods are a popular technique for solving systems of polynomial equations. In this paper, we will only consider zero dimensional systems. Given such a system

P (z) = 0, (1)
with P = (P 1 , ..., P n ) and z = (z 1 , ..., z n ), the idea is to find a suitable starting system

Q(z) = 0 (2)
of which all solutions are known, to introduce the homotopy

H(z , t) = (1 -t) P (z) + t Q(z), (3)
and to compute the solutions to P (z) = 0 by following the solutions of Q(z) = 0 from t = 1 to t = 0. Two main approaches exist:

Algebraic homotopies. In this setting, the polynomial equations have exact rational or algebraic coefficients. The homotopy continuation is done exactly using suitable resultants. At the end of the homotopy, the solutions of the system P (z) = 0 are again given exactly, as the solutions of simpler systems. The theory was developed in [GHMP95, GHH+97, Lec01, Dur08] and a concrete implementation is available in the Kronecker system [START_REF] Lecerf | Une alternative aux méthodes de réécriture pour la résolution des systm `es algébriques[END_REF].

Numeric homotopies. An alternative approach is to follow the solution paths using a numeric path tracking algorithm; see [START_REF]Solving polynomial systems using continuation for] engineering and scientific problems[END_REF][START_REF] Verschelde | Homotopy continuation methods for solving polynomial systems[END_REF][START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF] and references therein. This approach is usually faster, partly because most of the operations can be done at a significantly lower precision. However, the end result is only approximate. In particular, it cannot be used for the reliable resolution of overdetermined systems. Several implementations exist for numeric path tracking [START_REF] Verschelde | PHCpack: a general-purpose solver for polynomial systems by homotopy continuation[END_REF][START_REF] Bates | Bertini: software for numerical algebraic geometry[END_REF][START_REF] Leykin | Macaulay 2 package for numerical algebraic geometry[END_REF].

It is surprising that little effort has been undertaken so far in order to bring both approaches closer together. Particularly interesting challenges are how to make numeric homotopy as reliable as possible and how to reconstruct exact end results from the numeric output.

Part of this situation might be due to the fact that interval analysis [Moo66, AH83, Neu90, JKDW01, Kul08, MKC09, Rum10] is not so well-known in the communities where homotopy methods were developed, with the exception of one early paper [START_REF] Kearfott | An interval step control for continuation methods[END_REF]. The main objective paper is to systematically exploit interval analysis techniques in the context of homotopy continuation. We will show how to certify homotopy continuations as well as single and multiple solutions of the polynomial system P (z) = 0. Section 3 is devoted to preliminaries from the area of reliable computation. In section 3.1, we start by recalling the basic principles of ball arithmetic [START_REF] Van Der Hoeven | Ball arithmetic[END_REF], which is a more suitable variant of interval arithmetic for our purposes. In section 3.2, we pursue by recalling the concept of a Taylor model [START_REF] Makino | Remainder differential algebras and their applications[END_REF][START_REF] Makino | Suppression of the wrapping effect by Taylor model-based validated integrators[END_REF], which is useful in order to compute with reliable enclosures of multivariate analytic functions on polydisks. We also introduce a variant of Taylors models in section 3.3, which simultaneously encloses an analytic function and a finite number of its derivatives. In sections 3.4 and 3.5, we discuss the well known problem of overestimation which is inherent to ball arithmetic. We will provide some techniques to analyze, quantify and reduce overestimation.

Before attacking the topic of certified path tracking, it is useful to review the theory of numeric path tracking first. In section 4, we start with the case of non singular paths, in which case we use a classical predictor corrector approach based on Euler-Newton's method. The goal of a numeric path tracker is to advance as fast as possible on the solution path while minimizing the risk of errors. Clearly, the working precision has to be sufficiently large in order to ensure that function evaluations are reasonably accurate. In section 4.4, we show how to find a suitable working precision using ball arithmetic. We consider this approach to be simpler, more robust and more general than the one proposed in [START_REF] Bates | Adaptive multiprecision path tracking[END_REF]. In order to reduce the risk of jumping from one path to another path, we also need a criterion for checking whether our numeric approximations stay reasonably close to the true solution path. A numerically robust way to do this is to ensure that the Jacobian of H does not change to rapidly during each step; see section 4.5 and [START_REF] Bates | Adaptive multiprecision path tracking[END_REF] for a related approach. Another technique is to detect near collisions of paths and undertake special action in this case; see section 4.6.

In section 5, we turn our attention to homotopies (3) such that the end system (1) admits multiple solutions. We will see that Euler-Newton iterations only admit a linear convergence near multiple solutions. Therefore, it is useful to search for alternative itera-tions which admit a better convergence. Now the solution path near a multiple solution is given by a convergent Puiseux series in t. When letting t → e 2pi t turn around the origin, we thus fall on another solution path. The collection of paths which are obtained through repeated rotations of this kind is called a herd. In sections 5.2 and 5.3, we will describe a new path tracking method with quadratic convergence, which operates simultaneously on all paths in a herd. The remaining issue of how to detect clusters and herds will be described in sections 5.4, 5.5 and 5.6.

In section 6, we turn our attention to the certification of single roots of (1) and single steps of a path tracker. An efficient and robust method for the certification of solutions to systems of non linear equations is Krawczyk's method [START_REF] Krawczyk | Newton-algorithmen zur bestimmung von nullstellen mit fehler-schranken[END_REF], with several improvements by Rump [START_REF] Rump | Kleine Fehlerschranken bei Matrixproblemen[END_REF]. In section 6.1, we adapt this classical method to the setting of ball arithmetic. In section 6.2, we will see that an easy generalization of this method provides an algorithm for certified path tracking. An alternative such algorithm was given in [START_REF] Kearfott | An interval step control for continuation methods[END_REF], but the present algorithm presents similar advantages as Krawczyk's method with respect to other methods for the certification of solutions to systems of non linear equations. However, both methods still suffer from overestimation due to the fact that error bounds are computed on a polydisk which contains the solution path. Using the technique of Taylor models, we will show in section 6.3 that it possible to compute the error bounds in small tubes around the actual solution path, thereby reducing the problem of overestimation.

In section 7, we first consider the more difficult problem of certifying multiple roots in the univariate case. We will describe two methods based on Rouché's theorem and a third method which rather aims at certifying a local factorization. The last method also serves as an important ingredient for making the Weierstrass preparation theorem effective in an analytic context.

Before we turn our attention to the certification of multiple roots in the multivariate case, it will be convenient to have a general toolbox for effective complex analysis in several variables. In sections 8.1 and 8.2, we first propose two ways how to formalize "computable analytic functions" We next propose several basic algorithms for computations with such functions.

In section 9, we consider the more difficult problem of certifying multiple roots in the multivariate setting. Several algorithms have been proposed for this problem [OWM83, Lec02, GLSY05, GLSY07, LVZ06, LVZ08, RG10, MM11]. However, most existing strategies require the computation of a large number of derivatives of the system, which becomes prohibitive for large clusters of solutions. In the simplest case when the Jacobian matrix of the polynomial system has corank at most one at the singularity, our solution is based on the simultaneous consideration of all solution paths in a herd. The general case will essentially be reduced to this case using analytic elimination techniques which are similar to the geometric resolution method of polynomial systems [GHMP95, [START_REF] Giusti | Lower bounds for diophantine approximation[END_REF][START_REF] Durvye | Algorithmes pour la décomposition primaire des idéaux polynomiaux de dimension nulle donnés en évaluation[END_REF]. In section 9.8, we will also outline a bridge between numeric and algebraic solutions which provides an alternative way to certify solutions.

In section 10, we study generalizations to the resolution of systems of analytic equations in a given polydisk. In section 10.2, we first consider a system of analytic equations as a perturbation of a system of polynomial equations. If, for each of the solutions to the system of polynomial equations, we can control the speed with which the solution moves under perturbations, then we can solve the perturbed system. Unfortunately, this forces us to keep track of solutions which are far away from the region of interest. In section 10.3, we present an alternative strategy based on incremental resolution, as in [GHMP95, [START_REF] Giusti | Lower bounds for diophantine approximation[END_REF][START_REF] Durvye | Algorithmes pour la décomposition primaire des idéaux polynomiaux de dimension nulle donnés en évaluation[END_REF]. Although this strategy also may require to work outside the region of interest, it does stay closer.

Notations

Positive elements. Given a subset R ⊆ R ∪ {±∞}, we denote

R = {x ∈ R: x 0} R > = {x ∈ R: x = / 0}
Vector notation. Unless stated otherwise, we will use the L 1 -norm for vectors u ∈ C n :

u = |u 1 | + ••• + |u n |. (4)
This norm should not be confused with taking componentwise absolute values

|u| = (|u 1 |, ..., |u n |)
For u, v ∈ R n we also define

u v ⇔ u 1 v 1 ∧ ••• ∧ u n v n u < v ⇔ u 1 < v 1 ∧ ••• ∧ u n < v n max (u, v) = (max (u 1 , v 1 ), ..., max (u n , v n )) u • v = u 1 v 1 + ••• + u n v n
If z 1 , ..., z n are formal variables, then we write

z u = z 1 u 1 ••• z n u n
Matrix notation. We write K r ×c for the set of r × c matrices over a set K. The matrix norm of a matrix M ∈ C r ×c corresponding to the L 1 -norm (4) for vectors

M = sup z =1 M z = i max j |M i,j |.
Directed acyclic graphs. We recall that labeled directed acyclic graphs are often used for the representation of symbolic expressions with potential common subexpressions. For instance, ×

+

x y is a typical dag for the expression (x + y) 2 . We will denote by s f the size of a dag f . For instance, the size of the above dag is s f = 4.

Reliable arithmetic

Ball arithmetic

Let us briefly recall the principles behind ball arithmetic. Given a normed vector space K, we will denote by K or B(K, R) the set of closed balls with centers in K and radii in R . Given such a ball z ∈ B(K, R), we will denote its center by cen(z) and its radius by rad(z). Conversely, given z ∈ K and r ∈ R, we will denote by z + B(r) the closed ball with center z and radius r.

A continuous operation f : K d → K is said to lift into an operation f lift : K d → K on balls, which is usually also denoted by f , if the inclusion property

f (x 1 , ..., x d ) ∈ f (x 1 , ..., x d ) (5)
is satisfied for any x 1 , ..., x d ∈ K and x 1 ∈ x 1 , ..., x d ∈ x d . We also say that f (x 1 , ..., x d ) is an enclosure for the set {f (x 1 , ..., x d ): x 1 ∈ x 1 , ..., x d ∈ x d }, whenever (5) holds. For instance, if K is a Banach algebra, then we may take

x + y = cen(x) + cen(y) + B(rad(x) + rad(y))

xy = cen(x) -cen(y) + B(rad(x) + rad(y))

x y = cen(x) cen(y) + B(rad(x) (|cen(y)| + rad(y)) + |cen(y)| rad(x)).
Similar formulas can be given for division and elementary functions. Certified upper and lower bounds for |x| will be denoted by ⌈x⌉ = |cen(x)| + rad(x) and ⌊x⌋ = max {0, |cen(x)| -rad(x)}.

It is convenient to extend the notion of a ball to more general radius types, which only carry a partial ordering. This allows us for instance to regard a vector x = (x 1 , ..., x d ) ∈ B(K, R) d of balls as a "vectorial ball" with center cen(x) = (cen(x 1 ), ..., cen(x d )) ∈ K d and radius rad(x) = (rad(x 1 ), ..., rad(x d )) ∈ R d . If x = (x 1 , ..., x d ) ∈ K d , then we write x ∈ x if and only if x i ∈ x i for all i ∈ {1, ..., d}. A similar remark holds for matrices and power series with ball coefficients.

In concrete machine computations, numbers are usually approximated by floating point numbers with a finite precision. Let R ˜be the set of floating point numbers at a given working precision, which we will assume fixed. It is customary to include the infinities ±∞ in R ˜as well. The IEEE754 standard [START_REF]ANSI-IEEE Standard[END_REF] specifies how to perform basic arithmetic with floating point numbers in a predictable way, by specifying a rounding mode R ∈ {↓, ↑, } among "down", "up" and "nearest". A multiple precision implementation of this standard is available in the Mpfr library [START_REF] Hanrot | MPFR, a C library for multipleprecision floating-point computations with exact rounding[END_REF]. Given an operation f : R d → R, we will denote by f R : R ˜d → R ˜its approximation using floating pointing arithmetic with rounding mode R. This notation extends to the case when R and R ˜are replaced by their complexifications C and C ˜= R ˜[i].

Let K = R and K ˜= R ˜or K = C and K ˜= C ˜. We will denote by K ˜or B(K ˜, R ˜) the set of closed balls in K with centers in K ˜and radii in R ˜ . In this case, we will also allow for balls with an infinite radius. A continuous operation f : K d → K is again said to lift to an operation f : K ˜d → K ˜on balls if (5) holds for any x 1 , ...,

x d ∈ K ˜and x 1 ∈ x 1 , ..., x d ∈ x d .
The formulas for the ring operations may now be adapted to

x + y = cen(x) + cen(y) + B(rad(x) + ↑ rad(y) + ↑ ǫ +,x,y ) x -y = cen(x) -cen(y) + B(rad(x) + ↑ rad(y) + ↑ ǫ -,x,y ) x y = cen(x) × cen(y) + B(rad(x) × ↑ (|cen(y)| + ↑ rad(y)) + ↑ |cen(y)| × ↑ rad(x) + ↑ ǫ ×,x,y ),
where ǫ +,x,y , ǫ -,x,y and ǫ ×,x,y are reliable bounds for the rounding errors induced by the corresponding floating point operations on the centers; see [START_REF] Van Der Hoeven | Ball arithmetic[END_REF] for more details.

In order to ease the remainder of our exposition, we will avoid technicalities related to rounding problems, and compute with "idealized" balls with centers in K ∈ {R, C} and radii in R . For those who are familiar with rounding errors, it should not be difficult though to adapt our results to more realistic machine computations.

Remark 1. In classical interval analysis so called interval lifts of operations f : K d → K are sometimes required to satisfy the inclusion monotonicity property

x 1 ⊆ y 1 ∧ ••• ∧ x d ⊆ y d =⇒ f (x 1 , ..., x d ) ⊆ f (y 1 , ..., y d ),
for all x 1 , ..., x n , y 1 , ..., y n ∈ K, which clearly implies the usual inclusion property (5). For floating intervals, it is easy to ensure this stronger property using correct rounding. In the ball setting, the exact ring operations in R and C are clearly inclusion monotonic, but it seems cumbersome to preserve this stronger property for floating balls. For this reason, we systematically develop our theory without assuming inclusion monotonicity.

Taylor models

If we are computing with analytic functions on a disk, or multivariate analytic functions on a polydisk, then Taylor models [START_REF] Makino | Remainder differential algebras and their applications[END_REF][START_REF] Makino | Suppression of the wrapping effect by Taylor model-based validated integrators[END_REF] provide a suitable functional analogue for ball arithmetic. We will use a multivariate setup with z = (z 1 , ..., z d ) as our coordinates and a polydisk

D = B(ρ) = {z, |z | |ρ|} for a fixed ρ = (ρ 1 , ..., ρ d ) ∈ (R > ) d .
Taylor models come in different blends, depending on whether we use a global error bound on D or individual bounds for the coefficients of the polynomial approximation. Individual bounds are sharper (especially if we truncate up to an small order such that the remainder is not that small), but more expensive to compute. Our general setup covers all possible blends of Taylor models.

We first need some more definitions and notations. Assume that N d is given the natural partial ordering. Let e k denote the k-th canonical basis vector of N d , so that (e k ) k = 1 and (e k ) l = 0 for l = / k.

For every i ∈ N d , recall that i = |i 1 | + ••• + |i d |. A subset I ⊆ N d
is called an initial segment, if for any i ∈ I and j ∈ N d with j i, we have j ∈ I. In that case, we write I ˚= {i ∈ I: i + {e 1 , ..., e d } ⊆ I } and ∂I = I \ I ˚. In what follows, we assume that I and J are fixed initial segments of N d with J ˚⊆ I. For instance, we may take

I = T n = {i ∈ N d : i n} and J = T n+1 or J = T n or J = {0}. Let K = R or K = C. Given a series f = i∈N d f i z i ∈ K[[z]], we will write supp f = {i ∈ N d : f i = / 0} for its support. Given a subset S ⊆ K[[z]
] and a subset S ⊆ N d , we write f S = i∈S f i z i and S S = {g ∈ S: supp g ⊆ S }. If f is analytic on D, then we denote its sup-norm by

f D = sup z ∈D |f (z)|.
A Taylor model is a tuple P = (ρ, I , J , cen(P ), rad(P )), where ρ, I and J are as above, cen(P ) ∈ K[z] I and rad(P ) ∈ R[z] J . We will write

T = T D,I ,J = B D (K[z] I , R[z] J )
for the set of such Taylor models. Given P ∈ T and i ∈ N d , we will also denote P = cen(P ) + B D (rad(P )) and P i = cen(P ) i + B(rad(P ) i ). Given an analytic function f on D, we write f ∈ P , if there exists a decomposition

f = cen(P ) + i∈J ε i z i with ε i ∈ C[[z]] and ε i D rad(P ) i for all i. In particular, if f ∈ P , then f (z) ∈ i∈I ∪J P i z i ,
for any z ∈ D. Given two Taylor models P , Q ∈ T, we will say that P is included in Q, and we write P ⊆ Q if f ∈ Q for any f ∈ P . This holds in particular if P i ⊆ Q i for all i, in which case we say that P is strongly included in Q and write P ⊑ Q. We finally define ̟(P ) ∈ C by

̟(P ) = P 0 + i= / 0 P i B(ρ) i ,
so that f (z) ∈ ̟(P ) for all f ∈ P and z ∈ B(ρ).

Addition, subtraction and scalar multiplication are defined in a natural way on Taylor models. For multiplication, we need a projection π = π J : N d → J with π(i) i for all i and π(i) = i if i ∈ J . One way to construct such a mapping is as follows. For i ∈ J , we must take π(i) = i. For i ∈ / J , let k be largest such that i k = / 0. Then we recursively define π(i) = π(ie k ). Given P , Q ∈ T, we now define their product by

P Q = i,j ∈I P i Q j B D (ρ) i+j -π(i+j) z π(i+j) .
Using the observation that z i+j ∈ B D (ρ) i+j -π(i+j) z π(i+j) , this product satisfies the inclusion property that f g ∈ P Q for any analytic functions f ∈ P and Q ∈ Q on D.

D-stable Taylor models

For some applications, it is convenient to use Taylor models for enclosing both an analytic function and a certain number of its derivatives. Let us show how to incorporate this in our formalism. Throughout this section, we assume that I = J and that D is an initial segment with D ⊆ I.

Given a Taylor model P ∈ T D,I ,I and i ∈ D, we notice that ∂ i P /∂z i can be regarded as a Taylor model in T D,I ′ ,I ′ with I ′ = {j ∈ N n : i + j ⊆ I }. Let f ∈ D → C be an analytic function and Q ∈ T D,I ,I . We define the relations ∈ D and ⊆ D by

f ∈ D P ⇐⇒ ∀i ∈ D, ∂ i f ∂z i ⊆ ∂ i P ∂z i P ⊆ D Q ⇐⇒ ∀f ∈ D P , f ∈ D Q.
Clearly, P ⊑ Q ⇒ P ⊆ D Q for all P and Q.

Let ω: C d → C be an operation. Then ω is said to D-lift to T D,I ,I , if for all P 1 , ..., P d ∈ T D,I ,I and all

f 1 ∈ D P 1 , ..., f d ∈ D P d , we have ω • (f 1 , ..., f d ) ∈ D ω(P 1 , ..., P d ).
Addition, subtraction and scalar multiplication D-lift in the usual way. As to multiplication, we take

P × D Q = i,j ∈I i+j ∈I P i Q j z i+j + i,j ∈I i+j ∈ / I k∈∂I k i+j C i,j,k P i Q j B D (ρ) i+ j -k z k , with C i,j,k = max l∈D l k ∂ l z i+j z i+j -l ∂z l ∂ l z k z k -l ∂z l .
In order to see that × D satisfies the D-inclusion property, it suffices to check that

z i+j ∈ D z i × D z j for all i, j ∈ I. This is clear if i + j ∈ I. Otherwise, z i × D z j = k∈∂ I k i+j C i,j ,k B D (ρ) i+ j -k z k .
For any l ∈ D with l i + j, there exists a k ∈ ∂I with l k i + j. Hence,

∂ l z i+j ∂z l = ∂ l z i+ j z i+ j -l ∂z l z i+j -l ∈ ∂ l z i+ j z i+ j -l ∂z l B D (ρ) i+j -k z k-l = ∂ l ∂z l ∂ l z i+j z i+j -l ∂z l ∂ l z k z k-l ∂z l B D (ρ) i+ j -k z k ⊆ ∂ l (z i × D z j ) ∂z l .
In the particularly useful case when

I = J = T 1 = {i ∈ N d : i
1}, we notice that C e i ,e j ,e i = C e i ,e j ,e j = 1 for all i = / j and C e i ,e i ,e i = 2 for all i.

Overestimation

The major problem in the area of ball arithmetic is overestimation. For example, even though the expression xx evaluates to zero for any x ∈ R, its evaluation at any ball in R with a non zero radius is not identically equal to zero. For instance,

(1 + B(0.1)) -(1 + B(0.1)) = B(0.2).
Algorithms which rely on ball arithmetic have to be designed with care in order to avoid this kind of overly pessimistic error bounds. In particular, if we evaluate a dag using ball arithmetic, then a symbolically equivalent dag might lead to better error bounds.

Consider a continuous function f : K d → K with K as in section 3.1. We recall that f is said to lift into an operation f lift :

K d → K if the inclusion property f (x) ∈ f lift (x)
is satisfied for all x ∈ K d and x ∈ x. Clearly, such a lift is not unique: for any ε: K d → K with cen ε(x) = 0 for all x, the function f alt = f lift + ε is also a lift of f . If we require that cen f lift (x) = f (cen x), then the best possible lift is given by

f best (x) = f (cen x) + B( sup x ′ ∈x |f (x ′ ) -f (cen x)|).
In general, this lift may be expensive to compute. Nevertheless, its existence suggest the following definition of the quality of a lift. The overestimation χ f lift(x) of f lift at x is defined by

χ f lift(x) = rad f lift (x) rad f best (x) . ( 6
)
This quantity is easier to study if we let rad x tend to zero. Accordingly, we also define the pointwise overestimation function χ

f lift: K d → R by χ f lift(x) = limsup ε→0 χ f lift(x + ε). (7)
Here ε → 0 means that cen ε = 0 and rad ε → 0. If f lift is computed by evaluating a dag f , then it would be nice to have explicit formulas for the pointwise overestimation. For rad x → 0 and assuming that the lift f std is evaluated using the default ball implementations of +, -and × from section 3.1 , we claim that there exists a dag ∇ ¯f with

rad f std (x + B(ε)) = (∇ ¯f ) • |ε| + O(ε 2 ),
for ε → 0. Indeed, we may compute ∇ ¯f using the rules

∇ ¯c = (0, ..., 0) (c ∈ K) ∇ ¯Xk = (0, ... k -1 , 0, 1, 0, ..., 0) (k ∈ {1, ..., r}) ∇ ¯(f ± g) = ∇ ¯f + ∇ ¯g ∇ ¯(f g) = (∇ ¯f ) |g | + |f | (∇ ¯g),
where X k stands for the k-th coordinate function. Now we also have

rad f best (x + B(ε)) = |∇ f | • |ε| + O(ε 2 ), for ε → 0. Consequently, χ f std(x) = limsup ε= / 0 (∇ ¯f )(x) • |ε| |(∇f )(x)| • |ε| .
If d = 1, then this formula simplifies to

χ f std(x) = (∇ ¯f )(x) |f ′ (x)| .
Example 2. With d = 1, let us compare the dags f = X 2 -2 X + 1 and g = (X -1) 2 . We have ∇ ¯f = 2 X + 2 and ∇ ¯g = 2 |X -1|, whence

χ f std(x) = |x| + 1 |x -1| χ g std(x) = 1.
The example shows that we have an infinite amount of overestimation near double zeros, except if the dag is explicitly given as a square near the double zero. More generally, for the dag

f = X n -n X n-1 + n 2 X n-2 + ••• + (-1
) n with an n-fold zero, we obtain

χ f std(x) = (|x| + 1) n-1 |x -1| n-1 .
At a distance ε of the zero, ball arithmetic thus produces bounds which are (2 / ε) n-1 times too pessimistic.

Remark 3. An interesting problem is whether a good understanding of the pointwise overestimation also helps us to bound the overestimation on more general balls. One concrete question is whether we have

rad f std (x) sup x∈x χ f std(x) rad f best (x),
for all polynomial dags f and balls x. This inequality seems to hold in all easy cases that we have looked at, but we do not have a proof that it holds in general.

Reducing the overestimation

The example 2 shows that standard ball arithmetic generally produces an infinite amount of overestimation near double or multiple zeros. This raises the problem how to compute better ball lifts which do not present this drawback. One possible remedy is to systematically compute the ball lifts using Taylor models. Indeed, assume that we want to evaluate f at the ball x = c + B(ρ). Let D = B(ρ), I and J be as in section 3.

(D ¯F c)(ε) = 0 (c ∈ K) (D ¯F X k )(ε) = 0 if (0, ... k-1 , 0, 1, 0, ..., 0) ∈ I ε k otherwise (k ∈ {1, ..., r}) (D ¯F(f ± g))(ε) = (D ¯F f )(ε) + (D ¯F g)(ε) (D ¯F(f g))(ε) = ((D ¯F f ) |g| + |f | (D ¯F g))(ε) + i,j ∈I \{0} i+j ∈F ε i+ j i! j! |f (i) | |g (j) |
For ξ = (x 1 + ǫ 1 , ..., x n + ǫ d ) as above, we then have

f (ξ) ⊆ i∈I 1 i! f (i) (x) ǫ i + B D ((D ¯F f )(ρ)).
Now assume that I = J = T n = {i ∈ N d : i n} and let µ be the valuation of f at x. If µ < n, then we have

rad f best (x + B(ρ)) = sup ε∈B(r) i =µ 1 i! f (i) (x) ρ i + O(ρ µ+1 ) (8) rad f tay (x + B(ρ)) = sup ε∈B(r) i =µ 1 i! f (i) (x) ρ i + O(ρ µ+1 ) (9)
χ f tay(x) = 1. (10)
If µ = n, then we still have (8), but (9) and (10) become rad f tay (x + B(ρ))

(D ¯F f )(ρ) χ f tay(x) = limsup ε= / 0 (D ¯F f )(ε) i =n 1 i! f (i) (x) ε i .
If µ > n, then we generally have

χ f tay(x) = ∞,
although χ f tay(x) < ∞ may occur in lucky cases.

Numeric path tracking

General framework

Let Ω be an open subset of C n and H: Ω × C → C n an analytic function. We consider H as a function H(z, t) in z and the time t, where z ∈ Ω and t ∈ C, and also call H a homotopy.

Assuming that H(z 1 , 1) = 0 for some z 1 ∈ Ω and that we are not in a "degenerate" case, there exists a unique analytic function [0, 1] → Ω: t → z t with H(z t , t) = 0 for all t. We are interested in the value of z t when t → 0. More generally, given a vector

z 1 = (z 1 1 , ..., z 1 k ) ∈ Ω k of vectors, there exists a unique function [0, 1] → Ω k : t → z t with H(z t 1 , t) = ••• = H(z t k , t) = 0 for all t.
The goal of a numeric path tracker is to approximate the function t → z t as well and as quickly possible and, above all, to compute its value z 0 at the "end point" t = 0. In what follows, we will denote by R ˜= R p the set of floating point numbers with p bit mantissas. We also define C ˜= C p = R p [i], Ω ˜= Ω ∩ C ˜and assume that we have a program for computing a numeric approximation H ˜: Ω ˜× C ˜→ C ˜n of H. Given z 1 ∈ Ω ˜with H ˜(z 1 , 1) ≈ 0, we thus want to compute z 0 ∈ Ω ˜with H ˜(z 0 , 0) ≈ 0, by following the homotopy.

In many cases, we will be interested in homotopies for solving a system

P 1 (z) = ••• = P n (z) = 0 (11)
of polynomial equations. The number d of solutions to a generic system of this kind is given by the Bezout number d = d 1 ••• d n , where d i is the total degree of P i for each i. For suitable scaling parameters λ 1 , ..., λ n , σ 1 , ..., σ n ∈ C, we now define H:

C n+1 → C by H 1 (z, t) = (1 -t) P 1 (z, t) + t λ 1 (z d 1 -σ 1 d 1 ) • • • H n (z, t) = (1 -t) P n (z, t) + t λ n (z d n -σ n d n ). Let K = {0, ..., d 1 -1} × ••• × {0, ..., d n -1}
For any k ∈ K, the point

z 1 k = (σ 1 e 2pik 1 /d 1 , ..., σ n e 2pik n /d n )
clearly satisfies H(z 1 k ) = 0, whereas any z 0 k with H(z 0 k , 0) = 0 satisfies P (z 0 k ) = 0. If the system (11) is zero dimensional and the values λ 1 , ..., λ n , σ 1 , ..., σ n are complex and sufficiently random (we also say that the homotopy is in general position), then the system

H 1 (z , t) = ••• = H n (z, t) = 0 is also zero dimensional for every t ∈ [0, 1].
In what follows we will always assume that the homotopy has been chosen in such a way.

Solutions at infinity

One classical difficulty with homotopy methods for solving a polynomial system (11) is that many of the solution paths z t k may tend to infinity in the sense that (z t k ) i → ∞ for some i and t → 0. Computations which infinities can be avoided by rewriting the equations in projective coordinates. More precisely, setting z pr = (z 0 , ..., z n ), the projectivation

A pr ∈ C[z pr ] of a polynomial A ∈ C[z] is defined by A pr (z 0 , ..., z n ) = z 0 deg A A( z 1 z 0 , ..., z n z 0
).

Applying this to the system (11), we obtain a new system

P 1 pr (z pr ) = ••• = P n pr (z pr ) = 0 (12)
of homogeneous equations in z pr . For a random hyperplane

α 0 z 0 + ••• + α n z n = β , (13)
the composite system (12-13) is again zero dimensional, but without solutions at infinity. It is easy to reconstruct solutions to (11) from solutions to (12-13) and vice versa.

Predictor corrector methods

Assume that we have a way to approximate the Jacobian J H of H by J ˜H: Ω ˜→ C ˜n×(n+1) . For instance, if H is given by a dag, then a dag for J H can be computed using forward differentiation, and J ˜H just corresponds to the approximated evaluation of this dag. Assume that we are given y = H ˜(z, t) and J ˜H(z, t) at a certain point where H ˜(z , t) ≈ 0. We may write J ˜H(z, t) = (U , y˙) as the horizontal join of two matrices U ∈ C ˜n×n and y˙∈ C ˜n×1 . Given t ′ = t + d t close to t, we may find a z (1) for which H ˜(z (1) , t ′ ) ≈ 0 using Euler-Newton's method

z (1) = z -U -1 (y + y ˙(t ′ -t)).
The replacement (z, t) (z (1) , t ′ ) is called a prediction step. We may still apply the formula when t ′ = t, in which case z ′ is usually a better approximation than z to a genuine zero of H at t than z. In this situation, the replacement (z, t) (z (1) , t) is called a correction step.

From the computational point of view, the evaluation of the Jacobian J H ˜(z , t) is usually about n times more expensive than the evaluation of the function H ˜(z, t) itself (except for large n and sparse J H ˜). Instead of reevaluating the Jacobian after the prediction step at (z (1) , t ′ ), it may therefore be worth it to perform a few correction steps using the Jacobian at (z, t) instead:

z (2) = z (1) -U -1 H ˜(z (1) , t ′ ) • • • z (κ) = z (κ-1) -U -1 H ˜(z (κ-1) , t ′ ).
Since the convergence of z (1) , z (2) , ... is only linear, the number κ is typically chosen quite small (κ 3). One full prediction-correction cyclus now just consists of the replacement

(z, t) (z ′ , t ′ ) = (z (κ) , t ′ ).
From the complexity point of view, the evaluation of H ˜and J H ˜is usually far more expensive than the cost O(n 3 ) of linear algebra at size n, at least for the examples we will be interested in here. Therefore, it will not be necessary to device the linear algebra algorithms with special care (for instance, we may simply compute the inverse U -1 once and for all, instead of using LU decompositions). On the other hand, we typically want to increase the step size t ′t as much as possible, while trying to stay reasonably close to the true solution path.

Precision control

One obvious source of numeric errors is when the numeric precision being used is insufficient for producing sensible results. In [START_REF] Bates | Adaptive multiprecision path tracking[END_REF], a strategy has been proposed for selecting a sufficient precision for homotopy methods to be numerically reliable. We will now propose an alternative method for finding such a precision, whose justification is based on a simpler argument.

Let p be the current working precision. Our method is based on the following idea: when evaluating y = H ˜(z, t), the actual precision q of the result is usually smaller than p and of the form q = pc for some fixed constant. We will call q the effective precision and we may expect the numeric evaluations to be reliable as long as p is picked sufficiently large such that q τ 0 remains above a certain threshold τ 0 > 0 (e.g. τ 0 = 10).

We still need a more precise definition of the effective precision or a simple way to compute it. Assuming that H ˜admits a ball lift, we may evaluate y = H ˜(z, t) at the ball

(z, t) = (z + B(0), t + B(0)) ∈ B(C p , R p ) n+1 . Then q rel (z, t) = min i log 2 |cen(y i )| rad(y i )
provides an estimate for the relative precision of y. If H ˜(z, t) ≈ 0, then this precision is potentially quite low. In that case, we may also consider q rel (z, t ′ ) at the next time

t ′ = t + d t .
Instead of performing one extra ball evaluation, we may also use the following approximation of q rel (z, t ′ ):

q rel * (z, t ′ ) = min i log 2 |cen(y i ) + H ˜t(z, t) i d t | rad(y i ) .
We now take

q = q(z, t, t ′ ) = min i log 2 max {|cen(y i )|, |cen(y i ) + H ˜t(z , t) i d t |} rad(y i ) .
for the current effective precision at (z, t) and assuming a current step size d t .

Step size control

Since purely numeric homotopy methods are usually being designed for speed, the main focus is not on being 100% fool proof. Nevertheless, it remains worth it to search for cheap ways in order to detect errors and adapt the stepsize so as to avoid potential errors. Now assume that we perform one full prediction correction cyclus (z, t) (z ′ , t ′ ). We first need a criterion for when to accept such a step. The main problem with the design of numeric criteria is there is no way to decide whether a numeric quantity is small or large; such checks can only be performed with respect to other quantities. Instead of checking whether we remain close to the genuine solution path, it is therefore more robust to check that the Jacobian J ˜H does not change not change to quickly on the interval [t, t ′ ].

More precisely, let

y = H ˜(z, t), (U , y˙) = J ˜H(z, t), y ′ = H ˜(z ′ , t ′ ) and (U ′ , y˙′) = J ˜H(z ′ , t ′ ).
Then it is natural to only accept steps for which

U -1 U ′ -1 τ 1 , (14)
for a fixed threshold τ 1 < 1 (e.g.

τ 1 = 1 4
). Here we may use any matrix norm • , so it is most convenient to chose one which is easy to compute:

M = i max j |M i,j |.
The condition ( 14) is not fully satisfactory yet, since it relies on the expensive computation of a Jacobian U ′ . This is acceptable if the step has a good chance of being accepted (since we will need the Jacobian anyway for the next step), but annoying if the step is to be rejected. Before checking (14), it is therefore wise to perform a few cheaper checks in order to increase the probability that (14) will hold indeed. In particular, if κ 2, then we may verify that

U -1 (y ′ -y (1) ) -(z ′ -z (1) ) τ 2 z ′ -z (1) (15)
for the max-norm on vectors, where τ 2 τ 1 (e.g. τ 2 = 1 2 τ 1 ) and y (1) = H ˜(z (1) , t ′ ). This simplified check is linked to (14) by remarking that y ′y (1) ≈ U ′ (z ′z). The new check (15) should not be applied when z ′ and z (1) are too close for y ′ and y (1) to be computed with sufficient precision. More precisely, it should really be replaced by the check

U -1 (y ′ -y (1) ) -(z ′ -z (1) ) τ 2 z ′ -z (1) ∨ z ′ -z (1) 2 -τ 3 q(z,t,t ′ ) z (1) , (16)
where τ 3 is slightly smaller than one (e.g. τ 3 = 3 4 ) and q(z, t, t ′ ) stands for the "effective working precision" from section 4.4.

In addition to the above checks, one might wish to ensure that y ′ is reasonably small after each step. Unfortunately, there is no satisfactory reference with respect which smallness can be checked, except for y (1) , ..., y (κ-1) . The best we can do therefore consists of checking whether y (1) , y (2) , ... tend to 0 at some indicated rate:

y (i+1) τ 4 y (i) ∨ z (i+1) -z (i) 2 -τ 4 q(z,t,t ′ ) z (i) , (17)
for all i < κ, where τ 4 < 1 (e.g. τ 4 = 1 2 ). Again, we need to insert a safety exemption for the case when the convergence is exceptionally good.

Once that we have a criterion on whether a step (z, t) (z ′ , t ′ ) should be accepted, an algorithm for automatic stepsize control is easily implemented: assuming that we are walking from t = 1 to t = 0, we start by setting d t := -1. Given t and d t , we try a step (z, t) (z ′ , t ′ ) until t ′ := t + d t . If the step fails, then we set d t := λ fail d t with λ fail < 1 (e.g. λ fail = 1 2 ), and retry for the smaller stepsize. Otherwise, we accept the step t := t ′ and set d t := λ ok d t for the next step, where λ ok > 1 (e.g. λ ok = 2 √ ).

Near collisions

Another way to look at the numerical error problem is to investigate what can actually go wrong. Theoretically speaking, around each true solution path z t , there exists a small tube T t of variable polyradius r t , where Newton's method converges to the true solution z t .

As long as our current approximation z at time t remains in this tube T t , no errors will occur. Now the Newton iterations have a strong tendency of projecting back into the tubes, especially if we use the additional safeguard (17). Nevertheless, it might happen that we jump from one tube into another tube, whenever two solution paths come close together.

If we are considering a homotopy for solving a polynomial system P 1 = ••• = P n , then various solution paths will actually meet at t = 0 if the system admits multiple roots. Such multiple roots are an intrinsic difficulty and we will need dedicated "end game" strategies to ensure good numeric convergence in this case (see section 5 below).

For t > 0, and for suitably prepared functions H, the Lebesgue probability that two solutions paths meet at a point is zero. Nevertheless, we may have near collisions, which usually occur in pairs: the probability that more than two paths simultaneously pass close to a same point is extremely low.

So assume that we have a near collision of two solution paths. Then we have a true collision at (z * , t * ) for some complex time t * near the real axis. Locally around this collision point, the two paths are then given by

z t ± = z * ± u t -t * √ + O(t -t * ),
for some vector u. If we only know z t + at a few points, then we may try to compute z * , t * and u, and also check whether the second path z t -indeed exists. Now assume that we have approximated z t + and derivative z˙t + = dz t + /dt at two times

t 1 > t 2 . Denote these approximations by z ˜1 = z ˜1 + ≈ z t 1 + , z˙1 ≈ z˙t 1 + , z ˜2 = z ˜2 + ≈ z t 2 + and z˙2 ≈ z˙t 2 + . Then z˙i ≈ u 2 t i -t * √
for i ∈ {1, 2}, whence we may use the following approximations for z * , t * and u:

t ˜ * := (z˙2) 2 t 2 -(z˙1) 2 t 1 (z˙2) 2 -(z˙1) 2 u ˜:= 2 z˙2 t 2 -t ˜ * z ˜ * := z ˜2 -u ˜t2 -t ˜ * .
We next perform several safety checks. First of all, we obtained t ˜ * as the division of two vectors; we may use the mean value of the componentwise divisions and check that the variance remain small. We next verify that z ˜1u ˜t1t ˜ * and z ˜ * are reasonably close.

We also verify that the Newton iteration starting at z ˜2 -= z ˜ *u t 2t ˜ * converges to a solution close to z ˜2 -. We finally verify that the same thing holds for z

˜2 ± = z ˜ * ± u t 2 ¯-t ˜ * instead of z ˜2 ± , where t 2 ¯= Re (2 t ˜ * -t 2 ).
We will not go into technical details on the precise numerical checks here, since section 5.3 below contains a similar discussion for the case of multiple roots at t = 0. We may also adapt the herd iteration from section 5.2 below to near collisions, which allows for the simultaneous continuation of z t + and z t -. Contrary to the case when t → 0, we also need to recompute better estimations of t * at every step, which can be done via the simultaneous computation of z t ± and the two "conjugate" paths z t ± with t ¯= Re (2 t ˜ *t). Indeed, using the higher order expansion

z t ± = z * ± u t -t * √ + v (t -t * ) + w (t -t * ) 3/2 + O((t -t * ) 2 ),
we get

z t + + z t -= 2 z * + 2 v (t -t * ) + O((t -t * ) 2 ) z t + + z t -= 2 z * + 2 v (t ¯-t * ) + O((t ¯-t * ) 2 ) z˙t + + z˙t -= 2 v + O(t -t * ),
from which we may deduce high quality approximations of t * and z * . As soon as t ¯t is small with respect to Im t * , then the junction between paths and their conjugates occurs and we know how to traverse the near collision.

Multiple roots

Straightforward Euler-Newton type methods

Consider a homotopy induced by a polynomial system (11) with a zero dimensional set of solutions. It frequently occurs that some of the solutions are multiple roots, in which case the predictor corrector algorithm slows down significantly when t approaches 0. This is due to the fact that Newton's method only has a linear convergence if we are approaching a multiple root, whereas the convergence is quadratic for single roots.

In order to get a better understanding of this phenomenon, it is instructive to quantify the slow down in the case of an r-fold root of a univariate polynomial P , which is more or less representative for the general case. In the neighbourhood of the root α, we have

P +α (z) := P (α + z) = c z r + O(z r+1 ), with c = 1 r! P (r) (α).
Hence, the Newton iteration becomes

z ′ = z - P +α (z) P +α ′ (z) = 1 - 1 r z + O(z 2 ).
In particular, we see that we need roughly r iterations in order to divide z by e. We also notice that P (α + z) is roughly divided by e at every iteration. For complexity measures, it is more reasonable to study the speed of convergence of P (α + z) rather than z itself. Indeed, the relative precision of an r-fold root is intrinsically r times smaller than the working precision.

If we are rather considering a homotopy H(z, t) = (1t) P (z) + t Q(z), then we usually have q = Q(α) = / 0. Locally, we may thus write

H(α + z, t) = c z r + q t + O(z r+1 ) + O(z t).
Assume that we have H(α + z , t) = 0 for small z and t > 0, so that

z r = - q c t + O(z r+1 ).
Then the Euler-Newton iteration for step size d t yields

z ′ = z - q d t r c z r -1 = 1 - d t r t z + O(z 2 ).
Following our criterion (14), we should have

1 - d t r t r -1 -1 τ 1 .
Roughly speaking, this means that d t τ 1 t. Hence, t is multiplied by 1τ 1 at every step and z is multiplied by 1τ 1 every r steps.

The herd iteration

For high precision computations, it would be nice to have an algorithm with quadratic convergence in t. Before we give such an algorithm, let us first introduce some terminology and study the behaviour of the solutions paths when t → 0.

By assumption, we are given a system (11) with an r-fold root α ∈ Ω. Consider a solution path z t for the homotopy with lim t→0 z t = α. Since z t is algebraic in t, we may expand

z t = α + c 1 t 1/p + c 2 t 2/p + •••,
as a Puiseux series in t for a certain ramification index p (which we assume to be taken minimal). Now letting t turn around 0 once, we have

z e 2p i t = α + c 1 ω t 1/p + c 2 ω 2 t 2/p + •••,
where ω = e 2pi/p . When turning repeatedly, we thus obtain p pairwise distinct solutions paths z t k := z e 2p ik t with k ∈ {0, ..., p -1}. We will call such a family of solution paths a herd. Contrary to the homotopy methods from section 4, which operate on individual paths, the iteration that we will present now simultaneously operates on all paths in a herd. Consider a solution path z t with lim t→0 z t = α as above and the corresponding herd z t k = z e 2p ik t with k ∈ {0, ..., p -1}. We assume that both z ˜t k ≈ z k and z˙t k ≈ z˙t k are known for a given t > 0 and all k ∈ {0, ..., p -1}. Let (F 0 , ..., F p-1 ) and (F ˙0, ..., F ˙p-1 ) denote the FFTtransforms of the vectors (z ˜t 0 , ..., z ˜t p-1 ) and (z˙t 0 , ..., z˙t p-1

) with respect to ω -1 . Then we have

F k = i=0 p-1 z ˜t i ω -ik = p t k/p (c k + c k+p t + O(t 2 )) t F ˙k = t k/p (k c k + (k + p) c k+p t + O(t 2 )).
for all k. We now compute c ˜0, ..., c ˜2p-1 using the formulas

c ˜k+p := 1 p t 1+k/p t F ˙k - k p F k = c k+p + O(t) c ˜k := 1 p t k/p F k -c ˜k+p t = c k + O(t 2 ).
For t ′ > 0 of the order of t 2 , we now have

z t ′ k = z ˜t′ k + O(t 2 ) z ˜t′ k := c ˜0 + c ˜1 ω k (t ′ ) 1/p + ••• + c ˜2p-1 ω (2 p-1)k (t ′ ) (2p-1)/p , (18)
for all k ∈ {0, ..., p -1}. We call (18) the herd prediction. This prediction may be corrected using κ conventional Newton iterations at time t ′ , for a fixed constant κ ∈ N \ {0}.

A complete cyclus of this type will be called a herd iteration.

Step size control for herd iterations

Several technical details need to be settled in order to obtain a robust implementation of herd iterations. First of all, we need a numeric criterion for deciding when the approximations z ˜t k ≈ z t k and z˙t k ≈ z˙t k are of a sufficient quality for starting our herd iteration. Clearly, the error of the approximation should be in O(t 2 ).

We may first ensure ourselves that the approximation can not substantially be improved using Newton iterations: let (z ˜t k ) ′ be the result of applying one Newton iteration to z ˜t k at time t. Then we check whether

relerr(z ˜t k , (z ˜t k ) ′ ) := |z ˜t k -(z ˜t k ) ′ | |z ˜t k | τ 5 t 2 , ( 19
)
for some threshold τ 5 , such as τ 5 = 1 2 (although this check becomes unstable if z ˜t k ≈ 0, we notice that this situation cannot arise systematically for t → 0).

The check (19) for k ∈ {0, ..., p -1} does not yet guarantee that the z ˜t k correspond to approximate evaluations of the Puiseux expansions. In order to check that this is indeed the case, we first compute the c ˜k as described in the previous section. Defining

c ˜(t) = c ˜0 + c ˜1 t 1/p + ••• + c ˜2p-1 t (2 p-1)/p , we next evaluate z ˜t k+1/2 = c ˜(e 2pik+pi t)
for all k ∈ {0, ..., p} and apply one Newton iteration at time t to the results, yielding (z ˜t k+1/2 ) ′ . We now check whether relerr(z

˜t k+1/2 , (z ˜t k+1/2 ) ′ ) τ 6 t 2 , ( 20
)
for some threshold τ 6 , such as τ 6 = 1, and all k. Of course, this second check is more expensive than the first check (19). The thresholds should therefore be adjusted in such a way that the second check is likely to succeed whenever the first one does.

The above criteria can also be used for deciding whether a proposed herd iteration from t to t ′ should be accepted or not. We still have to decide how to chose t ′ . For a fixed constant γ > 1 and a positive integer s which may change at every step, we will take

t ′ = 2 -γ s t.
If a step is accepted, then we increase s by one or a larger integer smaller than 1/ log 2 γ. If a step is not accepted, then we decrease s by one and repeat the same procedure until acceptance or s = 0. If s = 0, then we have either reached the best possible accuracy for the current working precision, or our p paths did not really converge to the same point α. The first case occurs whenever the effective precision from section 4.4 drops below a given threshold. In the latter case, we revert to individual homotopies for further continuation.

Detection of clusters

Let us now go back to the initial polynomial system (11) and assume that we have computed numerical approximations of all

d = d 1 ••• d n individual homotopies (z t k
) k∈K up till a certain time t > 0. We need a way to partition the individual paths into herds. One obvious way is to follow all solution paths from t to e 2pi t and deduce the corresponding permutation of K. However, this computation is quite expensive, so it would be nice to have something faster.

A first step towards the detection of herds is to find all clusters, i.e. all groups of paths which tend to the same limit α. Here we notice that one cluster may contain several herds, as in the example

x 2 = t y 2 = t,
where all four solution paths (x t , y t ) = (ǫ x t √ , ǫ y t √ ) with ǫ x , ǫ y ∈ {-1, 1} tend to the quadruple root (0, 0) of x 2 = y 2 = 0. This cluster contains two herds

(x t , y t ) = (± t √ , ± t √ ) and (x t , y t ) = (± t √ , ∓ t √ ). Now let z ˜t k ≈ z t k and z˙t k ≈ z˙t k for all k ∈ K.
For each k ∈ K, we consider the ball

z t k = z ˜t k + B(2 t z˙t k ).
The radii of these balls has been chosen with care, such that, with high probability, any two paths which belong to the same herd are also in the same connected component of

Z := k∈K z t k
. This is best verified on the case of path

z t = α + c t 1/p + •••. Then the next path in the cluster is z e 2p i t = α + c ω t 1/p + ••• and 1 2 |z e 2p i t -z t | ≈ c 2 |ω -1| t 1/p 2 c p t 1/p ≈ 2 t z˙t.
An efficient way to separate different connected components of Z is via projection. Let λ ∈ R 2n be a random vector of real numbers of length λ = 1. Then any point z ∈ C n may be projected to the vector product π λ (z) = λ • (Re z, Im z) ∈ R. Applying this projection to our balls z t k , we obtain intervals x k . We may sort the x k (and the corresponding z t k ) on their centers in time O(d log d) and compute the various connected components of X := k∈K x k using a linear pass. Whenever x k and x l are in different connected components, then so are z t k and z t l . Assuming that t is sufficiently small, application of this procedure for 2 n random vectors λ results with probability one in the separation of all connected components corresponding to different clusters.

Detection of herds

Let K ′ ⊆ K be a set of indices such that the z k with k ∈ K ′ form a cluster with limit α. We still need a way to find the various herds inside the cluster. In a similar way as in section 5.3, we may improve the quality of our approximations z ˜k and z˙k via Newton iteration until z ˜t k = z t t + O(t 2 ) and z˙t k = z˙t k + O(t). From now on, we assume that we have done this. For each k ∈ K ′ and i ∈ {1, ..., n}, we may write

(z t k ) i = α i + c i k t β i k + •••, for some c i k ∈ C \ {0} and β i k ∈ Q > . We obtain a good approximation A ≈ α + O(t) using α ˜= 1 |K ′ | k∈K ′ z ˜t k . (21)
If |K ′ | is not too large (so that β i k has a small numerator and denominator), then we also obtain reasonably accurate approximations

β ˜i k ≈ β i k and c ˜i k ≈ c i k by β ˜i k = t (z˙t k ) i (z ˜t k ) i -α i c ˜i k = (z ˜k -α) t -β ˜k.
and check whether z e 2p i t k ≈ α ˜+ c ˜k e 2piβ ˜k t β ˜k is indeed close to some z ˜k′ with k ′ ∈ K ′ . Doing this for all k ∈ K ′ , we thus obtain a candidate permutation σ: K ′ → K ′ with z e 2p i t k = z t σ(k) for all k ∈ K ′ . Each cycle in this permutation induces a candidate herd. Using the criteria from 5.3, we may next check whether the quality of the candidate herd is sufficient. If not, then we may always resort to the more expensive computation of the solution path from t to e 2pi t.

Synchronization

Our algorithms for the previous sections for cluster and herd detection rely on the availability of approximations z ˜t k ≈ z t k on all paths at the same time t. Usually the individual homotopies are launched in parallel and advance at different speeds. Consequently, the synchronization of all paths at the same time t is a non trivial matter.

Strictly speaking, we notice that it is not necessary to synchronize all paths, but rather those paths which belong to the same cluster or herd. In particular, we will concentrate on those paths which tend to multiple roots.

So consider a path z t k which tends to a multiple root α. As long as z t k is approximated using an individual continuation, we have seen that the convergence to t → 0 is linear. For a fixed γ < 1 (such as γ = 1 2 ), the computation of z t k at all "checkpoints" t = γ , γ 2 , γ 3 , ... thus only requires a constant overhead. At every checkpoint, we may now launch the algorithm for the detection of clusters. For every candidate cluster K ′ , we next determine the checkpoint γ i with highest i at which z γ i k is available for all k ∈ K ′ . We launch our algorithm for the detection of herds at this checkpoint t = γ i .

In addition, it is a good practice to check that we still have points on all d = d 1 ••• d n paths at every checkpoint. For paths z t k which tend to a single root, we may approximate z γ i k for large i using a single step continuation from t = 0 to t = γ i . For the approximation of α using (21), we notice that it important that no paths of the cluster are missing or counted twice. Indeed, in the contrary case, we only have A = α + O(t β ) with β i = min k∈K ′ β i k for all i, which is insufficient for the computation of accurate approximations of β i k and c i k .

Certified homotopies

Certification of Newton's method

Consider an analytic function f : Ω → C n on some open subset Ω of C n and assume that f admits a ball lift. Given an isolated root z of f , it is well known that Newton's method converges to z in a small neighbourhood of z. It is a natural question to explicitly compute a ball neighbourhood for which this is the case (where we notice that a "ball" in C n is really a compact polydisk). One method which is both efficient and quite tight was proposed by Krawczyk [START_REF] Krawczyk | Newton-algorithmen zur bestimmung von nullstellen mit fehler-schranken[END_REF]. Recall that J f denotes the Jacobian of f .

Theorem 4. Let u ∈ C n , u = cen u and let g: u → C n be an analytic function. Let J g (u) ∈ C n×n be a ball enclosure of the set im J g . If

g(u) -J g (u) B(rad u) ⊆ u,
then g admits a fixed point z ∈ u.

Proof. For any z ∈ u, we have

g(z) = g(u) + 0 1 J g (u + (z -u) t) (z -u) dt.
Since J g (u) is convex, we also have

0 1 J g (u + (z -u) t) dt ∈ J g (u).
Hence

g(z) ∈ g(u) + J g (u) (u -u) ⊆ u.
It follows that g is an analytic function from the compact ball u into itself. By Brouwer's fixed point theorem, we conclude that there exists a z ∈ u with g(z) = z.

Corollary 5. Let u ∈ C n , u = cen u and let V ∈ C n×n be an invertible matrix with

V J f (cen u) ≈ 1. If Ω ⊇ u and u -V f (u) + (1 -V J f (u)) B(rad u) ⊆ u,
then the equation f (z) = 0 admits a root z ∈ u.

Proof. We apply the theorem for g(z

) = z -V f (z).
The above method is still a bit unsatisfactory in the sense that it does not guarantee the uniqueness of the solution. Denoting by int(X) the interior of a subset X of R n , the following sharpening of the method is due to Rump [START_REF] Rump | Kleine Fehlerschranken bei Matrixproblemen[END_REF].

Theorem 6. With the notations from theorem 4, if g(u) -J g (u) B(rad u) ⊆ int(u), then g admits a unique fixed point in u.

Proof. Let us first show that the spectral norm (i.e. the norm of the largest eigenvalue) of any M ∈ J g (u) is <1. Indeed, our assumption implies rad (J g (u) B(rad u)) < rad u.

Now consider the norm

v = max (|v 1 | / rad u 1 , ..., |v n | / rad u n ) on C n . Then, for any M ∈ J g (u) and v with v = 1, we have |M v | rad (M B(|v |)) rad (M B(rad u)) rad (J g (u) B(rad u)) < rad u, whence M v < 1. This is only possible if the spectral norm of M is <1. Now consider ϕ(z) = z -g(z)
. By what precedes, any matrix M in J ϕ (u) = 1 -J g (u) is invertible. For any two distinct points z, z ′ ∈ u, we have

ϕ(z ′ ) -ϕ(z) = 0 1 J ϕ (z + (z ′ -z) t) (z ′ -z) dt. Since J ϕ (u) is convex, there exists a matrix M ∈ J ϕ (u) with M = 0 1 J ϕ (z + (z ′ -z) t) dt.
By what precedes, it follows that ϕ(z ′ )ϕ(z) = M (z ′z) = / 0. We conclude that g(z) = / z or g(z ′ ) = / z ′ . The existence of a fixed point follows from theorem 4.

Corollary 7. With the notations of corollary 5, if

u -V f (u) + (1 -V J f (u)) B(rad u) ⊆ int(u),
then the equation f (z) = 0 admits a unique root z ∈ u.

Proof. Application of theorem 6 for g(z) = z -V f (z).

Assuming that we have computed a numeric approximation z ˜to a root z of f , a second question is how to find a suitable ball z ∋ z ˜for which the corollaries apply. Starting with z 0 := z ˜+ B(0), a simple solution is consider the sequence defined by

z i+1 = cen z i + B(max (rad z i , rad (K(z i ) -cen z i ))) (22)
⊇ z i ∪ K(z i ),
where

K(u) = cen u -V f (cen u) + (1 -V J f (u)) B(rad u)
Whenever K(z i ) ⊆ int(z i ), then we are done. In order to ensure the convergence of this method, we need to tweak the recurrence ( 22) and replace it by

z i+1 = cen z i + B((1 + ε) max (rad z i , rad (K(z i ) -cen z i )) + η), (23)
for suitable small positive constants ε and η. We refer to [START_REF] Rump | Kleine Fehlerschranken bei Matrixproblemen[END_REF] for more details on this technique, which is called ε-inflation.

Certification of a numeric homotopy continuation

Assume that the polynomial system (11) admits only simple roots and that we have obtained numeric approximations z ˜k = z ˜0 k for all these roots using a numeric path tracker. Then theorem 5 suffices for the joint certification of the numeric approximations {z ˜k} k∈K . Indeed, using the above technique, we first compute balls z k ∋ z ˜k for which theorem 5 applies. To conclude, it then suffices to check that these balls are pairwise disjoint. This can be done using the same algorithm as for the detection of clusters, which was described in section 5.4.

In the case when two balls z k and z k ′ do intersect, then we recompute approximations for the paths z t k and z t k ′ using a smaller step size, that is, by lowering the constant τ 1 in (14). We keep doing so until none of the balls z k intersect; even if some of the paths z t k may have been permuted due to numerical errors, the final set of all z k is correct if none of the balls intersect. Indeed, each of the balls contains a solution and there can be no more solutions than the number predicted by the Bezout bound.

If z k and z k ′ intersect then, instead of recomputing the paths z t k and z t k ′ using smaller and smaller step sizes, we may also search for a way to certify the entire homotopy computations. This will be the topic of the remainder of this section. Let us first show how to adapt the theory from the previous section to certified path tracking. From now on, we assume that H: Ω × C → C n is an analytic function which admits a ball lift.

Theorem 8. Let (u, t) ∈ C n × C be such that u ⊆ Ω. Let J = (∂H /∂z)(cen u, cen t) and let V ∈ C n×n be an invertible matrix with V J ≈ 1.

If cen u -V H(cen u, t) + (1 -V ∂ H ∂z (u, t)) B(rad u) ⊆ int(u),
then the equation H(z, t) = 0 admits a unique root z ∈ u for each t ∈ t.

Proof. Let t ∈ t and consider the function g: u → C n ; z → z -H(z, t). Then u -V H(u, t) encloses im g and 1 -V ∂ H ∂z (u, t) encloses im J g , and we conclude by theorem 6. Clearly, for any t, t ′ ∈ t, theorem 8 ensures the existence of a unique solution path from t to t ′ in the tube u × [t, t ′ ]. As at the end of the previous section, the question again arises how to compute balls u and t for which the conditions of the theorem are likely to be satisfied. Since the computation of ∂H ∂z (u, t) is expensive, it is important to keep down the number of iterations of the type ( 22) or (23) as much as possible (say at most one iteration).

Now assume that we performed a numeric homotopy computation from (z, t) to (z ′ , t ′ ). Then a reasonable first guess is to take

u = / 1 2 (z + z ′ ) + B(c (z ′ -z)) t = / 1 2 (t + t ′ ) + B( / 1 2 (t ′ -t)),
for some c > / 1 2 , say c = 1. Unfortunately, if one of the components of z ′z tends to zero, then this guess turns out to be inadequate. Therefore, it is recommended to use an additional inflation proportional to the norm of z ′z:

u = / 1 2 (z + z ′ ) + B(c (z ′ -z) + c ′ z ′ -z ),
for some small c ′ > 0, say c ′ = / 1 10 . Another idea is to use the radius of the previous step as a reference (except for the very first step, of course). For instance, if our previous step went from ( z , t ) to (z, t), then we may take

u = / 1 2 (z + z ′ ) + B(c (z ′ -z) + c ′′ (z -z ) t ′ -t t -t
), for some small c ′′ > 0, say c ′′ = / 1 10 .

Certification via tubular models

One important disadvantage of the method from the previous section for the certification of one path tracking step is that we use global error bounds on the tube u × t. Consequently, the inaccuracy rad u of u is proportional to the step size 2 rad t, whence any overestimation in the evaluation of H or J H due to the inaccuracy in u requires a reduction of the step size.

For this reason, it is much better to follow the solution path as closely as possible instead of enclosing it in a "square tube". This can be achieved via the use of Taylor models. Using D-stable Taylor models, it is possible to simultaneously compute of accurate enclosures for H and J H on the tube.

More precisely, let r ǫ ∈ (R > ) n , r δ ∈ R and D = B(r ǫ ) × B(r δ ). For a fixed k in N \ {0}, let I = J be an initial segment of N n+1 of the form

I = 0 × {0, ..., k} ∪ {e 1 , ..., e n } × {0} and let D = T 1 = {i ∈ N n+1 : i 1}. A D-stable Taylor model in B D (C[ǫ, δ] I , R[ǫ, δ] I )
will also be called a tubular model. We will write T D,I for the set of tubular models. Given y ∈ T D,I n , we let y cst ∈ C n and y lin ∈ C n×n be such that Theorem

(y cst ) i = ̟(y i ) (y lin ) i,j = ̟
9. Let c = (c ǫ , c δ ) ∈ Ω × C, r = (r ǫ , r δ ) ∈ (R > ) n × R > , D = B(r) and let E(δ) = E 0 + ••• + E k 1 δ k be an approximation of the solution ε δ to H(c ǫ + ε δ , c δ + δ) = 0. For instance, if k = 1 and H(c) ≈ 0, then we may take E(δ) ≈ -VF δ, with V ≈ (∂H /∂z)(c) -1 and F ≈ (∂H /∂t)(c). Consider u ∈ T D,I n and v ∈ T D,I with
u i = (c ǫ ) i + ǫ i + E(δ) i v = c δ + δ.
Let g(z, t) = z -V H(z, t), x = g(c ǫ + E(δ), v), y = g(u, v). If

x cst + y lin B(r ǫ ) ⊆ int(c ǫ + B(r ǫ )), (24)
then the equation H(z, t) = 0 admits a unique solution z ∈ c ǫ + E(tc δ ) + B(r ǫ ), for every t ∈ c δ + B(r δ ).

Proof. For an illustration of the proof, see figure 1

. Let u = cen(u) ∈ C[ǫ, δ] n and v = cen(v) ∈ C[δ]
. By construction, and using the facts that ∂u / ∂ǫ = 1 and ∂v / ∂ǫ = ∂E / ∂ǫ = 0, we have

g(u(0, δ), v(δ)) ∈ x cst ∂g ∂ u (u(ǫ, δ), v(δ)
) ∈ y lin for any ǫ ∈ B(r ǫ ) and δ ∈ B(r δ ). For a fixed t ∈ c δ + B(r δ ), it follows that y lin encloses (∂g /∂u)(•, t) on the disk U := c ǫ + E(tc δ ) + B(r ǫ ). Our hypothesis (24) also implies that

g(c ǫ + E(t -c δ ), t) + y lin B(r ǫ ) ⊆ int(U ).
From theorem 6, we conclude that g(•, t) admits a unique fixed point z ∈ U.

In order to apply the theorem, it remains to be shown how to find a good tube, i.e. how to choose c ǫ , c δ , r ǫ , r δ and E(δ). For a fixed order k of the approximation, the idea is to adjust c ǫ and E(δ) such that r ǫ can be chosen minimal.

Let us first consider the first order case k = 1. Assume that we performed a numeric path continuation from (z t , t) to (z t ′, t ′ ) and that both z˙t and z˙t′ are approximatively known. Then there exists a unique curve z ˜s of degree three with z ˜t = z t , z ˜t′ = z t ′, z ˜t = z˙t and z ˜t′ = z˙t′. Let z ˜s be a linear curve which minimizes the maximum µ i of |(z ˜sz ˜s) i | on [t, t ′ ] for every i. Then we take c δ = (t + t ′ )/2, r δ = (t ′t)/2, c ǫ = z ˆct and E(δ) = z ˜ct δ. We may also take r ǫ = c µ for some fixed c > 1 such as c = 2. However, for better performance it is recommended to apply an additional inflation to r ǫ , similar to what we did in the previous section.

For higher orders k, we proceed in an essentially similar way. We first compute a high order numeric polynomial approximation z ˜s of z s . For orders >3, this may require the accurate approximation of additional points (z t ′′, t ′′ ) with t ′′ ∈ (t, t ′ ) on the solution path. We next find a k-th order polynomial z ˆs which approximates z ˜s as good as possible and choose our tube in a similar way as above. It should be noticed that the evaluation g(u, v) in theorem 9 is at least thrice as expensive as the numeric evaluation of J H . This makes it worth it to improve the quality of the numeric approximations of points z t , z t ′, z t ′′ on the curve using one or more additional Newton iterations. The use of higher order approximations makes it possible to choose r ǫ very small, thereby avoiding a great deal of the overestimation due to the use of ball arithmetic.

Numeric spearhead computations and certification

We insist once more on the importance of performing all certifications as late as possible rather than along with the numeric computations themselves. One should regard the numeric computation (the spearhead) as an educated guess of what is happening and the certification as an independent problem at a second stage. In particular, only the numeric results which interests us (i.e. the solutions of the polynomial system) need to be certified and not the way we obtained them (i.e. the homotopies).

Even in the case when we are interested in certifying the homotopies themselves, it is best to do so once the numeric computations have already been completed. This allows for instance for more parallelism in the computations. Indeed, we may cut the entire homotopy path in several pieces and certify these pieces in parallel. More numeric data may also be available once all numeric computations have been completed, which might be useful for guiding and accelerating the certification stage.

Similar remarks apply more generally for certified integration of dynamical systems. In that setting, one is often interested in the flow in the neighbourhood of some initial condition. The sequential part of such a computation resides in the numeric integration for a particular initial condition. Once the corresponding numeric trajectory is known numerically, we may again cut it in pieces which and compute the corresponding flows and certifications in parallel. Whenever the dependence on the initial conditions is very strong, the actual numeric integration should be done using accurate high order Taylor methods using a multiple working precision.

Certification of multiple univariate roots

In section 5.1, we have studied in detail the numeric determination of a multiple root of a univariate polynomial. It is instructive to take up this study and examine how we certify such multiple roots. Since the property of being an r-fold root is lost under small perturbations, this is actually impossible using ball arithmetic. The best we can hope for is to certify the existence of r roots in a small ball, or the existence of an r-fold root of a small perturbation of the polynomial (see also [START_REF] Rump | Verification methods: rigorous results using floating-point arithmetic[END_REF]).

Applying Rouché's theorem

So consider a polynomial P with an approximate r-fold root at c ∈ C and assume that we wish to certify that P admits exactly r roots in the ball c + B(ρ), for some ρ > 0. One first strategy is to make use of the Taylor series expansion of P at c. More precisely, let T = T D,I ,I be the set of univariate Taylor models in ǫ with D = B(ρ) and I = {0, ..., s} for some s r. Evaluating P at c = c + ǫ, we obtain a Taylor model Q = P (c) with the property that P (c + z) ∈ Q 0 + ••• + Q s z s for any z ∈ B(ρ). It remains to be shown that any Q ∈ Q admits r roots in B(ρ). We claim that this is the case if

⌈Q 0 + ••• + Q r-1 B(ρ) r -1 + Q r+1 B(ρ) r+1 + ••• + Q s B(ρ) s ⌉ < ⌊Q r ⌋ ρ n . (25)
Indeed, assume that we have (25) and let

Q ∈ Q. Then |Q(z) -Q r z r | < |Q r z n |
for all z with |z | = ρ. By Rouché's theorem [Lan76, page 158], it follows that Q(z) and Q r z r admit the same number of roots on B(ρ). If r becomes large, or if P admits other roots close to B(ρ), then the bound (25) often does not hold. In that case, one may use more sophisticated techniques from [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF][START_REF] Van Der Hoeven | Efficient root counting for analytic functions on a disk[END_REF] in order to certify that Q admits r roots in B(ρ). From the complexity point of view, the series expansion method requires O(M(r)) evaluations of P , where M(r) denotes the cost of multiplying two polynomials of degrees r.

Computing the winding number

Another approach is to apply Rouché's theorem in a more direct way by computing P on a path γ starting at ρ and which circles around the origin once. If the reliable image P • γ of this path avoids the origin, then the number of roots of P coincides with the number of times that P • γ turns around the origin. More precisely, let ω = e 2pi/R for a suitable R > r (see also below) and let z i = ω i + B(| ω √ -1|) for i ∈ {0, ..., R -1}. Then we evaluate y i = P (z i ) and check whether 0 ∈ / y i for all i. If this is the case, then

r ′ = 1 2 p i=0 R-1 arg cen(y i+1 mod R ) cen(y i )
yields the exact number of roots of P inside B(ρ). This method requires R evaluations of P , but R needs to be sufficiently large if we want to ensure a reasonable chance of success for the method. Let us investigate the choice of an appropriate R in more detail on the simplest example when c = 1 and

P (z) = z r -r z r-1 + r 2 z r-2 + ••• + (-1) r .
Consider the evaluation of P at z = 1 + ρ + B(ǫ). We have

P (z) = P (1 + ρ) + B(σ) = ρ r + B(σ) σ = k=0 r r k ((1 + ρ + ǫ) k -(1 + ρ) k ) = (2 + ρ + ǫ) r -(2 + ρ) r
For small ǫ, the condition 0 ∈ / P (z) thus implies

r (2 + ρ) r -1 ǫ ≈ (2 + ρ + ǫ) r -(2 + ρ) r < ρ r .
Roughly speaking, for ρ → 0, this means that

ǫ < 1 r ρ 2 r -1 ρ R > ρ p ǫ > r 2 ρ r-1 .
We recall from example 2 that (ρ/2) r -1 also corresponds to the punctual overestimation of the ball evaluation of P at 1 + ρ. If we want to reduce R to a quantity which does not depend on ρ, then it follows from the considerations in section 3.5 that we need to evaluate P using Taylor models of order at least r. However, in that case, we might just as well use the first method based on a direct series expansion of P at c.

Certifying a local factorization

Whenever a polynomial P (z) admits r roots α 1 , ..., α r in a disk, then the polynomial

Q = (z -α 1 ) ••• (z -α r ) = z r + Q r -1 z r -1 + ••• + Q 0 is
a monic factor of P . Instead of trying to determine and certify the roots α 1 , ..., α r individually, another idea is to determine and certify this monic factor Q of P .

Assuming that P (z) admits no other roots near {α 1 , ..., α r }, we will show in this section that Q forms an isolated zero of P rem Q = 0, when regarding this equation as a system of r equations

(P rem Q) 0 = ••• = (P rem Q) r -1 = 0 in r variables Q 0 , ..., Q r -1 .
On one hand, this approach has the advantage that we may apply Corollary 7 in order to certify the factorization. On the other hand, we may use fast polynomial arithmetic for the actual evaluation of P rem Q.

Theorem 10. Let α 1 , ..., α ℓ ∈ C be roots of a polynomial P of multiplicities µ 1 , ..., µ ℓ ∈ N > , and

r = µ 1 + ••• + µ ℓ . Then Q = (z -α 1 ) µ 1 ••• (z -α ℓ ) µ ℓ
is an isolated zero of the system P rem Q = 0, when considered as a system of r equations

(P rem Q) 0 = ••• = (P rem Q) r -1 = 0 in r variables Q 0 , ..., Q r -1 .

Proof. Given the euclidean division

P = A Q + B of P by Q, we have A d Q + Q dA + dB = dP = 0 for small perturbations d Q of Q, with deg (dB) < r. Consequently, d(P rem Q) = -((P quo Q) d Q) rem Q.
Computing in the quotient algebra C[x]/(Q), this equation can be reread as 26), the Jacobian matrix of

∂(P rem Q) ∂Q = -(P quo Q) rem Q. ( 26
) Since P = Q R with gcd(R, Q) = 1, the multiplication mapping with R = P quo Q in C[x]/(Q) is invertible. By (
(Q 0 , ..., Q r -1 ) -→ ((P rem Q) 0 , ..., (P rem Q) r -1 )
is precisely the matrix of this multiplication mapping with respect to the canonical basis (1, x, ..., x r-1 ).

Effective complex analysis

For the sequel, it will sometimes be convenient to consider the more general context of analytic functions instead of mere polynomials. We provide two formalizations of computable multivariate analytic functions: an abstract one relying on analytic continuation in section 8.1, as well as a more concrete evaluation based one in section 8.2. In the remainder of the section, we will discuss evaluation of analytic functions in zero dimensional algebras, effective Weierstrass preparation, and computable meromorphic functions.

Computable multivariate analytic functions

We recall [Wei00] that a real number x ∈ R is said to be left (resp. right) computable if there exists an increasing (resp. decreasing) computable sequence (x n ) ∈ Q > with x = lim n→∞ x n . We say that x ∈ R is computable if x is both left and right computable. We denote the sets of computable, left computable and right computable real numbers by R com , R lcom and R rcom . We define C com = R com [i] to be the set of computable complex numbers. The definitions also adapt in a straightforward way to extended real numbers x ∈ R ¯= R ∪ {-∞, ∞}. The theory of computable real numbers provides a suitable abstract framework for studying which analytic problems can be solved.

In [START_REF] Van Der Hoeven | Effective complex analysis[END_REF][START_REF] Van Der Hoeven | On effective analytic continuation[END_REF], we proposed a similar concept of computable analytic functions. Given an analytic function f at the origin, we say that f is computable if there exists methods for computing the power series expansion of f , a lower bound for its convergence radius, an upper bound for f on any closed disk on which f converges, and a method for the analytic continuation of f . Formally speaking, denoting by F com the set of such functions, this means that we may compute

• The computable power series expansion series(f ) ∈ C com [[z]] com of f (this means
that we have an algorithm for the computation of the coefficients of series(f )).

• A lower bound r f ∈ R ¯lcom,> for the radius of convergence r f of f .

• A computable partial function ⌈⌈f ⌉⌉ • : R com,> ⇀ R rcom , which yields an upper bound

⌈⌈f ⌉⌉ ρ f ρ = sup |z | ρ |f (z)| for every ρ < r f . • A computable partial function f +• : C com ⇀ F com , which yields the analytic contin- uation f +δ of f (with f +δ (z) = f (z + δ)) as a function of δ ∈ C com with |δ | < r f .
Given f ∈ F com , we call r f its computable radius of convergence. Usually, r f is smaller than the genuine radius of convergence of series(f ). This definition admits several variants. In practice, it is usually most convenient to provide a method for the computation of bounds ⌈⌈f ⌉⌉ ρ using ball arithmetic and allow for infinite bounds. In that case, we automatically obtain an algorithm for the computation of r f , by taking r f := max {ρ ∈ R com,> : ⌈⌈f ⌉⌉ ρ < ∞}. This definition is also convenient to extend to the case of multivariate analytic functions f in z 1 , ..., z n . In this case, we require algorithms for the computation of:

• The computable power series expansion series(f

) ∈ C com [[z 1 , ..., z n ]] com of f . • A computable partial function ⌈⌈f ⌉⌉ • : (R com,> ) n ⇀ R ¯rcom , which yields a possibly infinite upper bound ⌈⌈f ⌉⌉ ρ f ρ = sup |z | ρ |f (z)|. • A computable partial function f +• : (C com ) n ⇀ F com , which yields the analytic continuation f +δ of f as a function of δ ∈ C com with ⌈⌈f ⌉⌉ |δ | < ∞.
Recall that the function ⌈⌈f ⌉⌉ • is necessarily upper continuous (e.g. [Hoe07, Theorem 2.3]).

In particular, for every ρ with ⌈⌈f ⌉⌉ ρ < ∞ there exists an ε ∈ Q > with ⌈⌈f ⌉⌉ (1+ε)ρ < ∞.

Multivariate analytic functions as evaluable functions

In practice, multivariate analytic functions such as exp (x log(1y)) erf(x 2y 2 ) are often built up as dags from univariate analytic functions such as log, exp and erf. In that case, it would be very expensive to systematically use power series expansions in several variables in order to compute with such functions. Instead, it would be better to represent such multivariate analytic functions as objects which can be evaluated at analytic functions in an arbitrary number of variables, or even at points in more general algebras.

More precisely, let A com ⊇ C com be an effective Banach algebra over C com . This means that A com is the set of computable points in a Banach algebra A over C, and that the operations +, -, × and the norm • : A → R can be computed by algorithm. Recall that

a + b a + b a b a b ,
for all a, b ∈ A. We do not necessarily assume A to be commutative. Given a multivariate analytic function f in z 1 , ..., z n with ⌈⌈f ⌉⌉ ρ < ∞ and pairwise commuting a 1 , ..., a n ∈ A with a i ρ i for i = 1, ..., n, the evaluation

b = f (a 1 , ..., a n ) = k∈N n f k a 1 k 1 ••• a n k n (27) is well defined. What is more: if f ∈ F com and a 1 , ..., a n ∈ A com , then f (a 1 , ..., a n ) ∈ A com . Indeed, we start by computing ε ∈ Q > and M = ⌈⌈f ⌉⌉ (1+ε)ρ < ∞ such that M < ∞. For fixed K ∈ N and b ˜= k 1 ,...,k n <K f k a 1 k 1 ••• a n k n , we then have b ˜-b i=1 n k f k a 1 k 1 ••• a n k n a i K , n M 1 + ε ε n 1 1 + ε Kn .
By choosing K sufficiently large, we may thus make b ˜b as small as desired.

Conversely, assume now that, in the definition of computable multivariate analytic functions, we replace the method series: F com → C com [[z 1 , ..., z n ]] com by an evaluation method eval A com: F com × (A com ) n → A com for any effective Banach algebra A com over C com . In particular, given ρ ∈ (R com ) n , we may take A to be the algebra of all formal power series ϕ ∈ C[[z 1 , ..., z n ]] for which

ϕ = sup k∈N n ϕ k ρ k is finite. Given f ∈ F com with r f > ρ, it
follows that the evaluation eval A com(f , z 1 , ..., z n ) is well defined, and this evaluation yields the power series expansion of f at the origin. This shows that providing an evaluation method eval A com is essentially equivalent to providing a method series for series expansion.

Remark 11. In fact, analytic continuation and bound computation can also be regarded as evaluations in suitable "Banach algebras". Indeed, the analytic continuation at δ ∈ C com corresponds to the evaluation at the analytic function δ + z. The computation of the bound ⌈⌈f ⌉⌉ ρ can be done by evaluating f at the ball B(0, ρ) and taking ⌈⌈f ⌉⌉ ρ = f (B(0, ρ)) , where B(c, r) = |c| + r. Nevertheless, explicit methods for analytic continuation and bound computations are usually of a better quality. They may also be needed for the implementation of more general evaluation methods.

Following the same line of ideas, it may also be useful to consider the evaluation of analytic functions at broken line paths (or more general piecewise analytic paths) instead of ordinary points, thereby combining analytic continuation and ordinary evaluation in a single method. One might even consider evaluations at "paths" A 1 , ..., A k of successive Banach algebras of a similar type, such as

A i = C[z]/ (z 2 -α i z -β i )
with the α i and the β i sufficiently close, and α 1 = β 1 = 0.

Evaluation in commutative zero dimensional algebras

In order to simplify notations, we will now stop our digression on the abstract notions of computability and drop the superscripts "com". In view of what we have seen in section 9.2, it is particularly interesting to evaluate multivariate analytic functions in commutative zero dimensional algebras A over C. Therefore, we will now study this special case in more detail.

Let A be a finite dimensional C-algebra. Given any basis for A, the elements of A can be represented by matrices, so A may be identified with a commutative subalgebra of the algebra C k×k of k × k matrices for some k. In particular, any matrix norm on C k ×k induces a norm on A.

Given a multivariate analytic function f at the origin which is convergent on a polydisk of polyradius ρ ∈ (R > ) n , we may use (27) in order to evaluate f at any points a 1 , ..., a n ∈ A with a i ρ i for all i. Since the a commute, it is actually possible to do a bit better. Indeed, it is classical that there exists an invertible matrix U (corresponding to a base change), such that

U a i U -1 =   T i,1 • • • T i,l   =   λ i,1 + E i,1 • • • λ i,l + E i,l   ,
where λ i,j ∈ C, each E i,j is a nilpotent triagonal k j × k j matrix, and

k = k 1 + ••• + k l . We may thus compute f (a 1 , ..., a n ) using f (a 1 , ..., a n ) = U    f (T 1,1 , ..., T n,1 ) • • • f (T 1,l , ..., T n,l )   U -1 .
For any j with k j = 1, we notice that the evaluation f (T 1,j , ..., T n,j ) reduces to an evaluation f (λ 1,j , ..., λ n,j ) at an ordinary point. If n = 1, then we also notice that

f (T 1,j ) = i=0 r f (i) (λ 1,j ) E 1,j i ,
where r is minimal with E 1,j r+1 = 0.

Weierstrass preparation

Implicit functions

Let us start with a special case of Weierstrass preparation: the implicit function theorem. Consider a computable analytic function

f in z 1 , ..., z n in the neighbourhood of a point α = (α 1 , ..., α n ). Assume that f (α) = (∂f /∂z n )(α) = 0 but (∂f /∂z n )(α) = / 0.
Then the implicit function theorem states that there exists a unique analytic function g = solve(f , α) in z 1 , ..., z n-1 such that f (z 1 , ..., z n-1 , g(z 1 , ..., z n-1 )) = 0 in a neighbourhood of (α 1 , ..., α n-1 ) and g(α 1 , ..., α n-1 ) = α n . It is not hard to check that g is computable, for instance in the sense of section 8.1:

• The coefficients of g can be computed using the classical formulas for implicit functions or using relaxed power series evaluation [START_REF] Van Der Hoeven | Relax, but don't be too lazy[END_REF].

• Bound computations can be done using Corollary 7.

• Given sufficiently small δ 1 , ..., δ n-1 and δ n = g(α 1 + δ 1 , ..., α n-1 + δ n-1 )α n , the analytic continuation of g from (α 1 , ..., α n-1 ) to (α 1 + δ 1 , ..., α n-1 + δ n-1 ) is simply

g +δ = solve(f +δ , α + δ).
More generally, given k analytic functions f = (f 1 , ..., f k ) such that the Jacobian matrix J = ∂ f / ∂z has full rank k at α, then it can be shown in a similar way that the system of k equations in k unknowns g 1 , ..., g k

f (z 1 , ..., z n-k , g 1 (z 1 , ..., z n-k ), ..., g k (z 1 , ..., z n-k )) = 0 g i (α 1 , ..., α n-k ) = α n-k+i
admits a unique analytic solution at α 1 , ..., α n-k which is again computable. Notice that this more general case also follows through a k-fold application of the case of a single function.

Weierstrass preparation

Let us still consider a computable analytic function f in z 1 , ..., z n in the neighbourhood of a point α = (α 1 , ..., α n ). Assume now that

f (α) = (∂ f / ∂ z n )(α) = ••• = (∂ r-1 f /∂z n r-1 )(α) = 0 but (∂ r f / ∂ z n r )(α) = / 0.
Then the Weierstrass preparation theorem states that there exist unique analytic functions g 0 , ..., g r-1 in z 1 , ..., z n-1 and an invertible analytic function u in z 1 , ..., z n such that

f = (z n r + g r-1 z n r -1 + ••• + g 0 ) u. (28)
We may regard (28) as a local factorization of the analytic function f in z n , which depends on z 1 , ..., z n-1 as parameters. Now Theorem 10 for the computation of local factorizations readily generalizes to computable analytic functions. More precisely, g = z n r + g r-1 z n r -1 + ••• + g 0 is the solution of the system f rem z n g = 0, which we consider as a system of r analytic equations

(f rem z n g) 0 = ••• = (f rem z n g) r-1 = 0 in n -1 + r unknowns z 1 , ..., z n-1 , g 0 , .
.., g r -1 . It can be checked that this system of equations admits full rank, so that we may compute g 0 , ..., g r -1 by applying the effective implicit function.

Certifying the multiplicity in a small neighbourhood

Effective Weierstrass preparation as described in the previous section can only be applied on rare occasions. Indeed, the assumption that we have an exact multiple cancellation

f (α) = (∂ f / ∂ z n )(α) = ••• = (∂ r -1 f /∂z n r -1
)(α) = 0 is too ambitious, and needs to be replaced by the weaker condition that f has a multiplicity r in z n in a small neighbourhood of α. This condition is still enough for the existence of a unique solution to the system f rem z n g in a neighbourhood of α.

This naturally leads us to the question how to certify that f has multiplicity r in z n , in a small neighbourhood of α, and uniformly in z 1 , ..., z n-1 . In fact, paying some care, all methods from section 7 can be generalized to this setting. First of all, when using ball arithmetic, the parameters z 1 , ..., z n-1 are simply replaced by small balls α 1 , ..., α n-1 around α 1 , ..., α n-1 , which reduces the problem to a univariate one. Secondly, when working with Taylor models at any chosen truncation order T , the tail

f T (z n -α n ) T + f T +1 (z n -α n ) T +1 + ••• of an analytic function near α n is replaced by a ball of the form B α n (ε) or B α n (ε z n T )
. This reduction allows us to work with polynomials instead of series. We also notice that the above certification methods can actually prove something slightly stronger than a uniform multiplicity: they can usually be used to verify that the zero set Z = {z ∈ α:

f (z) = 0} does not intersect α 1 × ••• × α n-1 × ∂α n .
If such is the case, then we say that the equation f (z) = 0 is equisolvable in z n on α.

Analytic elimination on small balls

Given a ball α ∈ C n (and where we recall that such a "ball" is really a polydisk) such that we can certify a constant multiplicity r of f in z n on α, the unique (and computable) analytic function g = z n r + g r -1 z n r-1 + ••• + g 0 such that f / g is an analytic unit on α will be called the Weierstrass normalization of f on α.

Given two analytic functions f and ϕ on α, assume that we can certify a constant multiplicity for at least one of them, say for f , and let g be as above. Then Weierstrass division of ϕ by g yields unique analytic functions χ and ρ such that ϕ = χ g + ρ and ρ is polynomial in z n with deg z n ρ < r. We may obtain (the residue class of) ρ as a computable analytic function through evaluation of ϕ in the quotient algebra F com / (g) of computable analytic functions in z 1 , ..., z n by the relation g = 0.

Having computed the remainder ρ = ϕ rem g of the Weierstrass division of ϕ by g, we may finally form the resultant R = Res z n (ρ, g) of the polynomials ρ and g in z n . This resultant is a computable analytic function in z 1 , ..., z n-1 whose zeros are precisely the projections of the common zeros of f and ϕ on C n-1 . We will also call R = Res z n (f , ϕ) the analytic resultant of f and g on α. Whenever Res z n (ϕ, f ) is also defined, it can be checked that Res z n (ϕ, f ) = w Res z n (f , ϕ) for some analytic unit w on α 1 × ••• × α n-1 . Indeed, writing f = u g and ϕ = v ψ for analytic units u and v, Res

z n (f , ϕ) ≡ Res z n (g, ϕ) = Res z n (g, ϕ rem g) = Res z n (g, v rem g) Res z n (g, ψ rem g) = Res z n (g, v rem g) Res z n (g, ψ) and Res z n (g, v rem g) is a unit. Similarly, Res z n (ϕ, f ) = Res z n (ψ, u rem ψ) Res z n (ψ, g).

Simplifying computable meromorphic functions

A computable meromorphic function is simply the quotient of two computable analytic functions such that the denominator is non zero. The computable meromorphic functions on a given domain form an effective field (but without a zero test). Using an external argument, we sometimes know that a computable meromorphic function f / g is actually analytic on some domain. In this case, we would like to conclude that f / g is computably analytic on this domain.

A reliable version of l'Hôpital's rule

Let us first consider the case when f and g are both univariate computable analytic functions in z, in the neighbourhood of a single isolated zero α. We start with the subcase when g = zα on a ball α with α ∈ α. Let σ = f ′ (α) be the ball evaluation of f ′ at α. For all z ∈ α, we then we have

f (z) = α z f ′ (u) du ⊆ α z σ du ⊆ σ (z -α).
Consequently, (f / (zα))(α) ⊆ f ′ (α). For general g, we obtain in a similar way that (g /(zα))(α) ⊆ g ′ (α). We conclude that

f g (α) ⊆ f ′ g ′ (α). (29)

Higher multiplicities

It would be interesting to know whether there exist generalizations of (29) to higher multiplicities or involving the successive derivatives of f / g. In this section we will present a more ad hoc strategy for computing a reliable Taylor series expansion of f / g in the case when f , g and f / g are all analytic. Let α and A be all small and a larger ball with the same center c and radii ε and R. Assume that f admits exactly r zeros α 1 , ..., α r both on α and on A.

Let g = (z -α 1 ) ••• (z -α r ) and h = f / g. Writing f i = f (i) (α) / i!, g i = g (i) (α) / i! and h i = h (i) (α)/i! for each i, we have h i ⊆ f i -(h i+1 g r-1 + ••• + h i+r g 0 ). (30)
Let M be an upper bound for |f (A)|. Since h achieves its maximum on A at its border, we also have

|h(A)| M (R -ε) r .
From Cauchy's formula, it follows that

h i ⊆ B(0, M ) (R -ε) i+r (31)
for all i. Plugging this into (30), we obtain

h i ⊆ f i + B(0, M ) (R -ε) i+r 1 (1 -ε/R) r -1 . ( 32
)
For the small balls, we may thus use (32) for the computation of an enclosure for h i , and switch to (31) whenever this enclosure becomes better.

Recovering analytic functions from their values on a circle

More generally, it is in principle possible to recover any analytic function f on a closed disk D from its values on its boundary using Cauchy's formula. We will now describe an efficient method to do this effectively. When applied to computable meromorphic functions as above, this method has the advantage that we do not need to know the numerator and denominator explicitly; it is sufficient that we can compute its ball values on the boundary circle. Without loss of generally, we may assume that D is the unit disk.

For a given number N of evaluation points (preferably a power of two) and with δ = e pi/N -1 , we evaluate f at the disks ω k = B e 2pik/N , δ for k ∈ {0, ..., N -1}, yielding

v k = f (ω k ) = B(v k , ε k ).
We next compute ball enclosures P i for the coefficients of the unique polynomial P with deg P < n and P e 2pik/N = v k for all k using one reliable discrete Fourier transform. The remainder R = f -P is bounded by ε = max k ε k on the unit circle, whence on the unit disk. This yields a Taylor model P + B D (ε) for f . Taking N sufficiently large, the so computed Taylor model can be made as accurate as needed.

Multivariate meromorphic functions

Let us now consider a meromorphic function f / g in several variables z 1 , ..., z n which is actually analytic on some small ball. Modulo a linear change of coordinates and shrinking the domain to a smaller ball α, we may first ensure that g has constant multiplicity in z n on α. We may then use any of the methods from the previous section for computing f / g as an analytic function in z n , which depends analytically on the parameters z 1 , ..., z n-1 .

Certification of multiple multivariate roots

Jacobians of corank one

Consider a system (11) and assume that we have numerically isolated a solution α of multiplicity r. Based on numerical computations, assume also that the Jacobian J = ∂P /∂z is expected to have corank at most one at α. In a small ball α around α (and where we recall that α is really a polydisk) we may certify the corank assumption by computing J(α) using ball arithmetic and then checking that an (n -1) × (n -1) submatrix admits a non zero determinant (still using ball arithmetic).

where A t is monomic of degree r and deg U t,i < r for each i.

Solving the original system Σ for z 1 , ..., z n of this form now amounts to a new system Σ ˜of n r equations in the unknowns (A t ) j and (U t,i ) j . This system can also be considered as a system of n equations in n unknowns z 1 , ..., z n over the algebra A = C[u] / (A t (u)), augmented with one special unknown u r and one special equation z 1 = u. In analogy with Theorem 10, we may hope that this new system is non singular near t = 0, so that we can certify the herd using Corollary 7. It can be checked that this is always so in the corank one case from the previous section. Unfortunately, Theorem 10 fails to generalize in general. In particular, the new system remains singular whenever several herds converge to the same singularity.

Equisolvability

Let us now return to the general case when the corank of the Jacobian matrix at the solution may be larger than one. First of all, it will be useful to generalize the notion of equisolvability from section 8.4.3 to the case of several equations. Consider the system

P k (z) = ••• = P n (z) = 0 (34)
on α ∈ C n , where P k , ..., P n might actually be analytic functions instead of mere polynomials. We say that (34) is equisolvable in z k , ..., z n on α if the zero set

Z = V α (P k ) = {z ∈ α: P k (z) = ••• = P n (z) = 0} (35)
does not intersect α <k × ∂α k , where

α <k = α 1 × ••• × α k -1 and α k = α k × ••• × α n .
In this means that the number of solutions (z k , ..., z n ) ∈ α k of (34) as a function of (z 1 , ..., z k-1 ) ∈ α <k does not depend on z 1 , ..., z k -1 . Furthermore, the following lemma shows that there are no exceptional fibers with an infinite number of solutions.

Lemma 13. Assume that (34) is equisolvable in z k on α and that X ⊆ α <k is the zero set of a non zero analytic function on α <k . Then none of the components of Z can be contained entirely in the set X × α k .

Proof. Assume for contradiction that (X × α k ) ∩ Z has dimension k -1. Taking u ∈ X , the intersection ({u} × α k ) ∩ Z has dimension 1, whence it contains an analytic solution curve t → ϕ(t) = (u, ϕ k (t), ..., ϕ n (t)) ∈ Z whose image has dimension one. Since (34) is equisolvable in z k on α, the curve (ϕ k (t), ..., ϕ n (t)) is entirely contained in α k \ ∂α k . Consequently, the analytic functions ϕ i (t) are all bounded, whence constant. This contradicts our assumption that im ϕ has dimension one.

Example 14. The equation

P (x, y) = 4 y 2 -x = 0
is equisolvable in y on B(0, 1) 2 , since there are always two solutions, but Q(x, y) = y -2 x = 0 is not equisolvable, since the equation has one solution for x = 0, but no solutions for x = 1. Of course, Q is equisolvable on B(0, 1) 2 for y on a larger disk (or x on a smaller disk). Graphically speaking, an equation P (x, y) = 0 is equisolvable if and only if all solution curves y(x) are "nicely horizontal" in B(0, 1) 2 .

Kronecker representations

With the notations from the previous section, consider polynomials

Q = z k d + Q k-1 z k d-1 + ••• + Q 0 V i = V i,d-1 z k d-1 + ••• + V i,0 , i > k,
where the coefficients Q j and V i,j are computable analytic functions in z 1 , ..., z k-1 on α <k , such that Res z k (Q, Q ′ ) = / 0 (with Q ′ = ∂Q/∂z k ) and the system (34) is equivalent to

Q(z k ) = 0 Q ′ (z k ) z i = V i (z k ), i > k, (36)
for all z ∈ α \ (X × α k ), where

X = V α <k (Res z k (Q, Q ′ )).
Then we call (36) the Kronecker representation for (34) on α.

Lemma 15. Assume that the system P 1 (z) = ••• = P n (z) = 0 admits a finite number of solutions in α and that the z 1 -coordinates of these solutions are pairwise distinct. Then there exists a Kronecker representation for P 1 (z) = ••• = P n (z) = 0 on α.

Proof. Let z 1 , ..., z k ∈ C n be such that z 1 1 , ..., z 1 k are pairwise distinct. Then we define the Kronecker representation for Z = {z 1 , ..., z k } to be the unique n-tuple

(Q, V 2 , ..., V n ) = K Z = (Q Z , V 2 Z , ..., V n Z ) of univariate polynomials with deg Q = k, Q k = 1, deg V 2 < k, ..., deg V n < k, such that Z is annihilated by the system Q(z 1 ) = 0 Q ′ (z 1 ) z 2 = V 2 (z 1 ) • • • Q ′ (z 1 ) z n = V n (z 1 ).
We may compute such a Kronecker representation as follows. If k = 1, then we simply take Q(z 1 ) = z 1z 1 1 and V i = z i 1 for all i > 1. For k > 1, we decompose Z = X ∐ Y with |X | = ⌊|Z |/2⌋ and compute K Z in terms of K X and K Y using

Q Z = Q X Q Y V i Z = V i X Q Y + Q X V i Y , i > 1.
This yields an efficient dichotomic algorithm for the computation of K Z . The result follows by applying the method to the solutions of P 1 (z) = ••• = P n (z) = 0 in α.

Theorem 16. Assume that (34) is equisolvable in z k on α. For the zeroset X of some non zero analytic function on α <k , assume that for any u ∈ U := α <k \ X, the solutions of (34) with z <k = u admit pairwise distinct z k -coordinates. Then there exists a Kronecker representation for (34).

Proof. For any u ∈ U , the above lemma yields a Kronecker representation

Q(u, z k ) = 0 Q ′ (u, z k ) z i = V i (u, z k ), i > k
for the system (34) with z <k = u. Moreover, the equisolvability assumption implies that the degree d of Q u does not depend on u. Furthermore, since Q ′ (u, z k ) = / 0 for u ∈ U , the way we compute the functions Q, V k+1 , ..., V n depends analytically on u. We thus obtain a "Kronecker representation" for (34) on the subset U × α k instead of α.

Let M > 0 be such that |α k | M . Expanding Q(u, z k ) = z k d + Q d-1 (u) z k d-1 + ••• + Q 0 (u) and (z k -M ) d = z k d + Q ¯d-1 z k d-1 + ••• + Q ¯0 with Q ¯i = M d-i d i , we have |Q i (u)| Q ¯i for all u ∈ U . Consequently, Q(u, z k ) remains bounded on U × α k .
Since X is a removable isolated singularity, it follows that Q(u, z k ) admits a unique analytic continuation on U. It also follows that we may also choose X to be the zero set of Res z k (Q, Q ′ ).

It remains to be shown that the V i with i > k can also extended analytically from U × α k to α k . For this, we consider a new coordinate z k λ = z k + λ k+1 z k+1 + ••• + λ n z n where λ k+1 , ..., λ n are parameters. For (λ k+1 , ..., λ n ) sufficiently close to (0, ..., 0), we may apply the above result to the system P k (z) = ••• = P n (z) rewritten as equations

P k λ (z λ ) = ••• = P n λ (z λ ) in z λ = z <k , z k λ , z >k .
Notice that this new system is defined on a domain α λ which is slightly smaller than α. We thus obtain a Kronecker representation

Q λ z k λ = 0 (Q λ ) ′ z k λ z i = V i λ z k λ , i > k, which is equivalent to P k λ (z λ ) = ••• = P n λ (z λ ) for z λ with R λ z <k λ = / 0, where R λ = Res z k λ(Q λ , (Q λ ) ′ ).
The function Q λ is defined on α λ and analytic both in λ and z λ . Differentiating the equation Q λ z k λ = 0 with respect to λ i with i > k, we obtain

(Q λ ) ′ z k λ z i = - ∂Q λ ∂λ i z k λ .
Letting λ tend to 0, we obtain that

V i (z k ) = -(∂Q 0 /∂λ i )(z k ). Proposition 17. Assume that P k+1 (z) = ••• = P n (z)

Univariate representations

The univariate representation is a variant of the Kronecker representation. It consists of polynomials

Q = z k d + Q k -1 z k d-1 + ••• + Q 0 U i = U i,d-1 z k d-1 + ••• + U i,0 , i > k,
such that the coefficients Q j are computable analytic functions in z <k on α <k , the coefficients V i,j are computable meromorphic functions in z <k which are analytic on α <k \ X , where

X = V α <k (Res z k (Q, Q ′ )) and Res z k (Q, Q ′ ) = / 0. Finally, the system (34) is assumed to be equivalent to Q(z k ) = 0 z i = U i (z k ), i > k, (37)
for all z ∈ α \ (X × α k ). Given a Kronecker representation, we may compute a univariate representation as follows. We first compute the inverse

I = I d-1 z k d-1 + ••• + I 0 of Q ′ modulo Q,
whose coefficients are meromorphic and analytic on U = α <k \ X . Taking U i = I V i mod Q, we may then retrieve the univariate representation from the Kronecker representation.

Intersection with a new equation

Given a Kronecker representation satisfying the hypothesis of Proposition 17, let us now study the intersection problem with a new equation P k-1 (z) = 0. We first compute a univariate representation for (34), as detailed in section 9.5. Now consider the resultant

R(z <k ) = Res z k (Q(z k ), P k -1 (z k , U >k (z k ))).
In practice, we may compute R by evaluating P k-1 at (z k , U >k (z k )) in the algebra B = A[z k ] / (Q) where A is the set of computable analytic functions on U, and then take the norm of this evaluation. Let φ 1 (z <k ), ..., φ d (z <k ) denote the solutions of (34) in z k as a function of z <k , considered as analytic functions in z <k on some Riemann surface above U . By a classical property of resultants, we have

R(z <k ) = i=1 d P k -1 (z <k , φ i (z <k )) (38)
on U . Whereas the functions φ i (z <k ) generally admit a non trivial monodromy, the value of R(z <k ) is uniquely determined as a function on z <k on U. Moreover, since φ i (z <k ) ∈ α k for z <k ∈ α <k , the analytic function R is bounded on U . Since the complement X of U in α <k is isolated, this means that R extends uniquely an analytic function on α <k .

Lemma 18. Given a Kronecker representation (36) for P k (z) = 0 and assuming that P k (z) = 0 is equisolvable in z k on α, the resultant R vanishes on V α (P k-1 ).

Proof. If z ∈ V α (P k-1 ) ∩ (U × α >k ), then (38) directly yields R(z) = 0. For general z ∈ V α (P k-1 ), let u ∈ C k-1 be a vector such such that z <k + u ε ∈ U for all sufficiently small ε = / 0. Then f i (ε) = P k-1 (z <k + u ε, φ i (z <k + u ε)) is a bounded function, given by a Puiseux series of valuation 0 in ε. Since P k (z) is equisolvable in z k on α, one of the curves (z <k + u ε, φ i (z <k + u ε)) tends to z for ε → 0. Since P k -1 (z) = 0, it follows that f i (ε) and therefore f 1 (ε) ••• f d (ε) tend to 0 for ε → 0. We conclude that R(z) = 0, using the analyticity of R.

We have shown above how to find an analytic relation R for z k-1 which is implied by P k on α. We still need to show how to express z k as a function of z <k . For this we use a similar technique as in the proof of Theorem 16. Consider formal variables ǫ k , ..., ǫ n with ǫ i ǫ j = 0 for all i, j and new coordinates z ǫ such that

z k -1 = z k -1 ǫ -ǫ k • z k ǫ z i = z i ǫ , i = / k -1.
With respect to these new coordinates, the univariate representation for (34) may be rewritten as

Q ǫ (z k ǫ ) = 0 z i ǫ = U i ǫ (z k ǫ ), i > k.
In a similar way as above above, we next evaluate the resultant

R ǫ (z <k ǫ ) = Res z k ǫ (Q ǫ (z k ǫ ), P k-1 (z <k -1 ǫ , z k -1 ǫ -ǫ k z k ǫ -ǫ >k • U >k ǫ (z k ǫ ), z k ǫ , U >k ǫ (z k ǫ ))). = R(z <k ) + ǫ k S k (z <k ) + ••• + ǫ n S n (z <k )
and prove that S k , ..., S n are analytic functions in z <k on α <k . The common zeros of (34) and P k-1 correspond to the zeros of this resultant with respect to the new coordinates.

Setting

z k-1 ǫ = z k-1 + ǫ • z k and R ǫ (z <k ǫ ) = 0 leads to R(z <k ) = 0 R ′ (z <k ) z i = S i (z <k ), i k,
where R ′ = ∂R/∂z k-1 . Assuming that Res z k -1 (R, R ′ ) is not the zero function on α <k -1 , this almost yields a Kronecker representation for P k-1 (z) = 0.

In order to obtain an actual Kronecker representation, we finally compute the Weierstrass normal form Q ˜of R with respect to z k -1 , together with an invertible analytic function W in z <k such that R = W Q ˜. Setting U ˜i = S i /W , we then obtain the Kronecker representation

Q ˜(z <k ) = 0 Q ˜′(z <k ) z i = U ˜i(z <k ), i k for P k-1 (z) = 0 on α, where Q ˜′ = ∂ Q ˜/∂z k -1 . If P k (z) = 0 is equisolvable in z k on α,
then Lemma 18 also implies that Q ˜vanishes on V α (P k-1 ).

Incremental multiplicity certification

We are now in a position to certify the multiplicity of the system (11) on a sufficiently small ball. We first construct a homotopy H for the system (11) which is in sufficiently general position. For notational reasons, it will be convenient to use z 0 instead of t for the time parameter. We take

H i (z 0 , ..., z n ) = (1 -z 0 ) P i (z 1 , ..., z n ) + z 0 I i (z 1 , ..., z n ) I i (z 1 , ..., z n ) = (z i + ε z i+1 + ••• + ε n-i z n ) d i -1, where ε = 1/(2 d 1 ••• d n ).
For k = n, ..., 1 our main objective is to certify that the system H k (z 0 , ..., z n ) = ••• = H n (z 0 , ..., z n ) = 0 is equisolvable in z k , ..., z n on a sufficiently small ball α with 0 ∈ α 0 around a numerical solution. For k = n, we proceed as in section 8.4.3. Assume now that we have proved equisolvability in z >k for the system H >k (z) = 0 on α and that we also have a Kronecker representation

Q ˆ(z k+1 ) = 0 Q ˆ′(z k+1 ) z i = V ˆi(z k+1 ), i > k + 1
for the same system, where Q ˆ′ = ∂ Q ˆ/∂z k+1 . Using the algorithms from section 8.4.3, we may again ensure that Q ˆis equisolvable in z k on α k . Using the algorithm from section 9.6, we next compute a Kronecker representation

Q(z k ) = 0 Q ′ (z k ) z i = V i (z k ), i > k for the intersection V α (H k ), where Q ′ = ∂ Q / ∂ z k . We will show below that Res z k (Q, Q ′ ) = / 0 on α <k . Furthermore, Q vanishes on V α (H k -1
), by Lemma 18. We conclude that the system H k (z) = 0 is equisolvable in z k on α, by Proposition 17. By induction, this completes our algorithm to certify (almost) multiple roots α on a sufficiently small polydisk α around α.

We still have to check that Res z k (Q, Q ′ ) is not identical to zero on α. Assume the contrary, so that Res z k (Q, Q ′ ) vanishes everywhere, by analytic continuation. In other words, for any fixed z <k , the equation Q(z) = 0 admits a multiple solution in z k . In particular, taking z 0 = 1, the generator q(z k ) of the ideal (I k , ..., I n ) ∩ C[z k ] admits a multiple solution. Denoting by z 1 , ..., z s the solutions of the system I k

(z) = ••• = I n (z) = 0, we have q(z k ) = (z k -z k 1 ) ••• (z k -z k s
). Now ε was taken sufficiently small so as to ensure that z k i = / z k j for all i = / j, contradicting the assumption that q admits a multiple root.

Global algebraic certification

Even if the coefficients of the system (11) are all rational or algebraic, then the computed solutions are only numeric approximations. For some applications, it is useful to have exact representations of the solutions. This allows for instance to check whether a given other polynomial with rational coefficients vanishes on the solution set or on some points of the solutions set.

The Kronecker representation provides one useful exact representation for the set of solutions. Modulo a generic linear change of coordinates, we may assume without loss of generality that for any distinct solutions z, z ′ of (11), their first coordinates z 1 and z 1 ′ are also distinct. Using the algorithm in the proof of Lemma 15, we may then compute a Kronecker representation for the system (11).

Assume now that (11) has rational coefficients and that we have computed numeric approximations z ˜1, ..., z ˜k of z 1 , ..., z k with bit precision p. Using the above method, we may thus compute a numeric approximation of the Kronecker representation with an accuracy of approximately p bits. We apply rational number reconstruction [GG02, Chapter 5] in order to provide a guess for the Kronecker representation with rational coefficients. We may check whether this guess is correct by evaluating P at z

1 = u, z 2 = V 2 (u) / Q ′ (u), ..., z n = V n (u) / Q ′ (u) over A = Q[u] / (Q(u)).
If the check fails, then we double the bit precision, use Newton's method to improve the approximations z ˜1, ..., z ˜k, and keep iterating.

In the case that we suspect our system to admit multiple roots, which means that we have at least one herd of numerical solutions, the centers of herds of solutions are still known numerically with a precision of approximately p bits. Consequently, we may use the above method in order to compute the radical of the zero dimensional ideal I generated by the P i . By evaluating the P i over suitable algebras with nilpotent elements, we may check in a similar way as above whether the multiplicity of each root matches with the numeric multiplicity and obtain an exact algebraic representation for I. In the projective setting with no solutions at infinity, this strategy can be used to provide a full algebraically certification of all solutions.

Systems of analytic equations

It is possible to generalize the techniques of this paper to the local resolution of a system of analytic equations on a polydisk. In this section, we outline some ideas for generalizations of this kind. As a general rule, such generalizations usually work better when the local solutions in the polydisk are well separated from the other solutions outside the polydisk.

Univariate equations

Consider the equation

f (z) = 0,
where f is an effective analytic function on the closed unit disk B(0, 1), without any roots on the unit circle. Assume that we wish to find all roots of f in B(0, 1). Using the methods from section 8.4.3, the first step is to compute and certify the number d of roots in B(0, 1).

Having determined the exact number d of roots f inside B(0, 1), most standard numerical methods for finding the d roots of a polynomial of degree d can be mimicked. For instance, we may use the homotopy λ (z dε d ) t + f (z) (1t) = 0 for λ = max |f k | and 0 < |ε| < 1. In the unlucky case that we only find a subset u 1 , ..., u k of the roots with k < d, we set P (z) = (zu 1 ) ••• (zu k ) and keep repeating the algorithm using a homotopy of the form λ (z d-kε ˜d-k ) P (z) t + f (z) (1t) = 0, and for a new random choice of ε ˜with 0 < |ε ˜| < 1. This method will ultimately pick random ε ˜sufficiently close to any of the roots, thereby ensuring the termination of the method. We may also use Aberth's method, while resetting points that fall outside the unit disk to random points inside the disk.

Direct reduction to the polynomial case

Now consider a system

f 1 (z) = ••• = f n (z) = 0 (39)
of equations, where f 1 , ..., f n are computable analytic functions on the closed unit polydisk B(0, 1) n . From now on, our aim is to determine the solutions of this system in B(0, 1) n . For any choice of degrees d 1 , ..., d n , we may consider the truncated polynomials

P i = k d i (f i ) k z k ,
and consider the truncated system

P 1 (z) = ••• = P n (z) = 0. ( 40
)
Using the homotopy methods from this paper, we may compute the solutions of this system and keep only those ones which are in B(0, 1) n . At a second stage, we form the homotopy

H i (z, t) = P i (z) t + f i (z) (1 -t)
and follow the solutions of the truncated systems from t = 1 to t = 0. This yields an uncertified candidate for the set of solutions to (39).

Remark 19. The above method admits several variants. For instance, if we want to avoid the explicit computation of truncated polynomials, then we may directly use the homotopy

H i (z, t) = z i d i -σ i d i t + λ i f i (z) (1 -t)
for suitable σ i ∈ C and λ i ∈ R > with |σ i | = 1. While following the homotopies, we may also decide to drop any paths which lead "too far" outside B(0, 1) n .

Given i, let us now investigate how to pick d i . For any d ∈ N, we define

µ i,d = max k =d |(f i ) k |.
Typically, we now take d i to be the degree d for which µ i,d is maximal (and in any case larger than this value). We may determine this degree by first computing a bound M i for |f i | on a polydisk B(0, ρ) n with ρ > 1, so that |f k | M i ρ -k for all k, whence µ i,d M i ρ -d for all d. Starting with δ := 0 and d := 1, we now repeat the following: if µ i,d > µ i,δ , then δ := d. Else, if M i ρ -d < µ i,δ , then we stop and take d i := δ. Otherwise, set d := d + 1 and continue.

Let us now investigate how to certify that we found all solutions in B(0, 1) n . Since (39) is really a perturbation of (40), one idea would be ensure that for each solution of (40), the solutions remain either inside B(0, 1) n or outside B(0, 1) n for small perturbations. More precisely, for each i, we start by computing a bound B i for |f i -P i | on B(0, 1) n (for instance, we may take B i = k >d i M i ρ -k with the above notations, but is sometimes better to take more explicit coefficients for a sharper bound). We next consider the system P 1 (z) + B(0, B 1 ) = ••• = P n (z) + B(0, B n ) = 0.

For each solution z to (40), we now compute a ball enclosure B(z , η) of this solution such that any solution z ˜to P 1 (z ˜) + ε 1 = ••• = P n (z ˜) + ε n = 0 with |ε i | B i near z lies in B(z, η). This can be done using the ball version of Newton's method. We next check whether each of the obtained enclosures B(z, η) is either entirely contained in B(0, 1) n or in its complement. If this is the case, then we have obtained the desired certification. Indeed, consider a solution z ˜∈ B(0, 1) n to (39). Then setting ε i = f i (z ˜) -P i (z ˜), we have P i (z ˜) + ε i = 0 and |ε i | B i , for all i. Hence z ˜∈ B(z, η) ⊆ B(0, 1) n for one of the above enclosures.

Remark 20. One disadvantage of the above certification method is that we have to keep track of both the solutions inside and outside B(0, 1) n for small perturbations of (40). An alternative idea would be to consider the set S of all solutions to (39) inside B(0, 1) n as a generalized point which is the solution of a suitably deflated system, as in we did in section 9.2 for the certification of herd homotopies. This generalized solution is usually unique in a large polydisk, corresponding to large perturbations of the system of equations for S. With some luck, this polydisk actually contains all sets of |S | points in B(0, 1) n , thereby certifying that S is the set of all solutions to (39) inside B(0, 1) n .

Unfortunately, this idea only works in dimension 1. For instance, in dimension two, the ball of systems X r = {x 2a xb, yc x -d} with a, b, c, d ∈ B(0, r) and r ∈ R > does not contain a system which admits the set {(0, 0), (0, 1)} as its solutions. Nevertheless, it can be checked that any set of two elements in B(0, 1) 2 is the solution of a system in either X 4 or Y 4 , where Y r = {y 2a yb, xc y -d} with a, b, c, d ∈ B(0, r) for any r ∈ R > .

Incremental resolution

As noticed in remark 20, one major disadvantage of the certification method from the previous section is that we need to consider the behaviour of all solutions to (40) under small perturbations, and even of those which lie far outside B(0, 1) n .

The incremental geometric resolution technique from the section 9.7 for the certification of multiple roots can also be used for finding and certifying an arbitrary set of roots in a polydisk α. For this to work, it is necessary that the system (40) is equisolvable on α, which means in practice that we chose a sufficiently general position and that the other roots of (40) outside α are sufficiently far away from α.

Remark 21. In the analytic setting, it can be interesting to compute the numeric solutions to (40) in an incremental way as well. A heuristic way to do the intersection step numerically goes as follows. We assume general position and consider the last intersection (the other intersections are done similarly, in fibers with z 1 = ••• = z k-1 = 0). Suppose that we are given a Kronecker representation Q(0, z 2 ) = 0 Q ′ (0, z 2 ) z i = V i (0, z 2 ), i > 2 for the system f 2 (z) = ••• = f n (z) = 0, in the fiber where z 1 = 0. Let µ be an upper bound for the number of roots of f 1 (z 1 , 0, ..., 0). For each of the roots ω of the equation z 1 µ = 1, we continue our Kronecker representation to a new one in the fiber with z 1 = ω. We finally use the homotopy

(z 1 µ -1) t + f 1 (z) (1 -t) = 0 Q(z 1 , z 2 ) = 0 Q ′ (z 1 , z 2 ) z i = V i (z 1 , z 2 ), i > 2,
and perform the analytic continuation of the µ known solutions at t = 1 until t = 0.

In the case when our method fails to prove equisolvability, we need to decompose the system (39) into simpler systems which are equisolvable. This can be done using the following techniques:

• By covering B(0, 1) n by a finite number of polydisks, each on which the system is equisolvable. For instance in the case of the function Q from example 14, we may cover B(0, 1) 2 by the polydisk B(0, 1) × B(0, 3), or by the polydisk B(0, / 1 2 ) × B(0, / 3 2 ) and several other polydisks of the form B(c, r) × B(0, 1).

• By applying a permutation of the coordinates. For instance the function Q from example 14 is equisolvable in x ∈ B(0, 1) with respect to y ∈ B(0, 1). We may also consider other linear changes of coordinates modulo adjustment of the regions. For instance, the equation x y = 0 cannot be made equisolvable by one of the above means. Nevertheless, after setting x = u + v and y = uv, the equation becomes equisolvable on B(0, 2) 2 .

• By factoring the equation. Especially when an analytic function is restricted to a small region, it frequently occurs that it can be factored into simpler functions on this region. For instance, the equation x 4y 3 + 100 x y = 0 can be factored as 100 In principle, it is always possible to reduce to the equisolvable case using these techniques, but the number of required subdivisions may quickly grow out of hands. Indeed, for each solution (which we assumed to be single), there exists a sufficiently small neighbourhood where the f i are essentially linear functionals. The design of a good algorithm for keeping the number of subdivisions as small as possible is beyond the scope of this paper.

Bibliography

  2 and let T = B D (K[ǫ] I , R[ǫ] J ) be the corresponding domain of Taylor models in ǫ = (ǫ 1 , ..., ǫ d ). Let ξ = (x 1 + ǫ 1 , ..., x n + ǫ d ) ∈ T d and consider the Taylor model evaluation of f at ξ f (ξ) = P + B D (E). Then f tay (x) := P std (B(ρ)) + B(⌈E std (B(ρ))⌉) yields an enclosure of {f (x): x ∈ x}. Although the evaluation of f tay (x) is usually far more expensive than the evaluation of f std (x), let us now study how much the overestimation has been reduced. Let F = N d \ I and let us introduce the operator D ¯F: ε → D ¯F(ε), which generalizes the mapping ε → ε • ∇ ¯. The operator is defined by induction over the size of f :

Figure 1 .

 1 Figure 1. Illustration of a solution path z t in a tube.

  (x + •••) (y + •••) = 0 on B(0, 1) 2 into two equations x + ••• = 0 and y + ••• = 0 which are equisolvable in x resp. y.

Joris van der Hoeven

Reliable homotopy continuation

proof reading by colleagues. * . This work has been supported by the ANR-09-JCJC-0098-01 MaGiX project, the Digiteo 2009-36HD grant and Région Ile-de-France.

At the next stage, we would like to certify that there exist indeed exactly r solutions in a suitable neighbourhood of α, when counted with multiplicities. Now using our effective version of the implicit function theorem from section 8.4.1, we may analytically eliminate at least n -1 variables from the equations (11), say z 2 , ..., z n from the equations P 2 (z) = ••• = P n (z) = 0. This yields computable analytic functions Φ 2 (z 1 ), ..., Φ n (z 1 ) such that P i (z 1 , Φ 2 (z 1 ), ..., Φ n (z 1 )) = 0 for all i ∈ {2, ..., n}. We next compute the computable analytic function R(z 1 ) = P 1 (z 1 , Φ 2 (z 1 ), ..., Φ n (z 1 )).

Using the techniques from section 8.4.3, we finally check that R has multiplicity r in a suitable neighbourhood α 1 of α 1 .

Remark 12. One of the simplest examples when the corank one hypothesis fails is the following system:

Certification of herd homotopies

Let us still consider the system (11) and assume that we are given a herd of numeric solution paths z t 1 , ..., z t r which all tend to the same limit α. In this section, we will investigate the problem of simultaneously certifying each of the paths in the herd.

We would like the follow the same strategy as in section 7.3, where we certified a single factor of order r instead of r separate roots. In the current setting, instead of viewing the z t i as distinct individual paths, this amounts to regarding the whole herd t → Z t = {z t 1 , ..., z t r } as a single multivalued path. From the algebraic point of view, this means that we need to consider the ideal I t which annihilates Z t . There are several ways to represent this ideal I t by a system Σ t of polynomial equations. We will use so called univariate representations. Now recall that each z t i can be considered as a vector z t i = (z t,1 i , ..., z t,n i ) of Puiseux series z t,j i in t of valuations 0. Setting

we notice that A t (z 1 ) is invariant if we turn t once around the origin. This means that A t (z 1 ) is really an analytic function in t at the origin. Assuming general position, A t (z 1 ) is actually the minimal annihilator of {z t,1 1 , ..., z t,1 r }. Furthermore, for t = / 0 sufficiently small, the roots z t,1 i are pairwise distinct, so there are unique interpolating polynomials U t,2 , ..., U t,n ∈ C[u] of degree <r with z t,i = U t,i (z t,1 ) for all i ∈ {2, ..., n}. In what follows, we will represent I t by the system of polynomials