A note on the complexity of univariate root isolation

Keywords:

Introduction

One of the most important procedures in computer algebra and algebraic algorithms in general is root isolation of univariate polynomials. The goal of this algorithms is to compute isolating intervals in the real case, or boxes in the complex case that isolate the roots of the polynomial and to compute one such interval, or box, for every root.

We restrict ourselves to exact algorithms, i.e algorithms that perform arithmetic with rational numbers of arbitrary size.

The most well known and frequently used algorithms are the subdivision algorithms, either based on Sturm sequences, (sturm) or on Descartes' rule of sign (descartes) or on Descartes' rule of sign and the properties of Bernstein basis representation (bernstein). The subdivision algorithms mimic the process of the binary search algorithm. They consider an initial interval that contains all the real roots and the repeatedly subdivide it until is certified that zero or one root is contained in the tested interval.

Quite recently it was proven that the complexity of sturm [START_REF] Du | Amortized bound for root isolation via Sturm sequences[END_REF] and descartes and bernstein [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF] in the worst case is O B (d 4 τ ), where d is the degree of the polynomial and τ the maximum coefficient it size. Moreover, the same proof for the number of steps that the algorithms perform holds for all the solvers, the polynomial need not be square-free and in the same complexity bound we can also compute the multiplicities of the real roots [START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF].

Another type of exact real root isolation solver is the cf solver [START_REF] Akritas | An implementation of Vincent's theorem[END_REF], which we will not mention in this work and is based on the continued fractions expansion of the real roots. Quite recently it was proven [START_REF] Tsigaridas | Univariate polynomial real root isolation: Continued fractions revisited[END_REF], using the metric theory of the continued fractions of the real numbers that the expected complexity of cf is O B (d 4 τ 2 ), thus matches the current known bounds of the subdivision based algorithms. Moreover, by spreading the roots away the expected complexity bound can be improved to O B (d 3 τ ) [START_REF] Tsigaridas | On the complexity of real root isolation using Continued Fractions[END_REF]. This algorithm we will denote it by mcf.

The direct competitors of the exact algorithms are the numerical algorithms that compute an approximation, up to a desired accuracy, of all the complex roots of a polynomial. These algorithms can be turned to isolation algorithm by requiring the accuracy to be equal to the separation bound of the polynomial. The ones with the best complexity up to now, namely O B (d 3 τ ) [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF][START_REF] Pan | Univariate polynomials: Nearly optimal algorithms for numerical factorization and rootfinding[END_REF], are recursively splitting the polynomial until they compute linear factor that approximate sufficiently the real roots.

Even though the worst-case complexity bounds of the exact algorithms are worse than those of the numerical ones, recent implementations behave in practice in a rather satisfactory manner, e.g. [START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF][START_REF] Johnson | Highperformance implementations of the Descartes method[END_REF][START_REF] Tsigaridas | Univariate polynomial real root isolation: Continued fractions revisited[END_REF][START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF]. On the other hand, optimal numerical algorithms are very difficult to implement.

If we want to isolate the complex roots of a polynomial, i.e to compute isolating boxes in the complex plane that contain one and only one complex root, then the exact algorithms that are known are based on Sturm sequences [START_REF] Wilf | A global bisection algorithm for computing the zeros of polynomials in the complex plane[END_REF][START_REF] Pinkert | An exact method for finding the roots of a complex polynomial[END_REF].

Main results and overview

In this paper we study the average complexity subdivision solvers for real and complex root isolation. The complexity bounds that we prove improve the ones of the worst case analysis and are a step towards undestanding the behavior of the exact subdivision based root isolation algorithms and possibly to reduce the gap with the numerical ones.

Table 1 shows the complexities of the solvers. The results for the Continued Fractions (cf and mcf) solvers appear in [START_REF] Tsigaridas | Univariate polynomial real root isolation: Continued fractions revisited[END_REF][START_REF] Tsigaridas | On the complexity of real root isolation using Continued Fractions[END_REF]. All the other results on the expected complexity are new and to the best of our knowledge this is the first time that such results appear for the sturm and the descartes/bernstein solver. In particular, the former has output-sensitive complexity since it depends on the number of real roots.

Solver

Average case Worst case

sturm O B (rd 2 τ 2 ) O B (d 4 τ 2 ) descartes O B (d 3 τ 2 ) O B (d 4 τ 2 ) bernstein O B (d 3 τ 2 ) O B (d 4 τ 2 ) cf O B (d 4 τ 2 ) mcf O B (d 3 τ ) Numeric O B (d 3 τ )
Table 1: Complexity bounds for univariate real root isolation.

For the problem of complex root isolation we simplify significantly the proof [START_REF] Du | Amortized bound for root isolation via Sturm sequences[END_REF] for the number of subdivision that the algorithm perfrom and we also present an average case analysis.

The rest of the paper is structured as follows. The next subsection concentrates on notation. The following section describes the general procedure of a subdivision solver. In Sec. 3 we present the average bit size of the separation bound and we study the average complexity of the sturm and descartes algorithm. In Sec. 4 we study the problem of complex root isolation and in Sec. 5 presents Kronecker's algorithm for completeness, but also in order to settle an open question by Collins and Loos on its average-case behavior. Finaly, we present some open questions.

Notation

In what follows O B means bit complexity and the O B -notation means that we are ignoring logarithmic factors. For a polynomial A = d i=1 a i X i ∈ Z[X], deg(A) denotes its degree. We consider square-free polynomials except if explicitly stated otherwise. By L (A) we denote an upper bound on the bit size of the coefficients of A (including a bit for the sign). For a ∈ Q, L (a) ≥ 1 is the maximum bit size of the numerator and the denominator. Let M (τ ) denote the bit complexity of multiplying two integers of bit size at most τ . Using FFT, M (τ ) = O B (τ lg c τ ) for a suitable constant c. ∆ is the separation bound of A, that is the smallest distance between two (complex) roots of A. Finally N = max {d, τ }.

INRIA

Algorithm 1: subdivisionSolver(A, I 0 )

Input: A ∈ Z[X], I 0 = [a, b] Output: A list of isolatins intervals initialization sm (A, I 0 ) 1 L ← ∅ 2 Q ← ∅ 3 Q ← push(Q, {A, I 0 }) 4 while Q = ∅ do 5 {f, I} ← pop(Q) 6 V ← count sm (f, I) 7 switch V do 8 case V = 0 continue 9 case V = 1 L ← add(L, I) 10 case V > 1 11 {f L , I L }, {f R , I R } ← split sm (f, I) 12 Q ← push(Q, {f L , I L }) 13 Q ← push(Q, {f R , I R }) 14 return L 15
In what follows we consider a polynomial A ∈ Z[X] such that deg(A) = d and L (A) = τ . The polynomial is square-free except if explicitly stated otherwise.

Subdivision Solvers

In this section we will desrcibe the general subdivision algorithm for isolating the real roots of a univariate polynomial. Recall that the subdivision-based algorithms mimic the process of the binary search algorithm.

The pseudo-code of the general subdivision algorithm is presented in Alg. 1. The input is a square-free polynomial A ∈ Z[X] and an interval I 0 which contains the real roots of A which we wish to isolate. Usually I 0 contains all real roots of A. The algorithm uses a stack Q that contains pairs of the form {f, I}. The semantics of the pair is that we want to isolate the real roots of f contained in the interval I.

In the pseudo-code of Alg. 1 we also use some external functions. To be more specific push(Q, {f, I}) inserts the pair {f, I} to the top of the stack Q and pop(Q) returns the pair that is in top of the stack and deletes it from it. Moreover, add(L, I) inserts the interval I to the list L of the isolating intervals.

Finally, there are three sub-algorithms with index sm which have different specialization with respect to the subdivision algorithm applied. The first one initialization sm does all the necessary pre-processing, count sm (f, I) returns the number (or an upper bound) of the real roots of f in I. Finally, split sm (f, I) splits the interval I to two equal subintervals and possibly modifies polynomial f .

Notice that the complexity of the subdivision algorithm depends on how many times the while-loop (Line 5 of Alg. 1) is executed and on the cost of the sub-algorithms count sm (f, I) and split sm (f, I). Moreover, at every step of the algorithm, since the tested interval is splitted to two equal sub-intervals, we may assume that the bit size of the endpoints of the intervals is augmented by one bit at every step of the algorithm. To be more specific, if we assume that the endpoints of the initial interval I 0 have bit size τ , then at the h step of the algorithm the bit size of the endpoints of the tested intervals I ⊆ J 0 is τ + h.

The various subdivision algorithms differ by the way they implement the three subalgorithms. The most well known ones are sturm, descartes and bernstein.

On the average case complexity

Let disc(A) be the discriminant and lc (A) the leading coefficient of A. Mahler's measure of

a polynomial A is M(A) = | lc (A) | d i=1 max {1, |γ i |}, where γ i are all the (complex) roots of A [2, 33, 18, 19]. Moreover M(A) ≤ 2 τ √ d + 1.
We prove the following theorem, which is based on a theorem by Mignotte [START_REF] Mignotte | Mathematics for computer algebra[END_REF], see also [START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF].

Theorem 1 (Davenport-Mahler-Mignotte) Let A ∈ Z[X], with deg(A) = d and L (A) = τ , where A(0) = 0.
Let Ω be any set of k couples of indices (i, j) such that 1 ≤ i < j ≤ d and let the non-zero (complex) roots of A be 0

< |γ 1 | ≤ |γ 2 | ≤ • • • ≤ |γ d |. Then 2 k M(A) k ≥ (i,j)∈Ω |γ i -γ j | ≥ 2 k-d(d-1) 2 M(A) 1-d-k disc(A). Proof: Consider the multiset Ω = {j|(i, j) ∈ Ω}, |Ω| = k. We use the inequality ∀ a, b ∈ C |a -b| ≤ 2 max {|a|, |b|} (1)
and the fact [START_REF] Mignotte | Mathematics for computer algebra[END_REF][START_REF] Mignotte | Polynomials: An algorithmic approach[END_REF] that for any root of A, 1

M(A) ≤ |γ i | ≤ M(A).
In order to prove the left inequality

(i,j)∈Ω |γ i -γ j | ≤ 2 k j∈Ω |γ j | ≤ 2 k max j∈Ω |γ j | k ≤ 2 k M(A) k . INRIA Recall [33, 18] that disc(A) = lc (A) 2d-2 i<j (γ i -γ j ) 2 .
For the right inequality we consider the absolute value of the discriminant of A:

|disc(A)| = | lc (A) | 2d-2 i<j |γ i -γ j | 2 = | lc (A) | 2d-2 (i,j)∈Ω |γ i -γ j | 2 (i,j) / ∈Ω |γ i -γ j | 2 ⇔ |disc(A)| = | lc (A) | d-1 (i,j)∈Ω |γ i -γ j | (i,j) / ∈Ω |γ i -γ j | ⇔ (i,j)∈Ω |γ i -γ j | = | lc (A) | 1-d ( (i,j) / ∈Ω |γ i -γ j |) -1 |disc(A)| (2)
We consider the product (i,j) / ∈Ω |γ iγ j | and we apply d(d-1)

2

k times inequality (1), thus

(i,j) / ∈Ω |γ i -γ j | ≤ 2 d(d-1) 2 -k |γ 1 | 0 |γ 2 | 1 • • • |γ d | d-1 ( j∈Ω |γ j |) -1 .
From the definition of Mahler's bound we have

|γ 1 | 0 |γ 2 | 1 • • • |γ d | d-1 ≤ |M(A)/ lc (A) | d-1 ,
and since ∀i, |γ i | ≥ M(A) -1 j∈Ω |γ j | ≥ |γ 1 | k ≥ M(A) -k ,
we conclude that

(i,j) / ∈Ω |γ i -γ j | ≤ 2 d(d-1) 2 -k M(A) d+k-1 | lc (A) | 1-d .
Combining the previous inequality with the last equation of (2) we finish the proof of the theorem.

2

A similar theorem but with more strict hypotheses on the roots first appeared in [START_REF] Davenport | Cylindrical algebraic decomposition[END_REF] and the conditions were generalized in [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF], see also [START_REF] Johnson | Algorithms for Polynomial Real Root Isolation[END_REF][START_REF] Krandick | Isolierung reeller nullstellen von polynomen[END_REF]. Th. 1 has a factor 2 d 2 instead of d d in [START_REF] Davenport | Cylindrical algebraic decomposition[END_REF][START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF][START_REF] Johnson | Algorithms for Polynomial Real Root Isolation[END_REF], which plays no role when d = O(τ ) or when notation with N is used. Possibly a more involved proof may eliminate this factor using the techniques in [START_REF] Mignotte | On the Distance Between the Roots of a Polynomial[END_REF]. Moreover, the loose hypotheses of Th. 1 dramatically simplify the proof for the number of steps of the sudivision-based solvers [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF][START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF].

Recall that the separation bound is the minimum distance between two complex roots of a polynomial, i.e ∆ = min i =j |γ iγ j |, where γ i are the (complex) roots of the polynomial and 1 ≤ i, j ≤ d. Let γ ci be a complex root of A that is closest to γ i . We denote by ∆ i the quantity

∆ i = |γ i -γ ci |.
A separating point for γ i and γ ci is of magnitude at most 1 2 ∆ i . Our goal is to bound the average bit size of the separation points, i.e the average bit size of lg ∆ i , which will denote by A[lg ∆ i ].

Lemma 2

The average bit size of a separation point is O(d + τ ).

Proof: Our goal is to bound A[lg ∆ j ], where 1 ≤ j ≤ d and

A[lg ∆ i ] = 1 d d i=1 lg ∆ i = 1 d lg d i=1 ∆ i . (3)
First we consider the quantity d j=1 ∆ j . In order to apply Th. 1 we should rearrange d i=1 |γ iγ ci | so that the requirements on the indices of roots are fulfilled. This can not be achieved when symmetric factors occur, i.e the quantity contains factors of the form |(γ jγ cj )(γ cjγ j )|, for some indices j. In order to avoid this case, we consider the positive integers k 1 and k 2 such that k 1 + k 2 = d and we factorize the quantity as

d j=1 ∆ j = k1 j=1 ∆ j k2 j=1 ∆ j , (4)
where the factors do not contain symmetric products. Notice that since A is square-free disc(A) ≥ 1.

Let k ≤ d be such that the hypotheses of Th. 1, for the quantity d i=1 ∆ i , are fulfilled. Thus

2 k M (A) k ≥ k j=1 ∆ j ≥ 2 k-d(d-1) 2 M (A) 1-d-k 2 k+kτ d k ≥ k j=1 ∆ j ≥ 2 k-d(d-1) 2 2 τ (1-d-k) d 1-d-k k + kτ + k lg d ≥ k j=1 lg ∆ j ≥ k - d(d -1) 2 + τ (1 -d -k) + (1 -d -k) lg d.
If in the last relation we replace k once by k 1 and once by k 2 , then we obtain the relations:

k 1 + k 1 τ + k 1 lg d ≥ k1 j=1 lg ∆ j ≥ k 1 - d(d -1) 2 + τ (1 -d -k 1 ) + (1 -d -k 1 ) lg d, k 2 + k 2 τ + k 2 lg d ≥ k2 j=1 lg ∆ j ≥ k 2 - d(d -1) 2 + τ (1 -d -k 2 ) + (1 -d -k 2 ) lg d.
INRIA and if we sum them, taking into account that k 1 + k 2 = d, then we have

dτ + d(1 + lg d) ≥ d j=1 lg ∆ j ≥ -d 2 -3dτ + 2τ + 2d + (2 -3d) lg d.
If we consider the last equation and (3) then we conclude that A[lg ∆ j ] = O(d + τ ), under the mild assumption that lg d = O(τ ). 2

The Sturm solver

The Sturm subdivision solver sturm counts the number of distinct real roots of A in an interval as follows. It evaluates a signed polynomial sequence of A and its derivate A ′ over the left endpoint of the interval and counts the number of sign variations. It does the same for the right enpoint and the difference of the sign variations is the number of real roots.

Obviously the complexity of the algorithm is the number of steps that it performs time the complexity of each step.

Lemma 3

The average number of subdivision steps of sturm is O(r(d + τ )), where r is the number of real roots.

Proof: Recall that initially all the real roots are contained in the interval [-2 τ , 2 τ ]. Since A has r real roots we must compute r separation points. The magnitude of the separation points is 1 2 ∆ i , for 1 ≤ i ≤ r, and we can compute each of them with ⌈lg 2 τ +1 ∆j ⌉ subdivision steps, by performing binary search in the initial interval.

Let #(T ) be the total number of subdivisions that we need in order to isolate all the real roots. Or in other words the number of steps that we need in order to compute all the separation points. Then

#(T ) = r j=1 lg 2 τ +1 ∆ j ≤ 2r + rτ - r j=1 lg ∆ j ,
where the r summand is for the r possible roundings. By Lem. 3, we have

A[lg ∆ j ] = O(d + τ ), hence #(T ) = 2r + rτ + r j=1 O(d + τ ) = r(τ + 2) + O(d + τ ) r j=1 1,
and thus #(T ) = O(r(d + τ )) bounds the average number of subdivision steps. 2

A similar proof holds for the worst case analysis [START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF], but in this case r should be replaced by d.

We need the following theorem:

Theorem 4 [START_REF] Lickteig | Sylvester-Habicht Sequences and Fast Cauchy Index Computation[END_REF][START_REF] Reischert | Asymptotically fast computation of subresultants[END_REF] The evaluation of a signed polynomial remainder sequence of A and A ′ , SR(A, A ′ ), over a number a, where a ∈ Q ∪ {±∞} and L (a) = σ, has bit complexity O B (d 2 max{τ, σ}).

Theorem 5

The average complexity of sturm is O B (r(d

4 + d 3 τ + d 2 τ 2 )), or O B (rd 2 τ 2 ) if d = O(τ ), or O B (rN 4
) where N = max{d, τ }.

Proof: The complexity of the sturm solver is the number of subdivision steps that it performs times the cost of each subdivision. From Lem. 3, the average number of subdivisions is in O(r(d + τ )). Since the average bit size of the separation points is O B (d + τ ) the cost of evaluating the signed polynomial sequence over them is O B (d 3 + d 2 τ ) (Th. 4). We conclude that the overall complexity is O B (r(d

4 + d 3 τ + d 2 τ 2 ). 2
Remark 6 If we consider polynomials of degree d with coefficients that are independent standard normals then the expected number of real roots is r ∼ 2 π lg d [START_REF] Kac | On the average number of real roots of a random algebraic equation[END_REF], see also [START_REF] Edelman | How may zeros of a random polynomial are real[END_REF]. In this case the average complexity of sturm becomes O B (N 4 ), thus matching the complexity of numerical solvers [START_REF] Pan | Univariate polynomials: Nearly optimal algorithms for numerical factorization and rootfinding[END_REF][START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF].

As another, perhaps more natural, definition of random polynomials we can consider random polynomials of degree d with independent normally distributed coefficients, each with mean zero, but with the variace of the i-th coefficient equal to d i . In this case the expected number of real roots is r ∼ √ d [START_REF] Edelman | How may zeros of a random polynomial are real[END_REF] and the average complexity of sturm is O B (N 4.5 ).

The Descartes solver

The Descartes subdivision solver descartes relies on Descartes' rule of sign, which states that the number of sign variations in the coefficient list exceeds the number of positive real roots of A by an even number. In general this rule provides an overestimation on the number of positive real roots. However, when the number of sign variations is zero or one, then we get the exact number of positive real roots. Initially the polynomial is transformed so that its roots are contained in (0, 1). In order to count the number real roots of A in an interval, we transform the polynomial to another polynomial, the real roots of which in (0, ∞) correspond to the real roots of A in the initial interval. If the number of real roots is greater than one then we subdivide it and the we continue the algorithm on each subinterval. Notice that the polynomial transformations needed can be performed by shifting appropriately the variable X. The complexity of the algorithm is the number of intervals that it considers times the time needed for shifting the polynomial by a number whose magnitude is, in the worst case, proportional to the separation bound. Descartes' rule of sign is independent of the basis that the polynomial is represented. This allows a variant of descartes, which is based on the Bernstein representation of polynomials [START_REF] Lane | Bounds on a polynomial[END_REF][START_REF] Mourrain | On the complexity of isolating real roots and computing with certainty the topological degree[END_REF][START_REF] Mourrain | Bernstein's basis and real root isolation[END_REF] and is amenable to very efficient implementations [START_REF] Johnson | Highperformance implementations of the Descartes method[END_REF]. In what follows we will refer only to descartes solver but the results also hold for bernstein.

INRIA

Lemma 7

The average number of subdivision steps of descartes is O(d 2 + dτ ).

Proof: We consider the subdivision tree, T , that is formed during the execution of the algorithm. Each node of the tree corresponds to an interval that the algorithm tests for roots and a polynomial. The root of the tree corresponds to the initial polynomial and interval. The number of nodes of the tree is the number of steps that the algorithm performs. Since the algorithm terminates all the leaves of the tree correspond to polynomials with either on or zero sign variations.

We consider all the nodes of the tree that are parents of two leaves. The polynomials that correspond to these nodes have ≤ 2 sign variations. We denote the set of these nodes by J and notice that |J| ≤ d, since there are at most d complex roots and at after a subdivision step the number of sign variations can not increase.

If I ∈ J, then the number of subdivisions that the algorithm performs in order to compute it, is lg 2 τ +1

|I| . If we consider the paths from the root of the subdivision tree to all I ∈ J, then the sum of all the nodes bounds bounds the number of nodes of the the subdivision tree, #(T ). Thus

#(T ) ≤ I∈J lg 2 τ +1 |I| .
It can be proven [START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF] that each I ∈ J contains two, possible complex roots, such that

|γ i -γ j | ≤ 2|I|. But since |γ i -γ j | ≥ ∆ i we have that ∆ i ≤ 2|I|. If we also take into account that |J| ≤ d, then #(T ) ≤ I∈J lg 2 τ +1 |I| ≤ d i=1 lg 2 τ +2 ∆ i ≤ d(τ + 2) - d i=1 lg ∆ i . By Lem. 3, we have A[lg ∆ j ] = O(d + τ ), hence #(T ) ≤ d(d + τ ) + O(d + τ ) d i=1 1.
and thus the average number of subdivision steps is O(d 2 + dτ ). 2

An important observation concerning descartes is that the algorithm depends on all the roots of the polynomial, real and complex, thus we are not able to prove a bound on the number of the subdivision steps that wil depend on the number of the real roots.

We should also mention that for the worst case analysis the proof of Lem. 7 is almost the same. The only difference is that we have to bound the quantity d i=1 lg ∆ i using Th. 1 instead of Lem. 7. This approach simplifies even further the proof appeared in [START_REF] Emiris | Real Algebraic Numbers: Complexity Analysis and Experimentation[END_REF].

To complete the analysis we will need the following theorem:

Theorem 8 (Fast Taylor shift) [30] Let A ∈ Z[X], with deg(A) = d and L (A) = τ and let a ∈ Z, such that L (a) = σ. Then the cost of computing B = A(a + X) ∈ Z[X] is O B (M d 2 lg d + d 2 σ + dτ ). Moreover L (B) = O(τ + dσ).
Theorem 9 The average complexity of descartes is O B (d

5 + d 4 τ + d 3 τ 2 ), or O B (N 5 )
where N = max{d, τ }.

Proof: The complexity of the algorithm is the number of subdivision steps that performs times the cost of each subdivision. The number of subdivisions is O(d 2 + dτ ). The cost of each subdivision step is dominated by the cost of translating the polynomial. Notice that the bit size of the initial polynomial is O(dτ ) since we transformed so that its roots to be in (0, 1). In the average case we have to translate the polynomial by a number of bit size

O(d + τ ) (Lem. 2). The cost of this computation is O B (d 3 + d 2 τ ). Thus the overall cost, in average case is O B (d 5 + d 4 τ + d 3 τ 2 ). 2
4 Complex root isolation

Theorem 1 allows us to prove, with simple arguments, a bound on the number of subdivisions needed to isolate the complex roots of a real or Gaussian polynomial. The algorithm is based on Sturm sequences, and on the computation of the Cauchy index; for a complete description, the reader may refer to [START_REF] Pinkert | An exact method for finding the roots of a complex polynomial[END_REF][START_REF] Wilf | A global bisection algorithm for computing the zeros of polynomials in the complex plane[END_REF] and for a more didactic approach to [START_REF] Yap | Fundamental Problems of Algorithmic Algebra[END_REF].

The main idea of the algorithm is the following. Suppose that we have a method to count the correct number of distince complex roots in a square in the complex plane. If the number of complex roots is > 1 then the square is subdivided to four (equal) squares using the center point of the original square and the algorithm continues on each box.

Suppose that initialy all the complex roots are in a square of side B. This is a bound for all the complex roots, it can easily be computed [START_REF] Mignotte | Polynomials: An algorithmic approach[END_REF][START_REF] Mignotte | Mathematics for computer algebra[END_REF][START_REF] Yap | Fundamental Problems of Algorithmic Algebra[END_REF] and typically is B = 2 τ +1 . At the h step of the algorithm we have to check for complex roots in squares that have sides equal to B/2 h .

We can consider the process of the algorithm as a ternary tree, where each node holds a square and the root of the tree holds the intial one. Each leaf of the tree contains contains a square that isolates a complex root of the polynomial and since there are at most d complex roots, this is also the number of the leaves of the tree. The squares that correspond to the leaves of the tree have sides of length equal to ∆ j and the number of nodes from a leaf to the root of the tree is log R ∆ j .

The number of subdivisions, #(T ), that the algorithm perfroms equals the number of nodes of the subdivision tree, which is

#(T ) = d j=1 lg B ∆ j = 2d + dτ - d j=1 lg ∆ j = 2d + dτ -lg d j=1 ∆ j . (5)
It remains to bound the quantity d j=1 ∆ j . For this we use Th. 1. Recall that the hypotheses of the theorem are not fullfilled when symmetric products occur, that we factorize

INRIA quantity as d i=1 ∆ i = k1 i=1 ∆ i k2 i=1 ∆ i ,
where k 1 + k 2 = d and the factors are such that no symmetric products occur.

By applying to each factor Th. 1 and by taking into account that disc(A) ≥ 1 we have that

d i=1 ∆ i ≥ 2 d 2 M (A) 2-3d ≥ 2 d 2 (2 τ √ d + 1) 2-3d since M (A) ≤ 2 τ √ d + 1.
If we combine the last equation with [START_REF] Davenport | Cylindrical algebraic decomposition[END_REF] we conclude that the total number of subdivision is O(d 2 + dτ ). Thus we can state the following lemma:

Lemma 10 The number of subdivisions for complex root isolation is O(d 2 + dτ ).

The proof of the previous theorem simplifies significantly the proof appeared in [START_REF] Du | Amortized bound for root isolation via Sturm sequences[END_REF] where an amortized-like argument is used. Moreover, our bound on the number of subdivisions has a factor d 2 instead of a factor d lg d, which as in the real root case plays no role when the initial polynomial is not square-free or when the notation with N is used.

The only thing that remains is in order to complete the algorithm is to describe how we count the number or complex roots inside a box in the complex plane. For this we will follow Wilf [START_REF] Wilf | A global bisection algorithm for computing the zeros of polynomials in the complex plane[END_REF].

Let SR(f, g) denote a signed polynomial remainder sequence of the polynomials f and g and VAR(SR(f, g) ; a)) the number of sign variation occur when we evaluate the sequence over a rational number a. The Cauchy index of the real function g/f in an interval [a, b] is

I b a (g/f ) = VAR(SR(f, g) ; a)) -VAR(SR(f, g) ; b)).
Consider a box in the complex plane with vertices, q 1 , q 2 , q 3 , q 4 (q 5 = q 1 ) in counterclockwise order, centered around a point p. We expand the polynomial A about the point Q k , i.e we apply the transformation X → a k + i k-1 X, where i = √ -1. Let the resulting polynomial be

A k (X) = d j=o (a j -i b j )X j = d j=0 a i X j + i d j=0 b i X j . Moreover, let A k1 (X) = d j=0 a i X j , A k2 (X) = d j=0 b i X j .
Then the number of complex roots of A inside the box is

1 2 d i=1 (VAR(SR(A k1 , A k2 ) ; |q k+1 -q k |) -VAR(SR(A k1 , A k2 ) ; 0)).
In order to compute the overall complexity of the algorithm, we have to take into account that, at each step the algorithm evaluates a signed polynomial remainder sequence of polynomials of degree d and coefficient bit size O(dτ ). This is so because at each step we have to translate the polynomial. In the worst case we have to translate it by a number of bit size equal to the bit size of the separation bound, i.e O(dτ ), thus by Th. 8 its bit size will be O(d 2 τ ). The computation of the Cauchy index corresponds to the evaluation of a signed polynomial remainder sequence over a rational number, which has bit size at most O(dτ ). Using Th. 4 this evaluation costs O B (d 4 τ ). If we multiply the previous cost by the number of subdivisions we conclude that the overall cost is O B (d 6 τ + d 5 τ 2 ), or O B (N 7 ).

On the average complexity of complex root isolation

In order to average complexity of complex root isolation we will use the same techniques that we used in real root isolation. Recall that the average bit size of a separation point is

O(d + τ ) (Lem. 2).
As we did in the proof of Lem. 3 we assume that the number of complex roots is c and by direct computation we can easily prove that the number of subdivisions is O(c(d + τ )). Using again Lem.

Kronecker's algorithm

Kronecker [START_REF] Kronecker | Über den zahlbergriff[END_REF], see also [START_REF] Collins | Real zeros of polynomials[END_REF], presented an algorithm, hereafter called kronecker, for real root isolation of square-free polynomials. The algorithm depends on Rolle's theorem: if a square-free polynomial has a unique real root in an interval, then it changes sign at the endpoints of the interval. The algorithm is exponential but we discuss it for two reasons: First, in order to answer an open question posed by Collins and Loos [START_REF] Collins | Real zeros of polynomials[END_REF] about the expected number of intervals that the algorithm tests. Second, for completeness.

The input of the algorithm is a square-free polynomial A and its output is a list of the isolating intervals of the real roots. First the algorithm computes an absolute bound on the real roots. All bounds are the same asymptotically, so B = 2 τ , without loss of generality. Next the algorithm computes the theoretical separation bound ∆ = 2 -O(dτ ) . The main loop of the algorithm consists in splitting the initial interval [-B, B] to sub-intervals of length ∆ and to check, using Rolles' theorem, which of them contain a real root of A.

The complexity of the algorithm depends on the number of times its main loop is executed. Or in other words on the number of intervals that the algorithm tests. Since the initial interval is [-2 τ , 2 τ ] its length is 2 τ +1 . Thus the algorithm tests 2 O(dτ ) intervals. At every step kronecker evaluates A over a rational number of bit size at most O(dτ ). The complexity of the evaluation is O B (d 3 τ ) using Horner's scheme and fast multiplication algorithms [START_REF] Zur Gathen | Modern Computer Algebra[END_REF][START_REF] Bini | Polynomial and Matrix Computations[END_REF]. This leads to the following known theorem:

INRIA

There is an obvious improvement to kronecker algorithm. At each step, instead of computing the evaluation of A over the two endpoints of one interval, one can compute simultaneously the evaluation of A over the endpoints of d 2 intervals, with a logarithmic overhead on the asymptotic complexity [START_REF] Yap | Fundamental Problems of Algorithmic Algebra[END_REF][START_REF] Zur Gathen | Modern Computer Algebra[END_REF]. Unfortunately, this improvement does not tackle the exponential behavior of the algorithm, since we still have to test 2 O(dτ ) /d = 2 O(dτ ) intervals.

Collins and Loos [START_REF] Collins | Real zeros of polynomials[END_REF] posed the question about the expected number of intervals that kronecker tests. We will answer this question under the assumption that each interval may contain a real root of A with equal probability.

If we assume that the number of real roots, say r, it is not known and that r < d, then the algorithm, always, must test all the intervals since it is not possible to know when to stop.

In what follows we assume that we know in advance the exact number of real roots, r ≤ d. Let K be the number of the intervals that we have to test. All the possible positions of the real roots are K r . The number of intervals that we have to test depends on the position of the real root with the biggest magnitude. The dominant real root can be anywhere from position r to position K. If it is in position ℓ, where r ≤ ℓ ≤ K, then the algorithm tests ℓ intervals. There are ℓ-1 r such cases. Thus the expected value of the tested intervals is If r = 1, then obviously there is no need to run kronecker algorithm since we already have an isolating interval for the real root, i.e [-2 τ , 2 τ ]. However, if we run this naive algorithm anyway, for example because to we want to find a 2 -O(dτ ) approximation of the real root, then the expected number of intervals that we have to test is K+1 2 . This number coincides with the expected number of comparisons needed in order to locate an element in an array using serial search [START_REF] Graham | Concrete Mathematics[END_REF].

Open questions

The number of steps of the sturm solver depends on the number of real roots. If there are no real roots, then sturm identifies this without performing any subdivision. This is not the case for descartes: the number of steps that it performs depends on all complex roots of the polynomial. A relevant open question is whether there exists a class of polynomials with no real roots, such that the number of steps of descartes is O(dτ ) or even O(τ )?

Is there a exact algorithm for real root isolation with complexity O B (N 4 ) or even O B (N 5 )?

  3 the cost of each step is O B (d 4 + dτ ) and thus the overall cost in the average case is O B (c(d 5 + d 4 τ + d 3 τ 2 )), or O B (cN 5 ).

Theorem 11 (

 11 Kronecker) [4] Let a square-free polynomial A red ∈ Z[X] such that deg(A red ) = d and L (A red ) = τ . Algorithm kronecker isolates all real roots of A red in O B (2 dτ d 3 τ ).

RR n 6043


Danny Munera
email: danny.munera@malix.univ-paris1.fr

Solving Stable

Solving Stable Matching Problems via Cooperative Parallel Local Search

Keywords: stable matching, local search, adaptive search, parallelism

Stable matching problems and its variants have several practical applications, like the Hospital/Residents problem, stable roommates problem or bipartite market sharing. An important generalization problem is the SMTI which allows for incompleteness and ties in the user's preference lists. Finding a maximal size stable matching for SMTI is computationally difficult. We developed a Local Search method to solve SMTI using the Adaptive Search algorithm and present experimental evidence that this approach is much more efficient than state-of-the-art exact and approximate methods (in terms of both computational effort required and quality of solution). We also tried a parallel version of our algorithm. For this we reused the Cooperative Parallel Local Search framework (CPLS) we designed. CPLS is a highly parametric framework for the execution in parallel of local search solvers allowing them to cooperate though communication. The cooperative parallel version of our local search algorithm improves performance so much that very large and hard instances can be solved quickly.

Introduction

Stable Marriage (SM) problem [START_REF] Gale | College admissions and the stability of marriage[END_REF] is the problem to find a stable matching between a set of n men and a set of n women, each of whom have ranked all members of the other set in a strict order of preference. This match is stable if there is no man-woman pair where both would rather marry each other than their current partner -such a pair is called a blocking pair. A generalization of this problem is the Stable Marriage with Incomplete List and Ties (SMTI) [START_REF] Manlove | Hard variants of stable marriage[END_REF]. In SMTI, the preference lists may include ties (to express indifference) and may be incomplete (to express that some partners are unacceptable). The goal now is to find the stable matching of maximal size (that is, with the smallest number of singles). This problem is NP-hard [START_REF] Iwama | Stable Marriage with Incomplete Lists and Ties[END_REF].

Stable matching problems have many practical applications such as car-sharing or bipartite market sharing, job markets and social networks. Another important application is the Hospitals/Residents problem (HR) [START_REF] Manlove | Hard variants of stable marriage[END_REF] which consists in assigning resident doctors to hospitals, based on their preferences. For HR, there are nation-wide programmes in various countries like the Scottish Foundation Allocation Scheme (SFAS), the Canadian Resident Matching Service (CARMS) or the the National Resident Matching Program (NRMP) in the USA. Many of these applications involve very large sets and it is thus a real challenge to design efficient algorithms (in terms of execution time and/or solution quality).

We address the SMTI problem using Local Search (LS) and parallelism. For this, we first proposed a Cooperative Parallel Local Search framework (CPLS) [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF] which allows the user to run several Local Search solvers in parallel with sophisticated forms of cooperation (through communication). This framework is highly parametric and can be customized for different situations (e.g. the user can modify the LS algorithm to use, the level of "intensification" and "diversification" in the search, the communication topology). Due to space limitation, we do not describe more the CPLS framework, the interested reader may refer to [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF].

We then proposed a Local Search algorithm based on the Adaptive Search (AS) method [START_REF] Munera | Solving Hard Stable Matching Problems via Local Search and Cooperative Parallelization[END_REF]. For this, we showed how to model SMTI problems as permutation problems to limit the neighbourhood to explore. Furthermore, we designed efficient heuristics based on the properties of the blocking pairs to improve the current assignment. Other heuristics are also proposed to escape local minima. The resulting algorithm (called AS-SMTI) can be implemented efficiently with a compact data representation. We present experimental evidence that this sequential algorithm is much more efficient than state-of-the-art exact and approximate methods. Moreover, we developed a parallel implementation of the AS-SMTI using the CPLS framework (written in the X10 language [START_REF] Charles | X10 : an object-oriented approach to non-uniform cluster computing[END_REF]). We show that the cooperative parallel version exhibits super-linear speedup on average and behaves particularly well on hard instances.

In AS-SMTI, we model SMTI as a permutation problem : the sequence of n (X 1 . . . X n ) takes on as values permutations of the values 1 . . . n (implementing an all-different constraint). X i = j is interpreted as either (m i , w j ) in a match M , or m i is single if w j is not on its preference list. The AS-SMTI algorithm starts from a random matching, then it explores a limited neighbourhood based on a heuristic which selects the "worst" blocking-pair (BP) to fix and/or a single man to marry. If a pair (m, w) forms a BP in M , the error associated to this BP is the distance between the woman w and the current partner of m (X m = w ) in the preference list of m. Thus, the further the assigned woman w is from w, the larger the error. To fix the "worst" BP, AS-SMTI swaps the values of the BP variable X m = w and the variable that contains the desired match, X i = w (i.e. two men exchange their partner).

The total cost function of a matching M measures both its stability (number of BPs) and its quality (number of singles). Hence : cost(M ) = #BP (M ) × n + #Singles(M ) where #BP (M ) is the number of BPs in M , and #Singles(M ) is the number of singles in M . The number of BPs is weighted with n to prioritize stable marriages over marriages with fewer singles.

Finally if a local minimum is reached, AS-SMTI executes a customized reset procedure which basically tries to fix the 2 worst BPs and/or to assign a woman to a single man. This procedure is stochastic ; it will also fix the second worst variable with a probability p : good results are obtained with a high value, e.g. p 0.98.

The algorithm stops when a perfect solution (PS) is found (a stable matching with no singles) or when a given timeout is reached (the best matching found so far is then returned).

Performance Evaluation

We compared our AS-SMTI algorithm to other approaches both from the point of view of solution quality and of pure performance. For this comparison we selected three different methods : another Local Search method (the LTIU algorithm of [START_REF] Gelain | Local Search Approaches in Stable Matching Problems[END_REF]), a very efficient 3/2-approximation algorithm (the McDermid's method [START_REF] Mcdermid | A 3/2 -approximation algorithm for general stable marriage[END_REF]) and a SAT approach (using the encoding proposed in [START_REF] Gent | SAT encodings of the stable marriage problem with ties and incomplete lists[END_REF]). For this evaluation, we used the random problem generator described in [START_REF] Gent | An Empirical Study of the Stable Marriage Problem with Ties and Incomplete Lists[END_REF] which takes three parameters : the size (n), the probability of incompleteness (p1) and the probability of ties (p2). We generated problems of size n = 100, with p1 ranging over [0.1, 0.9] and p2 over [0, 1], with step 0.1. We used an X10 [START_REF] Charles | X10 : an object-oriented approach to non-uniform cluster computing[END_REF] implementation of AS-SMTI, running sequentially on an AMD Opteron 6376 clocked at 2.3 GHz, i.e. using only one core.

Regarding solution quality, we measured the percentage of perfect solutions (PS) reached by all algorithms. In most of the cases, the solutions returned by AS-SMTI are better than solutions returned by its opponents. On average, AS-SMTI is slightly better than LTIU. However, LTIU reports a dramatic lost of performance when p2 = 1 obtaining less than 60% of PS while AS-SMTI reaches 100%. For McDermid (MD), the percentage of PS found by AS-SMTI is considerably higher than this obtained by MD, in particular using a probability of ties p2 ∈ [0.1..0.7]. The reference SAT encoding was restricted to the decision problem : is there a stable matching of size n ? which we answer by actually finding a perfect stable matching. AS-SMTI always found a PS for the tested cases.

Regarding the performance (which measures the computational effort), Figure 1 compares AS-SMTI with its opponents. AS-SMTI is much faster that LTIU and SAT (thus the log Y scale) : the LTIU method averages over 30s on a similar machine while AS-SMTI is two orders of magnitude faster (when p2 increases, the lead extends even further). Similarly, AS-SMTI outperforms the SAT version by a factor of about 50. From the raw performance point of view, the comparison with McDermid is particularly significant. Recall that an approximation algorithm gives priority to the speed (running in linear time) at the expense of the quality. Surprisingly, in many cases, AS-SMTI is up to an order of magnitude faster than MD. However, when the probability of incompleteness is very high (e.g. p1 = 0.9) MD outperforms AS-SMTI. It is worth noticing that MD always returns the same, single and (sub)optimal solution, while AS-SMTI will yield more than one solution, with observably better quality. Moreover, a solution quality vs. performance trade-off is always possible in AS-SMTI, by tweaking the timeout parameter.

Parallel Evaluation

We also experimented with a parallel version of AS-SMTI. For this, we used the Cooperative Parallel Local Search (CPLS) framework we previously proposed [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]. Due to space limitation we do not detail how we instantiated and tuned CPLS for the SMTI case (see [START_REF] Munera | Solving Hard Stable Matching Problems via Local Search and Cooperative Parallelization[END_REF]). We attacked problems of size 1000 generated with p1 = 0.95 and p2 = 0.8. These problems involve a large number of variables, a huge search space (1000! 10 2567 ) and are difficult to solve due to the high level of incompleteness in the preference lists. All selected problem admit a perfect solution. We generated two different sets of problems. The first set, called normal, is composed by 10 random SMTI problems of size n = 1000. The second set, called hard, was generated from 100 random problems of size n = 1000, we ran them sequentially and selected the 10 hardest instances.

Each problem was executed 50 times (the results are averaged), varying the number of cores from 1 (sequential) to 128 and using an unlimited timeout to force the solver to discover the perfect solution. The used architecture is a cluster of 4 machines, each with 4 × 16-core AMD Opteron 6376 CPUs running at 2.3 GHz and 128 GB of RAM. Figure 2 presents a log-log graph of execution times. The Ideal time corresponds to a linear speedup : time is halved when the number of cores is doubled. On the normal data set, our parallel AS-SMTI is rather efficient : using 80 cores the average time decreases from 44.5s to 0.519s which corresponds to a speedup factor of 86. Beyond this number of cores, the time tends to stabilize. This could mean that we have reached an incompressible limit to solve problems of size 1000, or that we are being limited by low-level factors.

Regarding the hard problems, the difficulty of the instances is clear from the sequential time : in the normal experiment, a problem was solved in about 44s while about 285s are needed for the hardest instances. Using the cooperative search, execution times are drastically improved and the best speedup is 492 with 128 cores which corresponds to a reduction of the execution time from 284.5s to 0.579s. It is worth noticing that this time is very similar to the best time (0.519s) obtained on the normal set. From a practical point of view, it appears that parallelism with cooperation neutralizes the relative difficulty of problems, as instances which are originally about 6 times harder get solved in approximately the same time. Of course, the problem retains its worst-case NP-hard complexity, and parallel search cannot change this.

We proposed the AS-SMTI algorithm which models SMTIs as permutation problems and solves them using a local search approach based on Adaptive Search. AS-SMTI outperforms state-of-the-art solvers for SMTIs. Moreover, the cooperative parallel version shows superlinear speedup and performs exceptionally well, particularly on very hard instances. We plan to experiment with very large instances on a massively parallel machine to test the scaling limitations. We also plan to study which features of the stable matching problem make it so suitable for cooperation.

Currently we are attacking the Hospital/Resident problem with Ties (HRT). We have already developed a solution based on model transformation : an HRT problem is transformed into SMTI which is handled by our AS-SMTI solver (slightly modified). The preliminary results on a sequential machine are very encouraging : execution times are much faster than the best Integer Programming approach while obtained solution are very close to the optimum 1 . We are now experimenting with parallelism on HRT. We also plan to improve the core solver to optimize features appearing when transforming HRT to SMTI (such optimizations will also benefit to other SMTI problems).

FIG. 1 -

 1 FIG.1-Execution time comparing AS-SMTI to LTIU, MD and SAT (times in s)

FIG. 2 -

 2 FIG. 2 -Execution time for normal and hard problems varying the number of cores.

Acknowledgments

I would like to thank all the co-authors of the papers included on this work for their contributions, help, and encouragements throughout this project.