Huu Phuoc Le
email: huu-phuoc.le@lip6.fr

Faster algorithms for computing real isolated points of an algebraic hypersurface

Keywords: Real algebraic geometry, Polynomial system solving, Computational geometry

Let R be the field of real numbers. We consider the problem of computing the real isolated points of a real algebraic set in R n given as the vanishing set of a polynomial system. This problem plays an important role for studying rigidity properties of mechanism in material designs. In this paper, we design two algorithms of whose complexities are respectively O D 8n and O D 6n for computing the real isolated points of real algebraic hypersurfaces in R n of degree D. We also propose several heuristic optimizations to avoid the most costly computation in our algorithms in most of the cases, which makes our complexity reduced to O D 3n for those cases. These algorithms lead to an implementation which is able to solve instances which were out of reach.

Introduction

Given f ∈ Q[x 1 , . . . , x n ], we denote by V ⊂ C n the algebraic hypersurface defined by f = 0 and by H the real trace of V, i.e., H = V ∩ R n . We aim at computing the isolated points of H, i.e. the set of points x ∈ H such that for some positive r, B(x, r) ∩ H = {x} where B(x, r) ⊂ R n the open ball centered at x of radius r. We denote this set of real isolated points by I (H).

This problem is a particular instance of a more general problem of computing the isolated points of a semi-algebraic set. Such problems arise naturally and frequently in the design of rigid mechanism in material design (see [START_REF] Oster | Reentrant tensegrity: A three-periodic, chiral, tensegrity structure that is auxetic[END_REF]). Those are modeled canonically with semi-algebraic constraints, and isolated points to the semi-algebraic set under consideration are related to rigidity properties of the mechanism. A particular example is the study of auxetic materials, i.e., materials that shrink in all directions under compression. These materials appear in nature (first discovered in [START_REF] Lakes | Foam structures with a negative poisson's ratio[END_REF]) e.g., in foams, bones or propylene; see e.g. [START_REF] Yang | Review on auxetic materials[END_REF], and have various potential applications. They are an active field of research, not only on the practical side, e.g., [START_REF] Grima | Auxetic behavior from rotating triangles[END_REF][START_REF] Gaspar | Novel honeycombs with auxetic behaviour[END_REF], but also with respect to mathematical foundations; see e.g. [START_REF] Borcea | Geometric auxetics[END_REF][START_REF] Borcea | Periodic auxetics: structure and design[END_REF]. On the constructive side, these materials are closely related to tensegrity frameworks, e.g., [START_REF] Roth | Tensegrity frameworks[END_REF][START_REF] Connelly | The stability of tensegrity frameworks[END_REF], which can possess various sorts of rigidity properties.

Hence, we aim to provide a practical algorithm for computing these real isolated points in the particular case of real traces of complex hypersurfaces first. This simplification allows us to improve the state-of-the-art complexity for this problem and to establish a new algorithmic framework for such computations which is more efficient in practice.

Prior works.. Let H be a hypersurface defined by f = 0 with f ∈ Q[x 1 , . . . , x n ] of degree D.

A first approach would be to compute a cylindrical algebraic decomposition adapted to H [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition: a synopsis[END_REF]. It partitions H into connected cells, i.e. subsets which are homeomorphic to ]0, 1[ i for some 1 ≤ i ≤ n. Next, one needs to identify cells which correspond to isolated points using adjacency information (see e.g. [START_REF] Arnon | A cluster-based Cylindrical Algebraic Decomposition algorithm[END_REF]). Such a procedure is at least doubly exponential in n and polynomial in D.

A better alternative is to encode real isolated points with quantified formula over the reals. Using e.g. [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]Algorithm 14.21], one can compute isolated points of H in time D O(n 2 ) . Note also that [START_REF] Vorobjov | Complexity of computing the local dimension of a semialgebraic set[END_REF] allows to compute isolated points in time D O(n 3 ) .

A third alternative is to use [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]Algorithm 12.16] to compute sample points in each connected component of H and then decide whether spheres, centered at these points, of infinitesimal radius, meet H. Note that these points are encoded with parametrizations of degree D O(n) (their coordinates are evaluations of polynomials at the roots of a univariate polynomial with infinitesimal coefficients). Applying [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]Alg. 12.16] on this last real root decision problem would lead to a complexity D O(n 2 ) since the input polynomials would have degree D O(n) . Another approach would be to run [2, Alg. 12.16] modulo the algebraic extension used to define the sample points. That would lead to a complexity D O(n) but this research direction requires modifications of [2, Alg. 12.16] since it assumes the input coefficients to lie in an integral domain, which is not satisfied in our case. Besides, we report on practical experiments showing that using [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]Alg. 12.16] to compute only sample points in H does not allow us to solve instances of moderate size.

Since the topological nature of this problem is related to connectivity. Computing isolated points of H is equivalent to computing those connected components of H which are reduced to a single point (see [17, Lemma 1]). Hence, one can consider computing roadmaps: these are algebraic curves contained in H which have a non-empty and connected intersection with all connected components of the real set under study. The authors of [START_REF] El Din | A nearly optimal algorithm for deciding connectivity queries in smooth and bounded real algebraic sets[END_REF] designed a probabilistic algorithm that computes roadmaps for smooth bounded real algebraic sets. This algorithm runs in time O((nD) 12n log 2 n ), which makes it the state-of-the-art complexity for roadmap computation. These results makes plausible to obtain in [START_REF] Le | Computing the real isolated points of an algebraic hypersurface[END_REF] an algorithm running in time (nD) O(n log n) to compute the isolated points of H.

Recently, [START_REF] Basu | Quantitative curve selection lemma[END_REF] presents a new algorithm for computing local semi-algebraic paths of a complexity D O(n) . Their algorithm also allows to solve the algorithmic algorithm we target.

Main results.. We present two probabilistic algorithms which take as input a polynomial f ∈ Q[x 1 , . . . , x n ] and compute the set of isolated points I (H) of H = V( f ) ∩ R n . Our output consists of a zero-dimensional parametrization C = (w(t), v 1 (t), . . . , v n (t)) ∈ Q[t] n+1 such that {(v 1 (t), . . . , v n (t)) | w(t) = 0} is a finite set containing I (H) and a set B = (I 1 , . . . , I s ) of intervals in R that satisfies:

• Each I i has rational endpoints and contains exactly one real root of w(t), namely t i .

• The set of isolated points of H is exactly

{(v 1 (t i ), . . . , v n (t i )) | 1 ≤ i ≤ s}.
These data represent symbolically the set of isolated points of H in the sense that one can derive from the pair (C , B) numerical values of the elements of I (H) with any required precision.

We sketch now the geometric ingredients which allow us to compute the real isolated points of an algebraic hypersurface defined by f = 0.

Assume that f is non-negative over R n (if this is not the case, we can replace it by its square) and let x ∈ I (H). Since x is isolated and f is non-negative over R n , the intuition is that for a small enough e > 0, the real solution set to f = e looks like a ball around x, hence a bounded and closed connected component C x . Then the restriction of every projection on the x i -axis to the algebraic set in C n defined by f = e intersects C x . When e tends to 0, these critical points in C x "tend to x". This first process allows us to compute a subset of H containing I (H); we call the elements of this set the candidates. Of course, one would like that this set of candidates is finite and this will be the case up to some generic linear change of coordinates, using e.g. [START_REF] El Din | Polar varieties and computation of one point in each connected component of a smooth real algebraic set[END_REF].

At this stage, we dispose of a zero-dimensional parametrization C = (w(t), v 1 (t), . . . , v n (t)) encoding the candidates. One would naturally check whether a ball of infinitesimal radius centered at each candidate intersects H. More precisely, let η = (η 1 , . . . , η n ) ∈ C∩R n be a candidate, we decide whether the system

f (x 1 , . . . , x n ) = n i=1 (x i -η i ) 2 -ε = 0 (1)
has at least one solution in R n . However, this direct approach faces immediately a complexity issue. Writing a quantified formula to solve this decision problem leads to a complexity of D O(n 2 ) since those points are encoded by a zero-dimensional parametrization of degree D O(n) . To bypass this difficulty, we carry out the computation over the quotient ring Q[t]/⟨w(t)⟩.

This computation relies on the results of [25, Appendix J], which provide an adaptation of the geometric resolution [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF] to polynomial systems with coefficients in Q[t]/⟨w(t)⟩. Using this version of geometric resolution, the resolution of polynomial systems over Q[t]/⟨w(t)⟩ induces only an additional cost of O (deg(w)) arithmetic operations of Q comparing to the classic geometric resolution. We will see that the degree of w(t) is actually bounded by 2D(D -1) n-1 . This allows us to obtain an algorithm that uses D O(n) arithmetic operations in Q.

On the algorithmic side, we go further exploiting the geometry of the problem to avoid using infinitesimals. We apply the algorithm to compute a rational number e 0 > 0 that replaces the infinitesimals in the system above.

These ingredients allow us to obtain the complexity result below.

Theorem I. Let f ∈ Q[x 1 , . . . ,
x n ] of degree D and H ⊂ R n be the real algebraic set defined by f = 0. There exists a probabilistic algorithm which, on input f , computes a data (C , B) encoding I (H) in case of success using O 64 n D 8n arithmetic operations in Q and one call of real root isolation on a univariate polynomial of degree bounded by 2 2n+2 D 3n .

Furthermore, we propose an alternative variant that leads to a more efficient algorithm in practice. Once the rational number e 0 is computed, this variant replaces the candidates by their approximations of coordinates in Q and solves a similar decision problem as the one given by the system (1) for these approximations. Such a strategy allows us to avoid the computation over Q[t]/ ⟨w(t)⟩. In Subsection 4.4, we will define rigorously this notion of approximations.

Since this variant makes use of univariate real root isolating algorithms, a complete complexity analysis would require a bound on the bit-size of polynomials given to real root isolating 3 algorithms. However, we observe that in practice these real root isolating steps are negligible compared to the computation over Q[t]/ ⟨w(t)⟩, this variant is therefore much more efficient in practice. A complexity estimate of this variant is given as below.

Theorem II. Let f ∈ Q[x 1 , . . . , x n ] of degree D and H ⊂ R n be the real algebraic set defined by f = 0. There exists a probabilistic algorithm which, on input f , computes a data (C , B) encoding I (H) in case of success using O D 6n+3 arithmetic operations in Q and two real root isolating calls on a univariate polynomial of degree bounded by 2D(D -1) n-1 .

Our complexity results above both lie in D O(n) as the algorithm in [START_REF] Basu | Quantitative curve selection lemma[END_REF] (up to some root isolation calls). However, we provide an explicit constant in the exponent and the algorithms are designed to enable practical implementation, which constitutes another contribution of this paper.

In Section 6, we present two heuristic subroutines for our implementation. The goal of these subroutines is to avoid as much as possible the most costly computations in our algorithm which include the computation over quotient rings. Taking into account these optimizations, we implement our algorithms in Maple using the libraries FGb [START_REF] Faugère | FGb: A Library for Computing Gröbner Bases[END_REF], RAGlib [START_REF] El Din | Real alebraic geometry library[END_REF] and msolve [START_REF] Berthomieu | Msolve: A library for solving polynomial systems[END_REF]. In Section 7, we report on practical experiments showing that they already allow us to solve nontrivial problems which are actually out of reach of [2, Alg. 12.16] which computes sample points in H only. Unfortunately, the real-life applications coming from material designs still remain intractable.

Organization of the paper.. This paper is structured as follows. Section 2 reviews some fundamental notions in real algebraic geometry that we use further for the geometry analysis of our problem. In Section 3, we prove the auxiliary results which coin the needed theoretical ingredient. Section 4 is devoted to describe our algorithms and their complexities. The heuristic optimizations for implementations are presented in Section 6. Finally, Section 7 reports on the practical performances of our implementation.

Acknowledgment.. The author would like to thank Mohab Safey El Din for helpful discussions and Timo de Wolff for introducing this problem to me.

Preliminaries

In what follows, we recall respectively the notions of Puiseux series and critical points. Puiseux series.. Let ε be an infinitesimal, i.e., a transcendental element over R such that 0 < ε < r for any positive element r ∈ R. The field of Puiseux series over R is defined as all the series with rational exponents

R⟨ε⟩ = i≥i 0 a i ε i/q | i ∈ N, i 0 ∈ Z, q ∈ N -{0}, a i ∈ R .
By, e.g., [2, Theorem 2.91], R⟨ε⟩ is a real closed field.

One can also define C⟨ε⟩ as for R⟨ε⟩ but taking the coefficients of the series in C. By [2, Theorem 2.17], the field C⟨ε⟩ is an algebraic closure of R⟨ε⟩.

Consider σ = i≥i 0 a i ε i/q ∈ R⟨ε⟩ with a i 0 0. Then, a i 0 is called the valuation of σ. When i 0 ≥ 0, σ is said to be bounded over R and the set of bounded elements of R⟨ε⟩ is denoted by R⟨ε⟩ b . One defines the function lim ε : R⟨ε⟩ b → R that maps σ to a 0 (which is 0 when i 0 > 0) and writes lim ε σ = a 0 ; note that lim ε is a ring homomorphism from R⟨ε⟩ b to R. All these definitions extend to R⟨ε⟩ n componentwise. For a semi-algebraic set S ⊂ R⟨ε⟩ n , we naturally define the limit of S as lim ε S = {lim ε x | x ∈ S and x is bounded over R}.

We refer to [2, Chap. 2] for more details on infinitesimals and real Puiseux series.

Critical points.. Let S ⊂ R n be a semi-algebraic set defined by a semi-algebraic formula Φ. We denote by ext(S, R⟨ε⟩) the semi-algebraic set of points which are solutions of Φ in R⟨ε⟩ n . Let K be an algebraically closed field, let ϕ ∈ K[x 1 , . . . , x n ] which defines the polynomial mapping (x 1 , . . . , x n ) → ϕ(x 1 , . . . , x n ) and V ⊂ K n be a smooth equidimensional algebraic set. We denote by crit(ϕ, V) the set of critical points of the restriction of ϕ to V. If c is the codimension of V and (g 1 , . . . , g s ) generates the vanishing ideal associated to V, then crit(ϕ, V) is the subset of V at which the Jacobian matrix associated to (g 1 , . . . , g s , ϕ) has rank less than or equal to c (see e.g., [START_REF] El Din | A nearly optimal algorithm for deciding connectivity queries in smooth and bounded real algebraic sets[END_REF]Subsection 3.1]).

Geometric results

The candidates

As above, let f ∈ Q[x 1 , . . . , x n ], V ⊂ C n be the algebraic hypersurface defined by f = 0 and H = V ∩ R n . We recall that I (H) denotes the set of isolated points of H.

By, e.g., [START_REF] Le | Computing the real isolated points of an algebraic hypersurface[END_REF]Lemma 1], the set I (H) is the finite union of the connected components of H which are singletons. A natural start for computing those components is to consider the critical points of the restrictions to H of the projections

π i : (x 1 , . . . , x n ) → x i .
Since we do not assume that V is smooth, to apply the critical point method, we retrieve the smoothness through deformation techniques by introducing an infinitesimal ε. More precisely, we consider the hypersurface V ε in C⟨ε⟩ n defined by the equation f 2 = ε 2 . Note that V ε is the union of two disjoint algebraic sets V( f -ε) and V( f + ε) in C⟨ε⟩ n . By e.g., [START_REF] Rouillier | Finding at least one point in each connected component of a real algebraic set defined by a single equation[END_REF]Lemma 3.5], V ε is a smooth algebraic set in C⟨ε⟩ n .

Let

H ε = V ε ∩ R⟨ε⟩ n .
x ε = x, let C ε be the connected component of H ε containing x ε . Then, C ε is bounded over R.
[17, Lemma 2] that we cite below gives explicitly a subset of V that contains I (H).

Proposition 3 ([17, Lemma 2]). Assume that I (H) is not empty and let x ∈ I (H).

There exists a semi-algebraically connected component

C ε of H ε such that C ε is bounded over R and lim ε C ε = {x}. Consequently, for 1 ≤ i ≤ n, there exists a point x ε ∈ crit(π i , V ε ) ∩ C ε such that lim ε x ε = x.
Then, we have

I (H) ⊂ ∩ n i=1 lim ε crit(π i , V ε ) ∩ R n .
To obtain a finite set of candidates, we need to introduce a generic linear change of variable to f . For any

A ∈ GL(n, Q), let f A (x) = f (A • x) and V A = V( f A ). For a generically chosen A, the set ∩ n i=1 lim ε crit(π i , V A ε ) is finite by [22, Theorem 1]. Since I (H)
is the image of the real isolated solutions of f A = 0 by the linear map induced by A -1 , we have the proposition below.

A ∩ GL(n, Q), the image C of ∩ n i=1 lim ε crit(π i , V A ε )
by the linear map induced by A -1 is finite and contains I (H).

Identification of the real isolated points

Once the set C of candidates is computed, we need to identify whether a candidate η = (η 1 , . . . , η n ) ∈ C ∩ R n is an isolated point of H. To do so, one can check whether a sphere centered at η of infinitesimal radius intersects H, which leads one to solve the system

f (x 1 , . . . , x n ) = n i=1 (x i -η i ) 2 -ε = 0
over R⟨ε⟩. This computation involves an infinitesimal ε, which would prevent a practically efficient algorithm. The objective of this subsection is to present a workaround to avoid the use of infinitesimals.

Let a = (a 1 , . . . , a n ) be an n-uple of positive rational numbers. We consider the function d, depending on a, defined by

d : R n × R n → R (x 1 , . . . , x n , y 1 , . . . , y n ) → n i=1 a i (x i -y i ) 2 .
The function d defines a metric in R n which can be extended to R⟨ε⟩ n . Further, the notations below denote respectively spheres, open balls and closed balls with respect to the metric d:

• S (x, r) = {y ∈ R n | d(x, y) = r}, • B(x, r) = {y ∈ R n | d(x, y) < r}, • B(x, r) = {y ∈ R n | d(x, y) ≤ r}.
For each x = (x 1 , . . . , x n ) ∈ C⟨ε⟩ n , we consider the function d x below:

d x : C⟨ε⟩ n → C⟨ε⟩ (y 1 , . . . , y n ) → n i=1 a i (y i -x i ) 2 .
Recall that the algebraic set V ε defined by

f 2 = ε 2 in C⟨ε⟩ n , is smooth. By an algebraic variant of Sard's theorem (see [25, Prop. B.2], [5, Theorem 9.6.2]), the critical values of the restriction of d x to V ε form a finite subset of C⟨ε⟩. Therefore, for every candidate x ∈ C ∩ R n , the set D(x) = {d x (lim ε y ε ) | y ε ∈ crit(d x , V ε ) ∩ R⟨ε⟩ n b , lim ε y ε x}
is a finite set of positive elements of R. When D(x) is not empty, Lemma 5 below allows us to identify whether the point x is isolated in H.

Lemma 5. Let x ∈ H and C x be the connected component of H containing x. Assume that the set D(x) defined as above is not empty. Let e x ∈ R such that

0 < e x < min D(x).
Then, the following statements are equivalent:

i) x is an isolated point of H.
ii) There exists e ∈]0,

e x [ such that H ∩ S (x, √ e) = ∅. iii) For every e ∈]0, e x [, H ∩ S (x, √ e) = ∅.
Moreover, if x is not an isolated point of H, then C x intersects S (x, √ e) for every e ∈]0, e x [.

Proof. By the definition of real isolated points, we immediately have that (i) implies (ii) and (iii) implies (i). It remains to demonstrate that (iii) is a consequence of (ii), which we separate into two statements: (ii) leads to (i) and then (i) leads to (iii). We now show that (ii) implies (i) by contradiction. We assume that the point x is not an isolated point of H.

If C x is not bounded, C x intersects S (x, √ e)
for every e > 0 and there is nothing to be proved. We now assume that C x is bounded.

By Lemma 1, there exists a point x ε such that x ε is bounded over R and lim

ε x ε = x. Let C ε ⊂ R⟨ε⟩ n be a connected component of the real algebraic set H ε containing x ε . By Lemma 2, C ε is bounded over R and, thus, its limit is connected in R n by [2, Proposition 12.49]. Hence, lim ε C ε is a connected subset of C x .
Moreover, as C ε is closed and bounded over R, the restriction of d x to C ε reaches its maximum at some point

y ε ∈ C ε (see [5, Theorem 2.5.8]). Note that y ε ∈ crit(d x , V ε ) ∩ R⟨ε⟩ n
b . So, we have that d x (lim ε y ε ) ≥ e x . Therefore, for any e ∈]0, e x [, the closed ball B(x, √ e) does not contain lim ε y ε . Since lim ε C ε is connected in R n and contains x and lim ε y ε , there exists a semi-algebraic continuous function γ : [0, 1] → lim ε C ε such that γ(0) = x and γ(1) = lim ε y ε . As d x (x) = 0 and d x (lim ε y ε ) ≥ e x , by the intermediate value property [2, Proposition 3.5], there exists t 0 ∈]0, 1[ such that d x (γ(t 0 )) = e for any e ∈]0, e x [. Therefore, the connected component C x intersects S (x, √ e) at γ(t 0 ). So, (ii) does not hold either when C x is bounded. We conclude that (ii) leads to (i). Finally, it remains to show that (i) implies (iii). Again, we prove this by contradiction. Assume that there exists e ∈]0, e x [ such that H ∩ S (x, √ e) ∅. Equivalently, there exists a point z ∈ H such that d(z) = e ∈]0, e x [. By Lemma 1, there exists a point

z ε ∈ H ε such that lim ε z ε = z. Let C z,ε be the connected component of H ε containing z ε .
In the closed and connected semi-algebraic set C z,ε , there exists a point z ′ ε at which the restriction of

d x to C z,ε reaches it minimum. So, z ′ ε belongs to crit(d x , V ε ) ∩ R⟨ε⟩ n b . Thus, we have that d x (lim ε z ′ ε ) ≤ d x (lim ε z ε ) < e x .
Using the definition of e x , we deduce that

d x (lim ε z ′ ε ) = 0, which is equivalent to lim ε z ′ ε = x. So, both z ε and z ′ ε lie in the connected component C z,ε of H ε and lim ε z ′ ε = x. If C z,ε
is not bounded over R, by Lemma 2, we have that C x is not bounded, which implies immediately that x is not an isolated point of H.

Otherwise, when C z,ε is bounded over R, by [2, Proposition 12.49], lim ε C z,ε is a connected subset of H that contains x and z. In this case, we also conclude that x is not isolated in H. Therefore, (i) leads to (iii), which finishes our proof. Lemma 6 handles the remaining case when a point x ∈ H such that D(x) is empty exists. Lemma 6. Assume that there exists x ∈ H such that D(x) is empty. Then, exactly one among the two statements below holds: i) H is connected and not bounded, so it does not have any isolated point.

ii) H is a single point x.

Proof. Assume that H has at least two connected components. Then, there exists a connected component of H ε that does not contain any point whose limit is x. Therefore, the restriction of d x to this connected component admits a critical point over this connected component. In consequence, D(x) is not empty, which contradicts the assumption of Lemma 6. Therefore, H has exactly one connected component.

Assume now if H is bounded and is not a single point x. As a consequence of Lemma 2, there exists a connected component C x,ε , that is bounded over R, of H ε such that {x} ⊊ lim ε C x,ε . Thus, the distance function δ x attains its maximum, which is non-zero, over H. This contradicts the assumption that D(x) is empty. Thus, we conclude the proof.

When D(x) is empty, we define by convention min

D(x) = +∞. Let e 0 ∈ R such that 0 < e 0 < min x∈C∩R n min D(x).
We deduce from Lemmas 5 and 6 the following proposition, which is the main criteria for designing our algorithms in Section 4.

Proposition 7. Let a = (a 1 , . . . , a n ) ∈ R n + and e 0 ∈ R defined as above. Then, for any candidate x = (x 1 , . . . , x n ) ∈ C ∩ R n , x is an isolated point of H if and only if the polynomial system f (y 1 , . . . ,

y n ) = n i=1 a i (y i -x i ) 2 -e 0 = 0 admits at least one solution in R n .
Proof. We first assume that D(x) is not empty. Using Lemma 5, we have that x is an isolated point of H if and only if H intersects the sphere S (x, √ e 0 ). Otherwise, if D(x) is empty, by Lemma 6, H is either a single point x or an unbounded connected component containing x. The similar conclusion follows immediately.

Algorithms

Overview of the algorithms

Our algorithms take as input f ∈ Q[x 1 , . . . , x n ] and compute:

• A zero-dimensional parametrization C = (w, v 1 , . . . , v n ) encoding the candidates such that I (H) ⊂ {(v 1 (t), . . . , v n (t)) | w(t) = 0};
• A finite set B of intervals isolating the real solutions of w(t) which correspond to I (H).

Both algorithms share the first step of computing the set of candidates C. This is done by the subroutine Candidates that takes as input the polynomial f and returns a zero-dimensional parametrization C encoding C. The design of this subroutine is given in [17, Subsection 3.2], which mimics the computation in [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]. Note that Candidates is probabilistic. As explained in Proposition 4, it computes the candidates for f A = f (A • x) where A is randomly chosen from GL(n, Q) then reverses the change of variables to obtain the candidates for f . Next, we pick randomly an n-tuple a = (a 1 , . . . , a n ) from Q n + , Proposition 7 requires us to compute a value e 0 ∈ Q such that for every

x ∈ C ∩ R n , 0 < e 0 < min{d x (lim ε y ε ) | y ε ∈ crit(d x , V ε ) ∩ R⟨ε⟩ n b , lim ε y ε x},
where d x is defined as

(y 1 , . . . , y n ) → a 1 (x 1 -y 1 ) 2 + . . . + a n (x n -y n ) 2 .
We call ComputeE0 a subroutine that takes as input f ∈ Q[x 1 , . . . , x n ] and a ∈ Q n + and returns such an e 0 ∈ Q + . The explicit description of ComputeE0 is given in Subsection 4.2. From this value e 0 , Proposition 7 identifies whether a candidate η ∈ C ∩ R n is isolated in H by deciding the emptiness of H ∩ S (η, √ e 0 ). This leads us to solve the following polynomial system over R n for each candidate η

= (η 1 , . . . , η n ) ∈ C ∩ R n : f (x 1 , . . . , x n ) = a 1 (x 1 -η 1 ) 2 + . . . + a n (x n -η n ) 2 -e 0 = 0. (2)
Therefore, Algorithm 1 below gives the outline of our algorithms in which the subroutine IsIsolated that takes as input f , C , a and e 0 and computes the isolating intervals B. It will be designed differently for each of our algorithms.

Algorithm 1: IsolatedPoints

Input:

A polynomial f ∈ Q[x 1 , . . . , x n ] Output: A zero-dimensional parametrization C and a set B of intervals in R 1 A chosen randomly from GL(n, Q) 2 C ← Candidates( f, A) 3 a chosen randomly in Q n + 4 e 0 ← ComputeE0( f, C , a) 5 B ← IsIsolated( f, C , a, e 0 ) 6 return (C , B)
Since the candidates are given by C , we cannot treat directly the system (2) and need some workarounds. The first variant of IsIsolated considers the system

f (x 1 , . . . , x n ) = a 1 (x 1 -v 1 (t)) 2 + . . . + a n (x n -v n (t)) 2 -e 0 = w(t) = 0. ( 3
)
We need to identify for which real roots of w(t), the above system has at least one real solution.

To do so, we basically compute a polynomial Q(t, z) ∈ Q[t, z] such that for any real root t 0 of w, Q(t 0 , z) admits real solutions if and only if the system (3) does too. The problem is hence reduced to a bivariate setting, which can be solved easily by classical real root counting algorithms. The design details of this variant is explained in Subsection 4.3.

The second variant of IsIsolated is explained in Subsection 4.4. The idea is to take advantage of the knowledge of e 0 , we replace the candidate η in the system (2) by an "approximation" η ∈ Q n of η and establish a similar result as Proposition 7 for these approximations (see Lemma 10). Briefly, we claim that a candidate η is an isolated point of H if and only if the system

f (x 1 , . . . , x n ) = a 1 (x 1 -η1 ) 2 + . . . + a n (x n -ηn ) 2 - e 0 4 = 0
has no real solution. Therefore, once those approximations are identified, one can apply classical real root finding algorithms to the above system.

Computing an appropriate value for e 0

In this subsection, we describe a subroutine that computes a value e 0 introduced in Proposition 7. Recall that, for each η = (η 1 , . . . , η n ) ∈ C n , the function d η is defined as

d η : C n → C, x = (x 1 , . . . , x n ) → n i=1 a i (x i -η i ) 2 .
To apply Lemma 5, we need to compute a value e 0 ∈ Q such that for every η ∈ C ∩ R n 0 < e 0 < min{d η (lim

ε x ε ) | x ε ∈ crit(d η , V ε ) ∩ R⟨ε⟩ n b , lim ε x ε η}.
Lemma 8 shows that, for a generic choice of a, every critical locus crit(d η , V ε ) is finite.

Lemma 8. Let η ∈ C be a candidate. Then there exists a non-empty Zariski open subset A of C n such that, for a = (a 1 , . . . ,

a n ) ∈ A ∩ Q n + , the critical locus crit d η , V ε is finite. Proof. Since V( f -ε) ⊂ C⟨ε⟩ n is smooth, the critical locus crit d η , V( f -ε) is defined by f -ε = y • ∂ f ∂x 1 -2a i (x 1 -η 1 ) = • • • = y • ∂ f ∂x n -2a n (x n -η n ) = 0. ( 4
)
Now we consider a = (a 1 , . . . , a n ) and y as indeterminates. Let φ : C n × C n × C be the polynomial mapping defined as

(x, a, y) → f -ε, y • ∂ f ∂x 1 -2a 1 (x 1 -η 1 ), . . . , y • ∂ f ∂x n -2a n (x n -η n ) .
Let X be the non-empty Zariski open subset of C defined as

(x 1 -η 1 ) • • • (x n -η n ) 0.
The Jacobian matrix of φ with respect to (x, a, y)

                    ∂ f ∂x 1 • • • ∂ f ∂x n 0 0 • • • 0 * • • • * ∂ f ∂x 1 -2(x 1 -η 1 ) • • • 0 . . . . . . * . . . . . . . . . . . . * • • • * ∂ f ∂x n 0 • • • -2(x n -η n )                    
has full rank when x ∈ X and

f -ε = y • ∂ f ∂x 1 -2a i (x 1 -η 1 ) = • • • = y • ∂ f ∂x n -2a n (x n -η n ) = 0.
By Thom's weak transversality theorem [8, Theorem 3.7.4], there exists a non-empty Zariski open subset A ∅ of C n such that, for a ∈ A ∅ , 0 is a regular value of the restriction of φ a to X. Thus, for a ∈ A ∅ , by Jacobian criterion, the restriction of the solutions of

f -ε = y • ∂ f ∂x 1 -2a i (x 1 -η 1 ) = • • • = y • ∂ f ∂x n -2a n (x n -η n ) = 0 to X is a finite set, i.e., crit(d η , V( f -ε)) ∩ X is finite. Now we study the restriction of crit(d η , V( f -ε)) to C n \ X. We choose a ∈ Q n + .
Let I be a non-empty proper subset of {1, . . . , n} and X I be the subset of C n defined by

x i = η i for i ∈ I and x i η i for i {1, . . . , n} \ I. Let x ∈ crit(d η , V( f -ε)) ∩ X I . As f (η) = 0, η V( f -ε) and x η.
Hence, y 0 in the system (4). Hence, y Since the set of candidates C is encoded by a zero-dimensional parametrization C , we do the whole computation at once through the function δ defined as

• ∂ f ∂x i -2a 1 (x i -η i ) = 0 implies that ∂ f ∂x i = 0. Since V( f -ε) is smooth, ∂ f ∂x i for i ∈ {1, . . . , n} \ I cannot vanish simultaneously at x. This means {x | x ∈ crit(d η , V( f -ε)), x i = η i for i ∈ I}
δ : C n × C → C, (x 1 , . . . , x n , t) → n i=1 a i (x i -v i (t)) 2 .
The following lemma is immediate.

Lemma 9. Let V ε,t ⊂ C⟨ε⟩ n+1 be the algebraic set defined by f 2 -ε 2 = w(t) = 0. Then, the set of critical values δ(crit(δ, V ε,t )) is the union of d η (crit(d η , V ε )) for η ∈ C.
Proof. The set crit(δ, V ε,t ) are defined by the points of V ε,t at which the matrix

           ∂ f 2 ∂x 1 . . . ∂ f 2 ∂x n 0 ∂δ ∂x 1 . . . ∂δ ∂x n ∂δ ∂t 0 . . . 0 w ′ (t)            has rank at most 2.
As w(t) is square-free, for every t 0 such that w(t 0 ) = 0, w ′ (t 0 ) is not zero. Therefore, the condition above restricted to V ε,t is equivalent to rank

       ∂ f 2 ∂x 1 . . . ∂ f 2 ∂x n ∂δ ∂x 1 . . . ∂δ ∂x n        ≤ 1.
For every complex root t 0 of w(t), let η 0 = (v 1 (t 0 ), . . . , v n (t 0 )). By fixing t = t 0 , the rank condition above is reduced to rank

        ∂ f 2 ∂x 1 . . . ∂ f 2 ∂x n ∂δ η 0 ∂x 1 . . . ∂δ η 0 ∂x n         ≤ 1, which defines the set crit(δ η 0 , H ε ). Thus, crit(δ, H ε,t ) = ∪ w(t 0 )=0 {(x, t 0 ) | η 0 = (v 1 (t 0 ), . . . , v n (t 0 )), x ∈ crit(δ η 0 , H ε )}.
This concludes the proof. Now we aim to compute the limit of the critical points and their corresponding values of the restriction of δ to the algebraic set H ε,t defined by f 2 -ε 2 = w(t) = 0. Note that Lemmas 8 and 9 imply that, for a generic a ∈ Q n + , the set of critical points crit(δ, H ε,t ) is finite. By [22, Theorems 1, 2], for generic values of a ∈ Q n , we have that

lim ε crit(δ, H ε,t ) ⊂ ⟨ f ⟩ + w(t), y • ∂ f ∂x i - ∂δ ∂x i for every 1 ≤ i ≤ n ∩ Q[x, t]
and the ideal on the right-hand side is zero-dimensional. From the above inclusion, one can follow a similar computation as in [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF] using the geometric resolution algorithm [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF]. However, as the degree of w(t) is bounded by 2D(D -1) n-1 (see Section 5), such a computation would lead to an arithmetic complexity D O(n 2 ) .

A workaround to obtain a better complexity is to use a variant of geometric resolution over the quotient ring A = Q[t]/⟨w(t)⟩ as explained in [25, Appendix J]. Note that w(t) is not necessarily irreducible, the extension A is only a product of fields and doing the computation over the ring A is not trivial. We will see in Subsection 5 that this approach allows us to obtain an algorithm with arithmetic complexity lying in D O(n) .

Our subroutine ComputeE0 is designed as follows.

a) First, we call a subroutine ParametricCurve that takes as input f , C , a ∈ Q n + and i ∈ {1, . . . , n} and computes a one-dimensional parametrization

J i over Q[t]/⟨w(t)⟩ of the system ∂δ ∂x j • ∂ f ∂x i - ∂δ ∂x i • ∂ f ∂x j = 0 j∈{1,...,n}\{i} and ∂ f ∂x i 0.
An explicit description of this subroutine can be found in [25, Appendix J.5].

b) Next, we compute a zero-dimensional parametrization E i of the intersection of H = V( f ) with the sets of solutions defined by the parametrizations J i above.

This is done by calling a subroutine IntersectCurve on the input J i and f , which is described also in [25, Appendix J.5].

c) We then call a subroutine Union that computes a zero-dimensional parametrization E that defines ∪ n i=1 Z(E i ).

d) Finally, taking as input the zero-dimensional parametrization E , we call a subroutine GetE0 that computes the required value e 0 . This can be done by calling FGLM algorithm [START_REF] Faugère | Efficient computation of zero-dimensional gröbner bases by change of ordering[END_REF] to compute a polynomial P(e) whose solutions encode the values e = i=1 a i (x iv i (t)) 2 for x ∈ Z(E ) and w(t) = 0. Next, we evaluate a lower bound of the minimal distance between the roots of P(e) using [2, Proposition 10.23].

Algorithm 2: ComputeE0 Input: f ∈ Q[x 1 , . . . , x n ], C = (w(t), v 1 (t), . . . , v n (t)) and a ∈ Q n + Output: e 0 ∈ Q 1 δ ← a 1 (x 1 -v 1 (t)) 2 + . . . + a n (x n -v n (t)) 2 2 for 1 ≤ i ≤ n do 3 J i ← ParametricCurve( f, a, C , i) 4 E i ← IntersectCurve(J i , f ) 5 E ← Union(E 1 , . . . , E n ) 6 e 0 ← GetE0(E ) 7 return e 0

The first variant of IsIsolated

In this subsection, we explain the details of the first variant of IsIsolated.

Using the value e 0 output by Algorithm 2, Proposition 7 allows one to identify the isolated points of H among the candidates by checking whether the polynomial system 2 e 0 = 0 admits real solutions for each candidate η = (η 1 , . . . , η n ) ∈ C ∩ R n . Again, one can consider the system

f (x 1 , . . . , x n ) = n i=1 a i (x i -η i )
f (x 1 , . . . , x n ) = n i=1 a i (x i -v i (t)) 2 -e 0 = w(t) = 0 (5)
to handle all the candidates at once. Let W t ⊂ C n+1 be the algebraic set defined the equation ( 5). Our strategy is to compute a finite subset of W t ∩ R n+1 that intersects every connected component of W t ∩ R n+1 . Then, all the real t-coordinates of those sample points correspond to the isolated points of H ∩ R n .

We consider the polynomial

F = f (x 1 , . . . , x n ) 2 +        n i=1 a i (x i -v i (t)) 2 -e 0        2 .
Note that F + w(t) 2 defines also the real algebraic set W t ∩ R n+1 . Therefore, the sample points above can be computed using the algorithm of [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF] on the input F +w(t) Let B be a matrix randomly chosen from GL(n, Q) and F B (x) = F(B • x). We apply the geometric resolution algorithm over A on the system of equations:

F B = 0, ∂F B ∂x j = 0, ∂F B ∂x 1 0.
This algorithm returns a zero-dimensional parametrization (U(z), V 1 (z), . . . , V n (z)) over the ring A, which means that V 1 , . . . , V n and U are elements of A[z], such that, for any real solution t 0 of w(t), the finite set defined by

{(V 1 (t 0 , z), . . . , V n (t 0 , z)) | z ∈ R, U(t 0 , z) = 0}
intersects every connected component of

f (x 1 , . . . , x n ) = n i=1 a i (x i -v i (t 0 )) 2 -e 0 = 0.
Hence, the isolated points of H ∩ R n are indeed

{(v 1 (t), . . . , v n (t)) | (t, z) ∈ R 2 : w(t) = U(t, z) = 0}.
Our problem boils down to solving the bivariate system w(t) = U(t, z) = 0 over R 2 .

In Algorithm 3 below, we introduce two subroutines:

• BivariatePolynomial takes as input the polynomials F and w(t) and returns the eliminating polynomial U(t, z). It uses the geometric resolution algorithm over A described in [25, Appendix J].

• BivariateSolve takes as input w(t) and U(t, z) and returns the set B of intervals that isolate the real roots of w(t) corresponding to I (H). Such a subroutine can be designed efficiently with resultants.

Algorithm 3: The first variant of IsIsolated Input: f , C = (w(t), v 1 (t), . . . , v n (t)), a and e 0 Output: A set B of isolating intervals

1 F ← f (x 1 , . . . , x n ) 2 + n i=1 a i (x i -v i (t)) 2 -e 0 2 2 U(t, z) ← BivariatePolynomial(F, w(t)) 3 B ← BivariateSolve(w(t), U(t, z)) 4 return B

Approximations of the candidates

This subsection describes the design of the second variant of IsIsolated. This variant does not require solving polynomial systems in the quotient ring Q[t]/⟨w(t)⟩ but is based mostly on isolating the candidates from the parametrization C . The main idea is to replace the candidate η in the criteria provided by Proposition 7 by a rational approximation η ∈ Q n . This allows the subroutine presented in Algorithm 4 to involve only the real points defining by C and not its whole complex solution set.

To compute the approximations, we first identify how close the points η and η need to be. The lemma below shows that requiring d a (η, η) < √ e 0 /2 is enough.

Lemma 10. Let η ∈ C ∩ R n and η be a point in R n satisfying d a (η, η) < √ e 0 /2. Then, η is an isolated point of H if and only if H does not intersect the sphere S ( η, √ e 0 /2).

Proof. If the set {x ε ∈ crit(δ η , V ε ) ∩ R⟨ε⟩ n b , lim ε x ε η} is empty, then, by Lemma 6, H is either a single point η or an unbounded connected set containing η. In either case, the conclusion of Lemma 10 is immediate. Thus, in what follows, {x ε ∈ crit(δ η , V ε ) ∩ R⟨ε⟩ n b , lim ε x ε η} is assumed to be non-empty.

We prove now the necessary implication. Assume that η is an isolated point of H. By Lemma Now we turn to the sufficient implication. Assume by contradiction that η is not isolated in H. By Lemma 5, the connected component C η of H containing η intersects the sphere S (η, √ e 0 ).

So, there exists a semi-algebraic continuous function γ : [0, 1] → C η such that γ(0) = η and γ(1) lying on the sphere S (η, √ e 0 ). We have that

d a (γ(1), η) ≥ d a (γ(1), η) -d a (η, η) > √ e 0 - √ e 0 /2 = √ e 0 /2.
As d a (γ(0), η) < √ e 0 /2 and d a (γ(1), η) > √ e 0 /2, by the intermediate value property [2, Proposition 3.5], there exists t 0 ∈]0, 1[ such that d a (γ(t 0 ), η) = √ e 0 /2. This implies that the intersection of H and S ( η, √ e 0 /2) is not empty, which concludes our proof.

Let t η be the real root of w(t) corresponding to η, i.e., η = v 1 (t η ), . . . , v n (t η ) . To apply Lemma 10, we need to choose t η ∈ Q such that the rational point η = (v 1 (t η), . . . , v n (t η)) satisfies that d a (η, η) < √ e 0 /2. This leads us to identify ρ > 0 such that |t ηt η| < ρ implies

a 1 (v 1 (t η ) -v 1 (t η)) 2 + . . . + a n (v n (t η ) -v n (t η)) 2 < e 0 /4.
Lemma 11 below allows us to compute explicitly an appropriate value for ρ.

Lemma 11. Let {t 1 , . . . , t ℓ } be the distinct real roots of w(t) = 0 and {η 1 , . . . , η ℓ } be the corresponding candidates. We consider a set of intervals (I j ) 1≤ j≤ℓ such that

• The intervals I j are pairwise disjoint.

• The interval I j contains only t j as a real root of w(t).

For each

1 ≤ i ≤ n, let K i = max ℓ j=1 max t∈I j |v ′ i (t)|.
Then, for any 1 ≤ j ≤ ℓ and t θ such that t θ ∈ I j and |t θ -

t j | < 1 K i • e 0
4na i , we have the following inequality:

|v i (t θ ) -v i (t j )| < e 0 4na i . Let ρ ≤ min n i=1 1 K i • e 0 4na i . For any real root t η of w(t) and t θ ∈ I j , |t θ -t η | < ρ implies d a (θ, η) < √ e 0 /2.
Proof. For 1 ≤ j ≤ ℓ and any t θ ∈ Q, we have that

v i (t θ ) -v i (t j ) = v ′ i ( t j )(t θ -t j ),
where t ∈ R lies between t θ and t j . For t ∈ I j =]r j , s j [, by the definition of K i , we have

|v ′ i (t)| ≤ K i . Then, for t θ ∈ I j such that |t θ -t j | < 1 K i • e 0 4na i , we have |v i (t θ ) -v i (t j )| = |v ′ i ( t j ) • (t θ -t j )| ≤ K i • |t θ -t j | < e 0 4na i . Now we take ρ ≤ min n i=1 1 K i • e 0 4na i . If t θ ∈ I j and |t θ -t j | < ρ, then we have d a (θ, η j ) = n i=1 a i (v i (t θ ) -v i (t j )) 2 < n i=1 e 0 4n = √ e 0 2 .
Lemmas 10 and 11 provides us the ingredients to design Algorithm 4. It requires us to introduce two subroutines Isolate and MaxOverInterval below.

• We need two versions of Isolate. The first one takes as input a polynomial p ∈ Q[t] and returns a set of disjoint intervals of rational extremities isolating the real roots of p.

Besides the polynomial p ∈ Q[t], the second version of Isolate requires a positive ρ ∈ Q as input and returns the intervals of length at most ρ that isolate the real roots of p.

The explicit descriptions of both of these real root isolating algorithms are given in [START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF].

• MaxOverInterval takes as input a polynomial p ∈ Q[t] and an interval [r, s] where r, s ∈ Q and returns an upper bound of max t∈[r,s] |p(t)|. Such a subroutine can be implemented using the following naive bound:

max t∈[r,s] |p(t)| ≤ deg(p) i=0 |c i | where p t-r s-r = c 0 • t deg(p) + • • • + c deg(p)
. Algorithm 4 proceeds through these following steps: a) We call Isolate on the input w(t) to obtain a set of intervals I j that isolate the real roots of w(t) and compute

K i = max ℓ j=1 max t∈I j |v ′ i (t)| using the subroutine MaxOverInterval on the input v ′ i (t) and each interval I j . b) We then compute ρ ∈ Q such that 0 < ρ ≤ min n i=1 1 K i • e 0
4na i and use Isolate on the polynomial w(t) and the precision ρ to obtain a set of intervals Ĩ j such that each Ĩ j contains exactly one real root of w(t) and | Ĩ j | < ρ. c) For 1 ≤ j ≤ ℓ, we choose a point t j in I j ∩ Ĩ j ∩ Q and evaluate η j = (v 1 ( t j ), . . . , v n ( t j )).

The set C of the approximations is taken as {( η j , I j | 1 ≤ j ≤ ℓ}. d) Finally, we decide whether the system 2 e 0 /4 = 0 has a real solution for each approximation η and return those which do not.

f (x 1 , . . . , x n ) = n i=1 a i (x i -ηi )
We summarize Section 4 in Algorithm 4 below, which is our second variant of IsIsolated. Note that the sphere in Lemma 10 can be replaced by any hypercube that contains the candidate η in its interior and is contained in the ball B(η, e 0 ). This leads one to decide the emptiness of semi-algebraic sets defined by f = 0 and some linear polynomial inequalities instead of the quadric defining the sphere. We observe in practice that this helps reduce the bit-size growth and accelerates the computation.

Algorithm 4: Algorithm IsIsolated-Approx Input: f , C = (w(t), v 1 (t), . . . , v n (t)), a ∈ Q n + and e 0 ∈ Q Output: A set of isolating interval B 1 {I 1 , . . . , I ℓ } ← Isolate(w(t)) 2 for i ∈ {1, . . . , n} do 3 K i ← max ℓ j=1 MaxOverInterval(v ′ i (t), I j ) 4 { Ĩ1 , . . . , Ĩℓ } ← Isolate w(t), ρ = min n i=1 1 K i • e 0 4na i 5 for j ∈ {1, . . . , ℓ} do 6 t j ∈ I j ∩ Ĩ j 7 η j ← (v 1 ( t j ), . . . , v n ( t j )) 8 C ← {( η j , I j ) | 1 ≤ j ≤ ℓ}, B ← ∅ 9 for ( η, I η ) ∈ C do

Complexity analysis

We establish now the complexity results for two variants of Algorithm 1, starting with the one using Algorithm 3 for IsIsolated.

Theorem I. Let f ∈ Q[x 1 , . . . , x n ].
Then, the variant of Algorithm 1 which uses Algorithm 3 computes the real isolated points of the algebraic hypersurface defined by f within O 64 n D 8n arithmetic operations in Q and one call of real root isolation on a univariate polynomial of degree bounded by 2 n+2 D 2n .

Proof. We start with the subroutine Candidates. Since crit(π i , H A ε ) is finite and defined by

( f A -ε) • ( f A + ε) = 0, ∂ f A ∂x j = 0 for all j i,
its degree is bounded by 2D(D -1) n-1 using Bézout bound. Consequently, the degree of the output zero-dimensional parametrization is bounded by 2D(D -1) n-1 . Using [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]Theorem 4] (which is based on the geometric resolution algorithm in [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF]), each zero-dimensional parametrization of crit(π i , H A ε ) is computed within O D 3n arithmetic operations in Q. The last step which takes intersections of those parametrizations is done using the algorithm in [25, Appendix J.1] ; it does not change the asymptotic complexity.

Hence, computing the parametrization C encoding the candidates can be done within O D 3n arithmetic operations in Q and the degrees of the polynomials w(t), v 1 (t), . . . , v n (t) are bounded by 2D(D-1) n-1 . It remains to estimate the arithmetic complexity of the subroutines ComputeE0 and IsIsolated.

Let κ be the degree of w(t). Algorithm 2 (ComputeE0) relies on computing the limit of crit(δ, H ε,t ) ∩ C⟨ε⟩ n b , where H ε,t is the algebraic set defined by

( f -ε) • ( f + ε) = 0, w(t) = 0 ( 6
)
and δ is the distance function

(x 1 , . . . , x n ) → n i=1 a i (x i -v i (t)) 2
We use the algorithm of [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF] on the function δ for the resolution of polynomial systems in the quotient ring Q[t]/⟨w(t)⟩. Using Bézout's bound on the system

f 2 -ε 2 = y • ∂ f ∂x 1 - ∂δ ∂x 1 = • • • = y • ∂ f ∂x n - ∂δ ∂x n = 0
defining crit(δ, H ε,t ) over A, the degree of crit(δ, H ε,t ) in C⟨ε⟩ n is bounded by 2D n+1 κ ≤ 4D n+2 (D-1) n-1 ≈ 4D 2n+1 . By [25, Appendix J.5], the arithmetic operations over A can be done using O (κ) operations in Q. Thus, applying [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]Theorem 5], we obtain the complexity bound O κ • D 3n+2 ≈ O D 4n+2 for obtaining the zero-dimensional parametrization E in Algorithm 2.

The call to GetE0 computes from the zero-dimensional parametrization E a univariate polynomial P(e) ∈ Q[e] whose solutions are the critical values of δ restricted to V ε . Since the degree of P(e) is bounded by 4D 2n+1 , this can be done using FGLM algorithm [START_REF] Faugère | Efficient computation of zero-dimensional gröbner bases by change of ordering[END_REF] within O D 6n+3 arithmetic operations over Q. Next, it computes the minimal distance between the real roots of P(e) using [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]Proposition 10.23]. The complexity of this computation is linear in the degree of P(E). Thus, it does not change the asymptotic complexity of Algorithm 2.

Therefore, Algorithm 2 can be done within O (D 6n+3 ) arithmetic operations in Q.

Algorithm 3 is basically computing sample points of the hypersurface

F = f (x 1 , . . . , x n ) 2 +        n i=1 a i (x i -v i (t)) 2 -e 0        2
over the quotient ring Q[t]/⟨w(t)⟩. Again, we follow the algorithm of [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF] on the input F with the extended version of geometric resolution to the quotient ring A. By [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]Theorem 6] with the overcost O (κ) of arithmetic operations over A, we obtain the complexity bound

O κ • (2D) 3n+2 ≈ O 8 n D 4n+2
for MinimalPolynomial.

The output polynomial U(t, z) has degree at most (2D) n in z and κ in t so its total degree is bounded by (2D) n + κ. Therefore, solving the bivariate system w(t) = U(t, z) = 0 can be done within

O ((2D) n + κ) 4 κ 2 ((2D) n + κ) 2 ≈ O 64 n D 8n
arithmetic operations in Q using geometric resolution. In the end, one needs to isolate the real roots of the eliminating polynomial output by the geometric resolution. That polynomial has degree bounded by κ((2D

) n + κ) ≤ 2 n+2 D 2n .
Adding up all the steps, we obtain the arithmetic complexity of Algorithm 1, which lies in O 64 n D 8n with a call to real root isolation on a polynomial of degree bounded by 2 n+2 D 2n .

Note that for implementing our algorithm, we would mainly rely on Algorithm 4 (IsIsolated-Approx). Hence, we dedicate the rest of this subsection to discuss its complexity. The complexity result of our algorithm using Algorithm 4 is stated as follows.

Theorem II. Let f ∈ Q[x 1 , . . . , x n ].
Then, the variant of Algorithm 1 which uses Algorithm 4 requires O D 6n+3 arithmetic operations in Q and two real root isolating calls on a univariate polynomial of degree bounded by 2D(D -1) n-1 .

Proof. Recall that Algorithm 4 computes an approximation η = ( η1 , . . . , ηn ) for each candidate η ∈ C ∩ R n and decides whether the system

f (x 1 , . . . , x n ) = n i=1 a i (x i -ηi ) 2 - e 0 4 = 0
has a real solution. The arithmetic complexity for solving each of those decision problems lies in O 8 n D 3n+2 using [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]. Since the cardinality of C is bounded by 2D(D -1) Since each of ℓ evaluations η j ← v 1 ( t j ), . . . , v n ( t j ) takes O (nD n ) arithmetic operations, the cost of getting the approximations is bounded by O nD 2n . lim ε α lim ε β. Then, there exist two points x α and x β in crit(π j , V ε ) ∩ C i such that π j (x α ) = α and π j (x β ) = β. Since π j (lim ε x α ) π j (lim ε x β ), lim ε x α lim ε x β . Then lim ε x α and lim ε x β are two distinct points of C ∩ lim ε crit(π j , V A ε ). Consequently, there exists two distinct points in C ∩ lim ε crit(π j , V A ε ).

Lemma 14. For every bounded connected component C of H A that is not a singleton, there exist at least two points in n i=1 lim ε crit(π i , V A ε ) that belong to C.

Proof. Since C is not a singleton, there exists a coordinate x j such that the projection π j (C) is not a point. Using Lemma 13, we conclude the proof.

Let C 2 be the image of

n i=1 lim ε crit(π i , V A ε )
by the linear map induced by A -1 . Note that the set of candidates C is a subset of C 2 .

Proposition 15. Let C 2 be defined as above. Let x ∈ C ∩ R n and B ⊂ R n be a ball such that C 2 ∩ B = {x} and x is contained in the interior of B. Then, if the intersection of the boundary of B and H is empty, x is an isolated point of H.

Proof. Let C be the connected component of H containing x. If C is unbounded, then x is not an isolated point and the intersection of the boundary B with H is not empty. We now assume that C is bounded and is not a singleton. By Lemma 14, there exists y ∈ C 2 ∩ R n such that y x and y ∈ C. As x and y lie on different sides of B, by intermediate value theorem, H intersects the boundary of B, which ends the proof.

Since the subroutine Candidates in [17, Subsection 3.2] computes all the zero-dimensional parametrizations encoding lim ε crit(π i , V A ε ), the union C 2 can be obtained easily by taking the union of those zero-dimensional parametrizations. Next, we isolate the candidates in C ∩ R n such that each isolating ball contains exactly on point of C 2 ∩ R n .

Note also that for a generic matrix A ∈ GL(n, Q), the projection of any non-singleton bounded component of H A by π 1 is not a singleton. In this case, by Lemma 13, it is enough to take C 2 = lim ε crit(π 1 , V A ε ). This can accelerate implementations with a small cost of randomness (which is already a part of the subroutine Candidates).

Algorithm 5 contains the description of the subroutine SimpleIdentification. We call to a subroutine BoxIsolate that takes as input a zero-dimensional parametrization encoding a subset of C n and computes isolating boxes for its real zeros.

Algorithm 5: SimpleIdentification

Input: A zero-dimensional parametrization C 2 Output: A set B 1 of intervals of R 1 B 1 ← ∅ 2 Boxes ← BoxIsolate(C 2 ) 3 for box ∈ Boxes do 4 if box ∩ H = ∅ then 5 B 1 ← B 1 ∪ {t-coordinate of box} 6 return B 1 21
By Proposition 15, for each x ∈ C such that the intersection of the ball isolating x with H is empty, we conclude that x is an isolated point of H. For the non-empty intersections, we cannot decide whether x is isolated yet. The problem arises when the isolating boxes are not small enough so that they intersect not only the connected component of H containing x but also some other connected component. When this happens, one could try a smaller size of isolating boxes.

Remark 16. When every candidate is an isolated point, running SimpleIdentification is enough to return a correct and certified output, which is the whole set of candidates. No further computation is required in this case.

As each candidate is contained in lim ε crit(π i , V ε ) for every projection π i , we conjecture that every candidate is actually isolated. We are not aware of any counter-example.

Limits of critical curves

To compute a set of candidates, we consider the critical loci crit(π i , V ε ) for 1 ≤ i ≤ n. Our second optimization considers the critical loci of the projections on the plane; especially, the limits of those critical loci are curves in R n whose real isolated points contain the isolated points of H. Thus, one can compute a superset of I (H) through computing the real isolated points of limits of critical curves.

More precisely, for 1 ≤ i < j ≤ n, we denote by π i, j the projection π i, j : (x 1 , . . . , x n ) → (x i , x j ).

Recall that H ε is a smooth algebraic set defined by

( f -ε) • ( f + ε) = 0.
Lemma 17. Let A ∈ GL(n, Q). For every 1 ≤ i < j ≤ n, the set of isolated points of H A is contained in set of real isolated points of lim ε crit(π i, j , H A ε ). Proof. Let x be an isolated point of H A . By Proposition 3, x ∈ lim ε crit(π i , V A ε ). Since crit(π i , V A ε ) ⊂ crit(π i, j , V A ε ), x ∈ lim ε crit(π i, j , V A ε ). Note that lim ε crit(π i, j , V A ε ) is a subset of V A . Thus, if x is isolated in H A , it is also an isolated point in lim ε crit(π i, j , V A ε ) ∩ R n .

Remark 18. Note that a real isolated point of lim ε crit(π i, j , V ε ) is not necessarily isolated in H. Take for example the degenerate torus, given by the equation (x 2 1 + x 2 2 + x 3 ) 2 -4(x 2 1 + x 2 2 ) = 0. The real trace of lim ε crit(π 1,2 , V ε ) is the union of the point (0, 0) and the circle given by x 2 1 + x 2 2 -4 = x 3 = 0. Hence, Lemma 17 allows us to obtain a superset of I (H) only.

By [START_REF] El Din | Polar varieties and computation of one point in each connected component of a smooth real algebraic set[END_REF]Theorem 2], for a generic change of variables A, the critical locus crit(π i, j , H A ε ) is an equidimensional algebraic set of dimension one defined by

( f -ε) • ( f + ε) = 0, ∂ f ∂x k = 0 for 1 ≤ k ≤ n and k i, j.
The computation of lim ε crit(π i, j , H A ε ) can be done using a similar subroutine of [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]. For each 1 ≤ i, j ≤ n, we denote by J i, j the ideal

∂ f ∂x k = 0 for 1 ≤ k ≤ n, k i, j .
Lemma 19. Let π i, j be defined as above. There exists a non-empty Zariski open subset A of GL(n, C) such that, for any A ∈ A ∩ GL(n, Q), the algebraic set C defined by

V f A + J k : ∂ f A ∂x i ∞ ∩ J k : ∂ f A ∂x j ∞
is equi-dimensional of dimension 1 and contains lim ε crit(π i, j , H A ε ). As a consequence, any isolated point of H A is also isolated in C ∩ R n .

Proof. The proof of the first statement follows a similar outline of the proof of [22, Theorem 1 and Theorem 2]. From the inclusion

I (H A ) ⊂ lim ε crit π i, j , V A ε ⊂ C ∩ R n ⊂ H A ,
we deduce that every real isolated point of H A is also an isolated point of C ∩ R n .

We define the subroutine CurveLimitCheck that takes as input f ∈ Q[x 1 , . . . , x n ] and A ∈ GL(n, Q) and returns a set of isolating boxes B 2 . It calls to two subroutines:

• CurveLimit that takes as input f , A and a pair of index (i, j) and returns the eliminating polynomial of a rational parametrization encoding lim ε crit(π i, j , H A ε ). The design of this subroutine follows Lemma 19.

• BivariateIsolated that takes as input a bivariate polynomial U i, j and computes the boxes isolating the real isolated points of V(U i, j ). This can be done by computing a cylindrical algebraic decomposition adapted to U i, j = 0.

Algorithm 6: CurveLimitCheck Input: f ∈ Q[x 1 , . . . , x n ], A ∈ GL(n, Q) Output: A set B 2 of intervals of R 1 for 1 ≤ i < j ≤ n do 2 U i, j ← CurveLimit( f, A, (i, j)) 3 boxes i, j ← BivariateIsolated(U i, j ) computation is similar to the one of test1 with the main difference coming from the fact that the isolating boxes computed in approx requires more precision. This results in linear polynomials, that define hypercubes, of larger bit-sizes, which makes the column raglib slower than test1.

At the moment, we do not dispose of a geometric resolution algorithm for Q[t]/ ⟨w(t)⟩. The implementation of ComputeE0 relies on available tools such as FGb, msolve that run over Q. The complexity of this subroutine is actually bounded by D O(n 2 ) and the timings show that it is not practical. The value e 0 in these examples is obtained since we know in advance that the given real algebraic sets are finite.

n

Proposition 4 .

 4 There exists a Zariski open dense subset of GL(n, Q) such that for any A ∈

  coincides with the critical locus crit(d η,I , V( f I -ε)) where d η,I : (x j ) j∈{1,...,n}\I → j∈{1,...,n}\I a j (x j -η j )2 and f I is the polynomial obtained from f by substituting x i = η i for i ∈ I. Therefore, we can use the same arguments as above to prove the following.There exists a non-empty Zariski open subset A I of C n such that for a ∈ A I ∩ Q n + , the restriction of crit(d η , V( f -ε)) to X I is finite.Let A + = ∩ I⊊{1,...,n} A I which is a non-empty Zariski open subset of C n . Given any a ∈ A + ∩ Q n + , crit(d η , V( f -ε) is a finite set. Similarly for V( f + ε), we obtain a non-empty Zariski open subset A -. Taking the intersection A = A + ∩ A -ends the proof.

10 if 11 B

 1011 HasRealSolutions( η, f, a, e 0 ) = false then ← B ∪ I η 12 return B Remark 12.

  Two lemmas below will be used regularly in this paper. For every x ∈ H, there exists a point x ε ∈ H ε such that x ε is bounded over R and lim ε x ε = x. ∈ H ε such that x ε is bounded over R and lim ε

	Lemma 1. [20, Lemma 3.6] Lemma 2. [2, Proposition 12.51] Given a point x lying in a bounded connected component of H
	and x ε

  2 ⊂ Q[t, x 1 , . . . , x n ]. Such an algorithm returns a zero-dimensional parametrization over Q that defines a finite set intersects every connected component of W t . Since the total degree of F + w(t) 2 can go up to O(D n ), this computation faces the same complexity issue as in Subsection 4.2. Again, we can bypass this problem by solving over A[x 1 , . . . , x n ] where A = Q[t]/⟨w(t)⟩.

  5, the intersection of H and S (η, √ e) is empty for every e ∈]0, e 0 [. So, η is the only point of H lying in the open ball B(η, √ e 0 ). Since d a (η, η) < √ e 0 /2, the candidate η does not lie on

	the sphere S ( η, S ( η, √ e 0 /2) ∩ H = ∅. √ e 0 /2). Moreover, S ( η,	√ e 0 /2) is contained in the open ball B(η,	√ e 0 ). Then,

  It remains to estimate the complexity of computing the approximations, whose main steps consist of calling MaxOverInterval and isolating the real roots of the eliminating polynomial w(t) in the zero-dimensional parametrization encoding C. The subroutine MaxOverInterval is called n times for the polynomials v ′ i (t); this would require O (deg(w)) ≈ O (D n ) arithmetic operations over Q.

	n-1 , Algorithm 4
	runs within
	O 8 n D 4n+2
	arithmetic operations in Q.
	Note that all the complexities above are dominated by the complexity
	O D 6n+3
	of ComputeE0 (Algorithm 2).

Table 1 :

 1 cand cand2 |real sols.|/ deg(w) test1 total Experimental timings of Algorithm 7

	e0	approx raglib

B 2 ← ∩ 1≤i, j≤n boxes i, j

return B 2 Summary.. To conclude this section, we show below the pseudo-code of our implementation. The subroutine Candidates is modified so that it returns, besides C encoding the candidates, a zero-dimensional parametrization C 2 encoding the finite set C 2 .

(Huu Phuoc Le ) 1 Huu Phuoc Le is supported by the ANR grants ANR-18-CE33-0011 Sesame, the ANR-19-CE40-0018 De Rerum Natura, the joint ANR-FWF ANR-19-CE48-0015 ECARP project, the PGMO grant CAMiSAdo, the European Union's Horizon 2020 research and innovative training network program under the Marie Skłodowska-Curie grant agreement N°8 13211 (POEMA) and the Grant FA8665-20-1-7029 of the EOARD-AFOSR.

The real root isolation is called twice in Algorithm 4 on the polynomial w(t).

Summing up the above discussion, we conclude that Algorithm 4 requires

arithmetic operations in Q and two calls of real root isolation on a univariate polynomial of degree bounded by 2D(D -1) n-1 .

Furthermore, the complexity of real root isolation algorithms depends on the degree of the input polynomial, which is w(t) in our case, and its bit-size of coefficients. For instance, using the algorithm of [START_REF] Sagraloff | On the complexity of the Descartes method when using approximate arithmetic[END_REF] leads to a bit complexity

where τ is the largest bit-size of coefficients of w(t). While the degree of w(t) is already bounded by 2D(D -1) n-1 , τ is not estimated yet in this thesis. To identify a bound of τ, one needs to estimate the bit complexity of the subroutine Candidates [17, Subsection 3.2], especially the algorithm for computing at least one point per connected component of a real algebraic set given in [START_REF] El Din | Computing sampling points on a singular real hypersurface using Lagrange's system[END_REF]. This topic will be studied in future research.

Heuristic optimizations

Even though computing the constant e 0 requires at most D O(n) arithmetic operations in Q, its performance depends heavily on an efficient implementation of the geometric resolution algorithm over Q[t]/⟨w(t)⟩, which remains challenging to obtain. Thus, we aim to avoid such computations as much as possible. In what follows, we present two subroutines which are launched to test whether it is necessary for computing e 0 . In most of the case, with theses subroutines, our algorithm will return the set of isolated points without doing any further computation.

Heuristic identification of real isolated points

The optimization described in what follows identifies efficiently a subset of isolated points of H from the candidates without computing e 0 . The idea is to compute for each candidate x a ball B ∈ R n such that if the boundary of B intersects H, x is isolated in H.

We consider the polynomial f A obtained from f by a linear change of variable A and denote

Lemma 13. Let C be a bounded connected component of H A and 1 ≤ j ≤ n such that π j (C) is not a single point. Then, there exist at least two points in

and the projections π j (lim ε C i ) are connected. Thus, there exists 1 ≤ i ≤ k such that π j (lim ε C i ) is not a singleton. Hence, the projection of C i by π j is a closed interval [α, β] ⊂ R⟨ε⟩ where Algorithm 7: Implementation of IsolatedPoints Input:

Experimental results

In this section, we report on practical performances of our algorithms. Computations were done on an Intel(R) Xeon(R) CPU E7-4820 2GHz and 1.5 TB of RAM. Timings are given in seconds (s.), minutes (m.), hours (h.) and days (d.).

We take sums of squares of n random dense quadrics in n variables (with a non-empty intersection over R); we obtain dense quartics defining a finite set of points. None of the considered examples can be solved using CAD implementations in Maple within 10 days.

Table 1 below reports on the timings of our implementation (Algorithm 7). Timings for the subroutine Candidates are given in the column cand below. We use FGb library for computing Gröbner bases in order to perform algebraic elimination in our algorithms. The subroutine Has-RealSolutions in Algorithm 4 is done by raglib. Solving of zero-dimensional systems in the whole algorithm is done by msolve and real root isolation is done by the command RootFinding[Isolate] in maple.

The column cand2 shows the timings for computing the zero-dimensional parametrization C 2 and isolates its zeros (see Subsection 6.1). The column |real sols.|/ deg(w) shows the number of real candidates among the total number of candidates. This motivates the use of approximations, which runs only on candidates in R n , instead of computing over Q[t]/ ⟨w(t)⟩ which takes into account all candidates.

The column test1 reports on the timings of the first optimization (Algorithm 5). Exploiting the fact that isolating boxes are given by linear inequalities, we tweak raglib for solving the associated decision problems. As explained in the end of Subsection 6.1, by isolating C 2 with a small enough boxes in the subroutine SimpleIdentification, one can also obtain a certified output without computing e 0 . In our examples, it is the case and we do not need to carry out further computations. Timings of other steps are given as an indication for further researches.

Timings for ComputeE0 are given in the column e0. The columns approx and raglib respectively give the timings for computing the approximations and solving the decision problem by raglib. Note that the implementation used for two columns approx and raglib checks the emptiness of intersections of H with hypercubes (as explained in the end of Subsection 4.4). This