Peggy Cénac

Brigitte Chauvin

Frédéric Paccaut

Nicolas Pouyanne

Cénac, Chauvin, Paccaut, Pouyanne: uncommon suffix tries Uncommon Suffix Tries

Keywords: MSC 2010 60J05, 37E05 variable length Markov chain, probabilistic source, mixing properties, suffix trie

Introduction

Trie (abbreviation of retrieval) is a natural data structure, efficient for searching words in a given set and used in many algorithms as data compression, spell checking or IP addresses lookup. A trie is a digital tree in which words are inserted in external nodes. The trie process grows up by successively inserting words according to their prefixes. A precise definition will be given in Section 4.1. As soon as a set of words is given, the way they are inserted in the trie is deterministic. Nevertheless, a trie becomes random when the words are randomly drawn: each word is produced by a probabilistic source and n words are chosen (usually independently) to be inserted in a trie. A suffix trie is a trie built on the suffixes of one infinite word. The randomness then comes from the source producing such an infinite word and the successive words inserted in the tree are far from being independent, they are strongly correlated. As a principal application of suffix tries one can cite the lossless compression algorithm Lempel-Ziv 77 (LZ77). The first results on the average size of suffix tries when the infinite word is given by a symmetrical memoryless source are due to Blumer et al. [START_REF] Blumer | Average sizes of suffix trees and dawgs[END_REF] and those on the height of the tree to Devroye [START_REF] Devroye | A note on the height of suffix trees[END_REF]. Using analytic combinatorics, Fayolle [START_REF] Fayolle | Compression de données sans perte et combinatoire analytique[END_REF] has obtained the average size and total path length of the tree for a binary word issued from a memoryless source (with some restriction on the probability of each letter).

Here we are interested in the height H n and the saturation level ℓ n of a suffix trie T n containing the first n suffixes of an infinite word produced by a source associated with a so-called Variable Length Markov Chain (VLMC) (see Rissanen [START_REF] Rissanen | A universal data compression system[END_REF] for the seminal work, Galves-Löcherbach [START_REF] Galves | Stochastic chains with memory of variable length[END_REF] for an overview, and [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] for a probabilistic frame). One deals with a particular VLMC source associated with an infinite comb, described hereafter. This particular model has the double advantage to go beyond the cases of memoryless or Markov sources and to provide concrete computable properties. The analysis of the height and the saturation level is usually motivated by optimization of the memory cost. Height is clearly relevant to this point; saturation level is algorithmically relevant as well because internal nodes below the saturation level are often replaced by a less expansive table .  All the tries or suffix tries considered so far in the literature have a height and a saturation level both growing logarithmically with the number of words inserted, to the best of our knowledge. For plain tries, when the inserted words are independent, the results due to Pittel [START_REF] Pittel | Asymptotic growth of a class of random trees[END_REF] rely on two assumptions on the source producing the words: first, the source is uniformly mixing, second, the probability of any word decays exponentially with its length. Let us also mention the general analysis of tries by Clément-Flajolet-Vallée [START_REF] Clément | Dynamical sources in information theory: a general analysis of trie structures[END_REF] for dynamical sources. For suffix tries, Szpankowski [START_REF] Szpankowski | Asymptotic properties of data compression and suffix trees[END_REF] obtains the same result, with a weaker mixing assumption (still uniform though) and the same hypothesis on the measure of the words. Our aim is to exhibit two cases when these behaviours are no longer the same. The first example is the "logarithmic comb", for which we show that the mixing is slow in some sense, namely non uniformly polynomial (see Section 3.2 for a precise statement) and the measure of some increasing sequence of words decays polynomially. We prove in Theorem 4.8 that the height of this trie is larger than a power of n (when n is the number of inserted suffixes in the tree). The second example is the "factorial comb", which has a uniformly exponential mixing, thus fulfilling the mixing hypothesis of Szpankowski [START_REF] Szpankowski | Asymptotic properties of data compression and suffix trees[END_REF], but the measure of some increasing sequence of words decays faster than any exponential. In this case we prove in Theorem 4.9 that the saturation level is negligible with respect to log n. We prove more precisely that, almost surely,

ℓ n ∈ o log n (log log n) δ , for any δ > 1.
The paper is organised as follows. In Section 2, we define a VLMC source associated with an infinite comb. In Section, 3 we give results on the mixing properties of these sources by explicitely computing the suitable generating functions in terms of the source data. In Section 4, the associated suffix tries are built, and the two uncommon behaviours are stated and shown. The methods are based on two key tools concerning pattern return time: a duality property and the computation of generating functions. The relation between the mixing of the source and the asymptotic behaviour of the trie is highlighted by the proof of Proposition 4.7.

Infinite combs as sources

In this section, a VLMC probabilistic source associated with an infinite comb is defined. Moreover, we introduce the two examples given in introduction: the logarithmic and the factorial combs. We begin with the definition of a general variable length Markov Chain associated with a probabilized infinite comb. The following presentation comes from [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF]. Let A be the alphabet {0, 1} and L = A -N be the set of left-infinite words. Consider the binary tree (represented in Figure 1) whose finite leaves are the words 1, 01, . . . , 0 k 1, . . . and with an infinite leaf 0 ∞ as well. Each leaf is labelled with a Bernoulli distribution, respectively denoted by q 0 k 1 , k 0 and q 0 ∞ . This probabilized tree is called the infinite comb. The VLMC (Variable Length Markov Chain) associated with an infinite comb is the L-valued Markov chain (V n ) n 0 defined by the transitions

P(V n+1 = V n α|V n ) = q ←- pref (Vn) (α)
where α ∈ A is any letter and ←pref (V n ) denotes the first suffix of V n (reading from right to left) appearing as a leaf of the infinite comb. For instance, if V n = . . . 1000, then ←pref (V n ) = 0001. Notice that the VLMC is entirely determined by the data Figure 1: An infinite comb q 0 ∞ , q 0 k 1 , k 0. From now on, denote c 0 = 1 and for n 1,

c n := n-1 k=0 q 0 k 1 (0).
It is proved in [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] that in the irreducible case i.e. when q 0 ∞ (0) = 1, there exists a unique stationary probability measure π on L for (V n ) n if and only if the series c n converges. From now on, we assume that this condition is fulfilled and we call S(x) :=

n 0 c n x n (1)
its generating function so that S(1) = n 0 c n . For any finite word w, we denote π(w) := π(Lw). Computations performed in [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] show that for any n 0,

π(10 n ) = c n S(1)
and π(0

n ) = k n c k S(1) . (2)
Notice that, by stationarity π(0 n ) = π(0 n+1 ) + π(10 n ) and by disjointness of events, π(0 n ) = π(0 n+1 ) + π(0 n 1) for all n 1 so that π(10 n ) = π(0 n 1).

(

) 3
If U n denotes the final letter of V n , the random sequence W = U 0 U 1 U 2 . . . is a right-infinite random word. We define in this way a probabilistic source in the sense of information theory i.e. a mechanism that produces random words. This VLMC probabilistic source is characterized by: p w := P(W has w as a prefix ) = π(w),

for every finite word w. Both particular suffix tries the article deals with are built from such sources, defined by the following data.

Example 1: the logarithmic comb

The logarithmic comb is defined by c 0 = 1 and for n 1,

c n = 1 n(n + 1)(n + 2)(n + 3)
.

The corresponding conditional probabilities on the leaves of the tree are

q 1 (0) = 1 24 and for n 1, q 0 n 1 (0) = 1 - 4 n + 4
.

The expression of c n was chosen to make the computations as simple as possible and also because the square-integrability of the waiting time of some pattern will be needed (see end of Section 4.3), guaranteed by

n 0 n 2 c n < +∞.
Example 2: the factorial comb

The conditional probabilities on the leaves are defined by

q 0 n 1 (0) = 1 n + 2
for n 0, so that c n = 1 (n + 1)! .

Mixing properties of infinite combs

In this section, we first precise what we mean by mixing properties of a random sequence. We refer to Doukhan [START_REF] Doukhan | Mixing : properties and examples[END_REF], especially for the notion of ψ-mixing defined in that book. We state in Proposition 3.2 a general result that provides the mixing coefficient for an infinite comb defined by (c n ) n 0 or equivalently by its generating function S. This result is then applied to our two examples. The mixing of the logarithmic comb is polynomial but not uniform, it is a very weak mixing; the mixing of the factorial comb is uniform and exponential, it is a very strong mixing. Notice that mixing properties of some infinite combs have already been investigated by Isola [START_REF] Isola | Renewal sequences and intermittency[END_REF], although with a slight different language.

Mixing properties of general infinite combs

For a stationary sequence (U n ) n 0 with stationary measure π, we want to measure by means of a suitable coefficient the independence between two words A and B separated by n letters. The sequence is said to be "mixing" when this coefficient vanishes when n goes to +∞. Among all types of mixing, we focus on one of the strongest type: ψ-mixing. More precisely, for 0 m +∞, denote by F m 0 the σ-algebra generated by {U k , 0 k m} and introduce for A ∈ F m 0 and B ∈ F ∞ 0 the mixing coefficient

ψ(n, A, B) := π(A ∩ T -(m+1)-n B) -π(A)π(B) π(A)π(B) = |w|=n π(AwB) -π(A)π(B) π(A)π(B) , ( 4
)
where T is the shift map and where the sum runs over the finite words w with length |w| = n.

A sequence (U n ) n 0 is called ψ-mixing whenever lim n→∞ sup m 0,A∈F m 0 ,B∈F ∞ 0 |ψ(n, A, B)| = 0.
In this definition, the convergence to zero is uniform over all words A and B. This is not going to be the case in our first example. As in Isola [START_REF] Isola | Renewal sequences and intermittency[END_REF], we widely use the renewal properties of infinite combs (see Lemma 3.1) but more detailed results are needed, in particular we investigate the lack of uniformity for the logarithmic comb.

Notations and Generating functions

• For a comb, recall that S is the generating function of the nonincreasing sequence (c n ) n 0 defined by (1).

• Set ρ 0 = 0 and for n 1,

ρ n := c n-1 -c n ,
with generating function

P (x) := n 1 ρ n x n .
• Define the sequence (u n ) n 0 by u 0 = 1 and for n 1,

u n := π(U 0 = 1, U n = 1) π(1) = 1 π(1) |w|=n-1 π(1w1), (5)
and let

U(x) := n 0 u n x n
denote its generating function. Hereunder is stated a key lemma that will be widely used in Proposition 3.2. In some sense, this kind of relation (sometimes called Renewal Equation) reflects the renewal properties of the infinite comb.

Lemma 3.1 The sequences (u n ) n 0 and (ρ n ) n 0 are connected by the relations:

∀n 1, u n = ρ n + u 1 ρ n-1 + . . . + u n-1 ρ 1
and (consequently)

U(x) = ∞ n=0 u n x n = 1 1 -P (x) = 1 (1 -x)S(x) .
Proof. For a finite word w = α 1 . . . α m such that w = 0 m , let l(w) denote the position of the last 1 in w, that is l(w) := max{1 ≤ i ≤ m, α i = 1}. Then, the sum in the expression (5) of u n can be decomposed as follows:

|w|=n-1 π(1w1) = π(10 n-1 1) + n-1 i=1 |w|=n-1 l(w)=i π(1w1).
Now, by disjoint union π(10 n-1 ) = π(10 n-1 1) + π(10 n ), so that

π(10 n-1 1) = π(1)(c n-1 -c n ) = π(1)ρ n .
In the same way, for w

= α 1 . . . α n-1 , if l(w) = i then π(1w1) = π(1α 1 . . . α i-1 1)ρ n-i , so that u n = ρ n + n-1 i=1 ρ n-i 1 π(1) |w|=i-1 π(1w1) = ρ n + n-1 i=1 ρ n-i u i , which leads to U(x) = (1 -P (x)) -1 by summation.
⊓ ⊔

Mixing coefficients

The mixing coefficients ψ(n, A, B) are expressed as the n-th coefficient in the series expansion of an analytic function M A,B which is given in terms of S and U.

The notation [x n ]A(x) means the coefficient of x n in the power expansion of A(x) at the origin. Denote the remainders associated with the series S(x) by

r n := k n c k , R n (x) := k n c k x k
and for a 0, define the "shifted" generating function

P a (x) := 1 c a n 1 ρ a+n x n = x + x -1 c a x a R a+1 (x). ( 6
)
Proposition 3.2 For any finite word A and any word B, the identity

ψ(n, A, B) = [x n+1 ]M A,B (x)
holds for the generating functions M A,B respectively defined by:

i) if A = A ′ 1 and B = 1B ′
where A ′ and B ′ are any finite words, then

M A,B (x) = M(x) := S(x) -S(1) (x -1)S(x) ; ii) if A = A ′
M A,B (x) := S(1) c a+b c a c b P a+b (x) + U(x) [S(1)P a (x)P b (x) -S(x)] ; iii) if A = 0 a and B = 0 b with a, b 1, then M A,B (x) := S(1) 1 r a r b n 1 r a+b+n-1 x n + U(x) S(1)R a (x)R b (x) r a r b x a+b-2 -S(x) ; iv) if A = A ′
M A,B (x) := S(1) 1 c a r b x a+b-1 R a+b (x) + U(x) S(1)P a (x)R b (x) c a r b x b-1 -S(x) ; v) if A = 0 a and B = 0 b 1B ′
where B ′ is any finite words and a, b 0, then

M A,B (x) := S(1) 1 r a c b x a+b-1 R a+b (x) + U(x) S(1)R a (x)P b (x) r a c b x a-1 -S(x) .
Remark 3.3 It is worth noticing that the asymptotics of ψ(n, A, B) may not be uniform in all words A and B. We call this kind of system non-uniformly ψmixing. It may happen that ψ(n, A, B) goes to zero for any fixed A and B, but (for example, in case iii)) the larger a or b, the slower the convergence, preventing it from being uniform.

Proof. The following identity has been established in [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] (see formula (17) in that paper) and will be used many times in the sequel. For any two finite words w and w ′ , π(w1w ′ )π(1) = π(w1)π(1w ′ ). ( 7)

i) If A = A ′ 1 and B = 1B ′ , then (7) yields π(AwB) = π(A ′ 1w1B ′ ) = π(A ′ 1) π(1) π(1w1B ′ ) = S(1)π(A)π(B) π(1w1) π(1) . So ψ(n, A, B) = S(1)u n+1 -1
and by Lemma 3.1, the result follows.

ii) Let A = A ′ 10 a and B = 0 b 1B ′ with a, b 0 and a + b = 0. To begin with,

π(AwB) = 1 π(1) π(A ′ 1)π(10 a w0 b 1B ′ ) = 1 π(1) 2 π(A ′ 1)π(10 a w0 b 1)π(1B ′ ).
Furthermore, π(A) = c a π(A ′ 1) and by (3), π(0 b 1) = π(10 b ), so it comes

π(B) = 1 π(1) π(0 b 1)π(1B ′ ) = π(10 b ) π(1) π(1B ′ ) = c b π(1B ′ ).
Therefore,

π(AwB) = π(A)π(B) c a c b π(1) 2 π(10 a w0 b 1). Using π(1)S(1) = 1, this proves ψ(n, A, B) = S(1) v a,b n c a c b -1 where v a,b n := 1 π(1) |w|=n-1 π(10 a w0 b 1).
As in the proof of the previous lemma, if w = α 1 . . . α m is any finite word different from 0 m , we call f (w) := min{1 ≤ i ≤ m, α i = 1} the first place where 1 can be seen in w and recall that l(w) denotes the last place where 1 can be seen in w. One has If i = j then w is the word 0 i-1 10 n-i-1 , else w is of the form 0 i-1 1w ′ 10 n-1-j , with |w ′ | = ji -1. Hence, the previous sum can be rewritten as

|w|=n-1 π(10 a w0 b 1) = π(1)ρ a+b+n + π(1) n-1 i=1 ρ a+i ρ n-i+b + 1 i<j n-1 w |w|=j-i-1 π(10 a+i-1 1w10 n-1-j+b 1).
Equation [START_REF] Flajolet | Analytic Combinatorics[END_REF] shows

π(10 a+i-1 1w10 n-1-j+b 1) = π(10 a+i-1 1) π(1) π(1w1) π(1) π(10 n-1-j+b 1) = ρ a+i ρ n-j+b π(1w1).
This implies:

v a,b n = ρ a+b+n + n-1 i=1 ρ a+i ρ n-i+b + 1≤i<j≤n-1 ρ a+i ρ n-j+b w,|w|=j-i+1 π(1w1) π (1)
.

Recalling that u 0 = 1, one gets

v a,b n = ρ a+b+n + 1≤i≤j≤n-1 ρ a+i ρ n-j+b u j-i
which gives the result ii) with Lemma 3.1.

iii)

Let A = 0 a and B = 0 b with a, b 1. Set v a,b n := 1 π(1) |w|=n-1 π(0 a w0 b ). First, recall that, due to (2), π(A) = π(1)r a and π(B) = π(1)r b . Conse- quently, ψ(n, A, B) = π(1)v a,b n+1 -π(A)π(B) π(A)π(B) = S(1) v a,b n+1 r a r b -1.
Let w be a finite word with

|w| = n -1. If w = 0 n-1 , then π(AwB) = π(0 a+n-1+b ) = π(1)r a+b+n-1 .
If not, let f (w) denote as before the first position of 1 in w and l(w) the last one in w. If f (w) = l(w), then

π(AwB) = π(0 a+f (w)-1 10 n-1-f (w)+b ) = 1 π(1) π(0 a+f (w)-1 1)π(10 n-1-f (w)+b ) = π(1)c a+f (w)-1 c n-1-f (w)+b . If f (w) < l(w), then writing w = w 1 . . . w n-1 , π(AwB) = π(0 a+f (w)-1 1w f (w)+1 . . . w l(w)-1 10 n-1-l(w) ) = 1 π(1) 2 π(0 a+f (w)-1 1)π(1w f (w)+1 . . . w l(w)-1 1)π(10 n-1-l(w)+b ). Summing yields v a,b n = r a+b+n-1 + n-1 i=1 c a+i-1 c n-1+b-i + n-1 i,j=1 i<j w, |w|=j-i-1 c a+i-1 π(1w1) π(1) c n-1+b-j = r a+b+n-1 + 1≤i≤j≤n-1 c a+i-1 c n-1+b-j u j-i ,
which gives the desired result. The last two items, left to the reader, follow the same guidelines. ⊓ ⊔

Mixing of the logarithmic infinite comb

Consider the first example in Section 2, that is the probabilized infinite comb defined by c 0 = 1 and for any n 1 by

c n = 1 n(n + 1)(n + 2)(n + 3)
.

When |x| < 1, the series S(x) writes as follows

S(x) = 47 36 - 5 12x + 1 6x 2 + (1 -x) 3 log(1 -x) 6x 3 (8)

Remark 3.5

The C A,B cannot be bounded above by some constant that does not depend on A and B, as can be seen hereunder in the proof. Indeed, we show that if a and b are positive integers,

ψ(n, 0 a , 0 b ) ∼ 1 3 S(1) r a r b - 1 r a - 1 r b + 1 S(1) 1 n 3
as n goes to infinity. In particular, ψ(n, 0, 0 n ) tends to the positive constant 13 6 .

Proof of Proposition 3.4.

For any finite words A and B in case i) of Proposition 3.2, one deals with U(x) = ((1x)S(x)) -1 which has 1 as a unique dominant singularity. Indeed, 1 is the unique dominant singularity of S, so that the dominant singularities of U are 1 or zeroes of S contained in the closed unit disc. But S does not vanish on the closed unit disc, because for any z such that |z| 1,

|S(z)| 1 - n 1 1 n(n + 1)(n + 2)(n + 3) = 1 -(S(1) -1) = 17 18 . Since M(x) = S(x) -S(1) (x -1)S(x) = S(1)U(x) - 1 1 -x ,
the unique dominant singularity of M is 1, and when x tends to 1 in the unit disc, (8) leads to

M(x) = A(x) - 1 6S(1) (1 -x) 2 log(1 -x) + O (1 -x) 3 log(1 -x)
where A(x) is a polynomial of degree 2. Using the classical transfer theorem (see Flajolet and Sedgewick [7, section VI]) based on the analysis of the singularities of M, we get

ψ(n -1, w1, 1w ′ ) = [x n ]M(x) = 1 3S(1) 1 n 3 + o 1 n 3 .
The cases ii), iii), iv) and v) of Proposition 3.2 are of the same kind, and we completely deal with case iii).

Case iii): words of the form A = 0 a and B = 0 b , a, b 1. As shown in Proposition 3.2, one has to compute the asymptotics of the n-th coefficient of the Taylor series of the function

M a,b (x) := S(1) 1 r a r b n 1 r a+b+n-1 x n + U(x) S(1)R a (x)R b (x) r a r b x a+b-2 -S(x) . ( 9
)
The contribution of the left-hand term of this sum is directly given by the asymptotics of the remainder

r n = k n c k = 1 3n(n + 1)(n + 2) = 1 3n 3 + O 1 n 4 .
By means of singularity analysis, we deal with the right-hand term

N a,b (x) := U(x) S(1)R a (x)R b (x) r a r b x a+b-2 -S(x) .
Since 1 is the only dominant singularity of S and U and consequently of any R a , it suffices to compute an expansion of N a,b (x) at x = 1. It follows from ( 8) that U, S and R a admit expansions near 1 of the forms

U(x) = 1 S(1)(1 -x) + polynomial + 1 6S(1) 2 (1 -x) 2 log(1 -x) + O(1 -x) 2 , S(x) = polynomial + 1 6 (1 -x) 3 log(1 -x) + O(1 -x) 3 ,
and

R a (x) = polynomial + 1 6 (1 -x) 3 log(1 -x) + O(1 -x) 3 .
Consequently,

N a,b (x) = 1 6 1 r a + 1 r b - 1 S(1) (1 -x) 2 log(1 -x) + O(1 -x) 2
in a neighbourhood of 1 in the unit disc so that, by singularity analysis,

[x n ]N a,b (x) = - 1 3 1 r a + 1 r b - 1 S(1) 1 n 3 + o 1 n 3 . Consequently (9) leads to ψ(n -1, 0 a , 0 b ) = [x n ]M a,b (x) ∼ 1 3 S(1) r a r b - 1 r a - 1 r b + 1 S(1) 1 n 3
as n tends to infinity, showing the mixing inequality and the non uniformity. The remaining cases ii), iv) and v) are of the same flavour. ⊓ ⊔

Mixing of the factorial infinite comb

Consider now the second Example in Section 2, that is the probabilized infinite comb defined by ∀n ∈ N, c n = 1 (n + 1)! .

With previous notations, one gets

S(x) = e x -1 x and U(x) = x (1 -x)(e x -1)
.

Proposition 3.6

The VLMC defined by the factorial infinite comb has a uniform exponential mixing of the following form: there exists a positive constant C such that for any n 1 and for any finite words A and B,

|ψ(n, A, B)| C (2π) n .

Proof.

i) First case of mixing in Proposition 3.2: A = A ′ 1 and B = 1B ′ .

Because of Proposition 3.2, the proof consists in computing the asymptotics of [x n ]M(x). We make use of singularity analysis. The dominant singularities of M(x) = S(x) -S(1) (x -1)S(x) are readily seen to be 2iπ and -2iπ, and

M(x) ∼ 2iπ 1 -e 1 -2iπ • 1 1 -z 2iπ .
The behaviour of M in a neighbourhood of -2iπ is obtained by complex conjugacy. Singularity analysis via transfer theorem provides thus that where P a+b is an entire function. In this last formula, the brackets contain an entire function that vanishes at 1 so that the dominant singularities of M a,b are again those of S -1 , namely ±2iπ. The expansion of M a,b (x) at 2iπ writes thus

[x n ]M(x) ∼ n→+∞ 2(e -1) 1 + 4π 2 1 2π n ǫ n where ǫ n = 1 if n is even 2π if n is odd.
M a,b (x) ∼ 2iπ -S(1)P a (2iπ)P b (2iπ) 1 -2iπ • 1 1 -x 2iπ
which implies, by singularity analysis, that

[x n ]M a,b (x) ∼ n→+∞ 2ℜ 1 -e 1 -2iπ • P a (2iπ)P b (2iπ) (2iπ) n .
Besides, the remainder of the exponential series satisfies n a

x n n! = x a a! 1 + x a + O( 1 a ) ( 10
)
when a tends to infinity. Consequently, by Formula (6), P a (2iπ) tends to 2iπ as a tends to infinity so that one gets a positive constant C 1 that does not depend on a and b such that for any n 1,

|ψ(n, A, B)| C 1 (2π) n .
iii) Third case of mixing: A = 0 a and B = 0 b . This time, one has to compute

[x n ]M a,b (x) with M a,b (x) := S(1) 1 r a r b n 1 r a+b+n-1 x n + U(x) S(1)R a (x)R b (x) r a r b x a+b-2 -S(x)
the first term being an entire function. Here again, the dominant singularities of M a,b are located at ±2iπ and

M a,b (x) ∼ 2iπ -S(1)R a (2iπ)R b (2iπ) (1 -2iπ)r a r b (2iπ) a+b-2 • 1 1 -x 2iπ
which implies, by singularity analysis, that

ψ(n -1, A, B) ∼ n→+∞ 2ℜ 1 -e 1 -2iπ • R a (2iπ)R b (2iπ) r a r b (2iπ) a+b-2 1 (2iπ) n .
Once more, because of [START_REF] Pittel | Asymptotic growth of a class of random trees[END_REF], this implies that there is a positive constant C 2 independent of a and b and such that for any n 1,

|ψ(n, A, B)| C 2 (2π) n .
iv) and v): both remaining cases of mixing that respectively correspond to words of the form A = A ′ 10 a , B = 0 b and A = 0 a , B = 0 b 1B ′ are of the same vein and lead to similar results. ⊓ ⊔

Height and saturation level of suffix tries

In this section, we consider a suffix trie process (T n ) n associated with an infinite random word generated by an infinite comb. A precise definition of tries and suffix tries is given in section 4.1. We are interested in the height and the saturation level of such a suffix trie.

Our method to study these two parameters uses a duality property à la Pittel developed in Section 4.2, together with a careful and explicit calculation of the generating function of the second occurrence of a word (in Section 4.3) which can be achieved for any infinite comb. These calculations are not so intricate because they are strongly related to the mixing coefficient and the mixing properties detailed in Section 3.

More specifically, we look at our two favourite examples, the logarithmic comb and the factorial comb. We prove in Section 4.5 that the height of the first one is not logarithmic but polynomial and in Section 4.6 that the saturation level of the second one is not logarithmic either but negligibly smaller. Remark that despite the very particular form of the comb in the wide family of variable length Markov models, the comb sources provide a spectrum of asymptotic behaviours for the suffix tries.

Suffix tries

Let (Y n ) n≥1 be an increasing sequence of sets. Each set Y n contains exactly n infinite words. A trie process (T n ) n≥1 is a planar tree increasing process associated with (Y n ) n 1 . The trie T n contains the words of Y n in its leaves. It is obtained by a sequential construction, inserting the words of Y n successively. At the beginning, T 1 is the tree containing the root and the leaf 0 . . . (resp. the leaf 1 . . . ) if the word in Y 1 begins with 0 (resp. with 1). For n ≥ 2, knowing the tree T n-1 , the n-th word m is inserted as follows. We go through the tree along the branch whose nodes are encoded by the successive prefixes of m; when the branch ends, if an internal node is reached, then the word is inserted at the free leaf, else we make the branch grow comparing the next letters of both words until they can be inserted in two different leaves. As one can clearly see on Figure 2 a trie is not a complete tree and the insertion of a word can make a branch grow by more than one level. Notice that an internal node exists within the trie if there are at least two words in the set starting by the prefix associated to this node. This indicates why the second occurrence of a word is prominent.

Let m := a 1 a 2 a 3 . . . be an infinite word on A = {0, 1}. The suffix trie T n (with n leaves) associated with m, is the trie built from the set of the n-th first suffixes of m, that is Y n = {m, a 2 a 3 . . . , a 3 a 4 . . . , . . . , a n a n+1 . . .}.

For a given trie T n , we are mainly interested in the height H n which is the maximal depth of an internal node of T n and the saturation level ℓ n which is the maximal depth up to which all the internal nodes are present in T n . Formally, if ∂T n -→ denotes the set of leaves of T n ,

H n = max u∈Tn\∂Tn |u| ℓ n = max j ∈ N| #{u ∈ T n \ ∂T n , |u| = j} = 2 j .
See Figure 3 for an example.

Duality

Let (U n ) n 1 be an infinite random word generated by some infinite comb and (T n ) n 1 be the associated suffix trie process. We denote by R the set of rightinfinite words. Besides, we define hereunder two random variables having a key role in the proof of Theorem 4.8 and Theorem 4.9. This method goes back to Pittel [START_REF] Pittel | Asymptotic growth of a class of random trees[END_REF].

Let s ∈ R be a deterministic infinite sequence and s (k) its prefix of length k. For n ≥ 1, where "s (k) is in T n \ ∂T n " stands for: there exists an internal node v in T n such that s (k) encodes v. For any k 1, T k (s) denotes the number of leaves of the first tree "containing" s (k) . See Figure 4 for an example. Thus, the saturation level ℓ n and the height H n can be described using X n (s):

X n (s) := 0 if s (1) is not in T n max{k 1 | the word s (k) is already in T n \ ∂T n }, T k (s) := min{n 1 | X n (s) = k},
ℓ n = min s∈R X n (s) and H n = max s∈R X n (s). (11)
Moreover, X n (s) and T k (s) are in duality in the following sense: for all positive integers k and n, one has the equality of the events

{X n (s) k} = {T k (s) n}. (12)
The random variable T k (s) (if k 2) also represents the waiting time of the second occurrence of the deterministic word s (k) in the random sequence (U n ) n≥1 , i.e. one has to wait T k (s) for the source to create a prefix containing exactly two occurrences of s (k) . More precisely, for k 2, T k (s) can be rewritten as Notice that T k (s) denotes the beginning of the second occurrence of s (k) whereas in [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF], τ (2) s (k) denotes the end of the second occurrence of s (k) , so that

T k (s) = min n 1 U n U n+1 . . . U n+k-1 = s (k) and ∃!j < n such that U j U j+1 . . . U j+k-1 = s (k) .
τ (2) s (k) = T k (s) + k. (13)
More generally, in [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] , for any r 1, the random return times τ (r) (w) is defined as the end of the r-th occurrence of w in the sequence (U n ) n 1 and the generating function of the τ (r) is calculated. We go over these calculations in the sequel.

Return time generating functions

Proposition 4.7 Let k 1. Let also w = 10 k-1 and τ (2) (w) be the end of the second occurrence of w in a sequence generated by a comb defined by (c n ) n 0 . Let S and U be the ordinary generating functions defined in Section 3.1. The probability generating function of τ (2) (w) is

Φ (2) w (x) = c 2 k-1 x 2k-1 U(x) -1 S(1)(1 -x) 1 + c k-1 x k-1 (U(x) -1) 2 .
Furthermore, as soon as n 1 n 2 c n < ∞, the random variable τ (2) (w) is squareintegrable and

E(τ (2) (w)) = 2S(1) c k-1 + o 1 c k-1 , Var(τ (2) (w)) = 2S(1) 2 c 2 k-1 + o 1 c 2 k-1 . (14)
Proof. For any r 1, let τ (r) (w) denote the end of the r-th occurrence of w in a random sequence generated by a comb and Φ (r) w its probability generating function. The reversed word of c = α 1 . . . α N will be denoted by the overline c := α N . . . α 1 We use a result of [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] that computes these generating functions in terms of stationary probabilities q (n) c . These probabilities measure the occurrence of a finite word after n steps, conditioned to start from the word c. More precisely, for any finite words u and c and for any n 0, let

q (n) c (u) := π U n-|u|+|c|+1 . . . U n+|c| = u U 1 . . . U |c| = c .
It is shown in [START_REF] Cénac | Context trees, variable length Markov chains and dynamical sources[END_REF] that, for |x| < 1, Φ (1) w (x) =

x k π(w) (1x)S w (x) and for r ≥ 1, Φ (r) w (x) = Φ (1) w (x) 1 -

1 S w (x) r-1 where S w (x) := C w (x) + ∞ n=k q (n) ←- pref (w) (w)x n , C w (x) := 1 + k-1 j=1 1 {w j+1 ...w k =w 1 ...w k-j } q (j) ←- pref (w) (w k-j+1 . . . w k ) x j .
In the particular case when w = 10 k-1 , then

←- pref (w) = w = 0 k-1 1 and π(w) = c k-1
S(1) . Moreover, Definition (4) of the mixing coefficient and Proposition 3.2 i) imply successively that

q (n) ←- pref (w) (w) = π U n-k-|w|+1 . . . U n+k = w U 1 . . . U k = ←- pref (w) = π(w) ψ n -k, ←- pref (w), w + 1 = π(w)S(1)u n-k+1 = c k-1 u n-k+1 ,
This relation makes more explicit the link between return times and mixing. This leads to

n k q (n) ←- pref (w) (w)x n = c k-1 x k-1 n 1 u n x n = c k-1 x k-1 U(x) -1 .
Furthermore, there is no auto-correlation structure inside w so that C w (x) = 1 and

S w (x) = 1 + c k-1 x k-1 U(x) -1 .
This entails

Φ (1) w (x) = c k-1 x k S(1)(1 -x) 1 + c k-1 x k-1 U(x) -1
and Φ (2) w (x) = Φ (1) w (x) 1 -

1 S w (x) = c 2 k-1 x 2k-1 U(x) -1 S(1)(1 -x) 1 + c k-1 x k-1 U(x) -1 2
which is the announced result. The assumption n 1 n 2 c n < ∞ makes U twice differentiable and elementary calculations lead to (Φ (1) w ) ′ (1) =

S(1) c k-1 -S(1) + 1 + S ′ (1) S(1) , (Φ (2) w ) ′ (1) = (Φ (1) w ) ′ (1) + S(1) c k-1 , (Φ (1) w ) ′′ (1) = 2S(1) 2 c 2 k-1 + o 1 c 2 k-1
and (Φ (2) w )

′′ (1) = 6S(1) 2 c 2 k-1 + o 1 c 2 k-1
, and finally to (14).

⊓ ⊔

Logarithmic comb and factorial comb

Let h + and h -be the constants in [0, +∞] defined by

h + := lim n→+∞ 1 n max ln 1 π (w) and h -:= lim n→+∞ 1 n min ln 1 π (w) , (15)
where the maximum and the minimum range over the words w of length n with π (w) > 0. In their papers, Pittel [START_REF] Pittel | Asymptotic growth of a class of random trees[END_REF] and Szpankowski [START_REF] Szpankowski | Asymptotic properties of data compression and suffix trees[END_REF] only deal with the cases h + < +∞ and h -> 0, which amounts to saying that the probability of any word is exponentially decreasing with its length. Here, we focus on our two examples for which these assumptions are not fulfilled. More precisely, for the logarithmic infinite comb, (2) implies that π(10 n ) is of order n -4 , so that

h - lim n→+∞ 1 n ln 1 π (10 n-1 ) = 4 lim n→+∞ ln n n = 0.
Besides, for the factorial infinite comb, π(10 n ) is of order 1 (n+1)! so that

h + lim n→+∞ 1 n ln 1 π (10 n-1 ) = lim n→+∞ n! n = +∞.
For these two models, the asymptotic behaviour of the lengths of the branches is not always logarithmic, as can be seen in the two following theorems, shown in Sections 4.5 and 4.6. Theorem 4.9 (Saturation level of the factorial infinite comb) Let T n be the suffix trie built from the n first suffixes of the sequence generated by a factorial infinite comb. Then, the saturation level ℓ n of T n satisfies: for any δ > 1, almost surely, when n tends to infinity,

ℓ n ∈ o log n (log log n) δ .
The dynamic asymptotics of the height and of the saturation level can be visualized on Figure 5. The number n of leaves of the suffix trie is put on the x-axis while heights or saturation levels of tries are put on the y-axis. Plain lines represent a logarithmic comb while long dashed lines are those of a factorial comb (mean values of 25 simulations). Short dashed lines represent a third infinite comb defined by the data c n =

1 3 n-1 k=1 1 3 + 1 (1+k) 2
for n 1. Such a process has a uniform exponential mixing, a finite h + and a positive h -as can be elementarily checked. As a matter of consequence, it satisfies all assumptions of Pittel [START_REF] Pittel | Asymptotic growth of a class of random trees[END_REF] and Szpankowski [START_REF] Szpankowski | Asymptotic properties of data compression and suffix trees[END_REF] implying that the height and the saturation level are both of order log n. Such assumptions will always be fulfilled as soon as the data (c n ) n satisfy lim n c 1/n n < 1; the proof of this result is left to the reader. One can notice the height of the logarithmic comb that grows as a power of n. The saturation level of the factorial comb, negligible with respect to log n is more difficult to highlight because of the very slow growth of logarithms.

These asymptotic behaviours, all coming from the same model, the infinite comb, stress its surprising richness.

Height for the logarithmic comb

In this subsection, we prove Theorem 4.8. Consider the right-infinite sequence s = 10 ∞ . Then, T k (s) is the second occurrence time of w = 10 k-1 . It is a nondecreasing (random) function of k. Moreover, X n (s) is the maximum of all k such that s (k) ∈ T n . It is nondecreasing in n. So, by definition of X n (s) and T k (s), the duality can be written ∀n, ∀ω, ∃k n , k n X n (s) < k n + 1 and T kn (s) n < T kn+1 (s).

(16)

Claim: lim n→+∞ X n (s) = +∞ a.s. (17)
Indeed, if X n (s) were bounded above, by K say, then take w = 10 K and consider T K+1 (s) which is the time of the second occurrence of 10 K . The choice of the c n in the definition of the logarithmic comb implies the convergence of the series n n 2 c n . Thus (14) holds and E[T K+1 (s)] < ∞ so that T K+1 (s) is almost surely finite. This means that for n > T K+1 (s), the word 10 K has been seen twice, leading to X n (s) K + 1 which is a contradiction. We make use of the following lemma that is proven hereunder.

With notations (16), because of (17), the sequence (k n ) tends to infinity, so that (T kn (s)) is a subsequence of (T k (s)). Thus, (18) implies that

For all η > 0, write

T k (s) k 4+η = T k (s) -E(T k (s)) k 4+η + E(T k (s)) k 4+η
. The deterministic part in the second-hand right term goes to 0 with k thanks to (19), so that we focus on the term T k (s) -E(T k (s)) k 4+η

. For any ε > 0, because of Bienaymé-Tchebychev inequality,

P |T k (s) -E[T k (s)]| k 4+η > ε V ar(T k (s)) ε 2 k 8+2η = O 1 k 2η .
This shows the convergence in probability in Lemma 4.10. Moreover, Borel-Cantelli Lemma ensures the almost sure convergence as soon as η > 1 2 . ⊓ ⊔ Remark 4.11 Notice that our proof shows actually that the convergence to +∞ in Theorem 4.8 is valid a.s. (and not only in probability) as soon as δ > 1 36 .

Saturation level for the factorial comb

In this subsection, we prove Theorem 4.9.

Consider the probabilized infinite factorial comb defined in Section 2 by ∀n ∈ N, c n = 1 (n + 1)! .

The proof hereunder shows actually that ℓn log log n log n n is an almost surely bounded sequence, which implies the result. Recall that R denotes the set of all rightinfinite sequences. By characterization of the saturation level as a function of X n (see [START_REF] Rissanen | A universal data compression system[END_REF]), P (ℓ n k) = P (∃s ∈ R, X n (s) k) for all positive integers n, k. Duality formula ( 12) then provides P (ℓ n k) = P (∃s ∈ R, T k (s) n)

P (T k ( s) n)
where s denotes any infinite word having 10 k-1 as a prefix. Markov inequality implies ∀x ∈]0, 1[, P (ℓ n k + 1) P τ (2) (10 k-1 ) < n + k Φ

(2)

10 k-1 (x) x n+k (21)
where Φ

(2) 10 k-1 (x) denotes as above the generating function of the rank of the final letter of the second occurrence of 10 k-1 in the infinite random word (U n ) n≥1 . The simple form of the factorial comb leads to the explicit expression U(x) = x (1-x)(e x -1) and, after computation, Φ

(2)

10 k-1 (x) = e x -1 e -1 • x 2k-1 1 -e x (1 -x) k! (e x -1) (1 -x) + x k-1 1 -e x (1 -x) 2 . (22)
In particular, applying Formula (22) with n = (k -1)! and x = 1 -1 (k-1)! implies that for any k 1,

P(ℓ

  10 a and B = 0 b 1B ′ where A ′ and B ′ are any finite words and a + b 1, then

  10 a and B = 0 b where A ′ is any finite words and a, b 0, then

|w|=n- 1 π

 1 (10 a w0 b 1) = π(10 a+n-1+b 1) + 1≤i≤j≤n-1 |w|=n-1 f (w)=i,l(w)=j π(10 a w0 b 1).

  With Proposition 3.2, the asymptotics of the mixing coefficient comes from singularity analysis of the generating functions M A,B . Proposition 3.4 The VLMC defined by the logarithmic infinite comb has a nonuniform polynomial mixing of the following form: for any finite words A and B, there exists a positive constant C A,B such that for any n 1, |ψ(n, A, B)| C A,B n 3 .

ii)

  Second case of mixing: A = A ′ 10 a and B = 0 b 1B ′ . Because of Proposition 3.2, one has to compute [x n ]M a,b (x) with M a,b (x) := S(1) c a+b c a c b P a+b (xP a (x)P b (x) -S(x) ,

Figure 2 :

 2 Figure 2: Last steps of the construction of a trie built from the set (000 . . . , 10 . . . , 1101 . . . , 001, . . . , 01110 . . . , 1100 . . . , 01111 . . .).

Figure 3 :

 3 Figure 3: Suffix trie T 10 associated with the word 1001011001110 . . .. Here, H 10 = 4 and ℓ 10 = 2.

Figure 4 :

 4 Figure 4: Example of suffix trie with n = 20 words. The saturation level is reached for any sequence having 1000 as prefix (in red); ℓ 20 = X 20 (s) = 3 and thus T 3 (s) 20. The height (related to the maximum of X 20 ) is realized for any sequence of the form 110101 . . . (in blue) and H 20 = 6. Remark that the shortest branch has length 4 whereas the saturation level ℓ n is equal to 3.

Theorem 4 . 8 (n n 1 4

 481 Height of the logarithmic infinite comb) Let T n be the suffix trie built from the n first suffixes of a sequence generated by a logarithmic infinite comb. Then, the height H n of T n satisfies ∀δ > 0, H -δ -→ n→∞ +∞ in probability.

Lemma 4 . 10

 410 For s = 10 ∞ ,

  since the height of the suffix trie is larger than X n (s), Proof of Lemma 4.10. Combining (13) and (14) shows thatE(T k (s)) = E(τ (2) (w))k = 19 9 k 4 + o(k 4 )(19)and Var(T k (s)) = Var(τ(2) (w)) = 361 162 k 8 + o(k 8 ).

1 .x log x-x 1 + O 1 x

 11 1)!+k . Consequently, P ℓ (k-1)! k + 1 = O(k -2 ) is the general term of a convergent series. Thanks to Borel-Cantelli Lemma, one gets almost surely lim n→+∞ ℓ n! n Let Γ -1 denote the inverse of Euler's Gamma function, defined and increasing on the real interval [2, +∞[. If n and k are integers such that (k + 1)! n (k + 2when x goes to infinity, leads to the equivalent Γ -1 (x) ∼ +∞ log x log log x , which implies the result. ⊓ ⊔

Acknowledgements

The authors are very grateful to Eng. Maxence Guesdon for providing simulations with great talent and an infinite patience. They would like to thank also all people managing two very important tools for french mathematicians: first the Institut Henri Poincaré, where a large part of this work was done and second Mathrice which provides a large number of services.


André Galligo

Improved Budan-Fourier Count for Root Finding

Keywords: real univariate polynomial, real root isolation, refinement, Budan-Fourier theorem, Descartes rule, virtual roots, Budan table, Newton process, multiple roots, discretization, separation bound

. The set of virtual roots provide a good real substitute to the set of complex roots; it depends continuously on the coefficients of f . We will describe a root isolation method by a subdivision process based on a generalized Budan-Fourier count, fast evaluation and Newton like approximations. Our algorithm will provide isolating intervals for all augmented virtual roots of f . For a polynomials with integer coefficients of length size τ = Õ(n), its bit cost is in Õ(n 5 ). We rely on a new connexity property of the Budan table of f which collects the signs of the iterated derivatives of f .

Introduction

Real or complex root finding of a univariate real polynomial is one of the most classical problem. It re-appears periodically since the 19th century and there is an extensive bibliography on that subject see [START_REF] Mcnamee | A bibliography on roots of polynomials[END_REF], and [START_REF] Pan | Solving a polynomial equation: some history and recent progress[END_REF]. During the last decade, in relation with the applications in Computer Aided Design, the attention focused on the subdivision methods inspired by Vincent's classical algorithm, see e.g. [START_REF] Mourrain | The Bernstein basis and real root isolation[END_REF], the use of Descartes's rule through homographies (Moebius transforms) and the corresponding representation of real numbers by continued fractions see e.g. [START_REF] Tsigaridas | I: On the complexity of real root isolation using continued fractions[END_REF]. Another novelty is the use of representation of real numbers by dynamic long dyadic approximations called bitstream see [START_REF] Rouillier | Efficient isolation of polynomials real roots[END_REF], and of a secant-like method for accelerating the subdivision process see [START_REF] Abbott | Quadratic interval refinement for real roots[END_REF]. While completing this article, we have seen on arXiv a paper not yet published by Sagraloff [START_REF] Sagraloff | When Newton meets Descartes: A simple and fast algorithm to isolate the real roots of a polynomial[END_REF] which, with these tools, makes important progresses. Sagraloff presents a subdivision algorithms with the same order of bit complexity bounds than the "considered complicated" almost optimal ones proposed by Schoenague and Pan [START_REF] Pan | Solving a polynomial equation: some history and recent progress[END_REF]; although the "complicated" algorithms compute all the complex roots of f .

Our approach is not only conceptually simple and "visual", is also promising and we believe that it has the potential to also meet the quasi optimal complexity. We suggest to avoid a systematic use of Moebius transforms (in the way Descarte's splitting condition is applied), to replace complex roots by "augmented virtual" roots (see below), to apply Newton-Raphson approximation schemes when the derivatives are available, to rely on fast Taylor shifts and evaluations for univariate polynomials. Indeed this last family of algorithms is nowadays well understood and is available as basic commands on long arithmetic packages of some computer algebra systems, see e.g. [START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF].

In the 19th century the Budan-Fourier theorem, which counts signs variations of a sequence and was followed by the invention of Sturm sequences, was considered as a breakthrough. Subdivision methods, which exploits the ordered structure of the real numbers, are widely applied for calculating good approximations of solutions of polynomial equations or intersections of surfaces in many applied sciences. However the analysis of their complexity, hence efficiency, relies on the algebraic nature of the inputs. The geometric dictionary in complex algebraic geometry between invariants readable on equations and features of varieties is ultimately based on the fact that a polynomial of degree n admits n roots. This is not the case for real roots, and makes real algebraic geometry more complicated. A natural strategy for studying properties of real algebraic varieties is to consider simultaneously roots of iterated derivatives of the input. A first conceptual progress was achieved by Gonzales-Vega, Lombardi, Mahé in [START_REF] Gonzales-Vega | Virtual roots of real polynomials[END_REF] when they introduced the concept of the n virtual roots (counted with multiplicities) of a degree n polynomial f , to provide a good "real" substitute to complex roots: The ordered sequence of virtual roots depend continuously on the coefficient of f .

The table containing the signs of all the derivatives of a polynomial f is called, in this paper, its Budan table. It is called after Budan de Boilorant [START_REF] De | Nouvelle méthode pour la résolution des équations numériques d'un degré quelconque[END_REF], who competed with the famous J. Fourier [START_REF] Fourier | Analyse des équations déterminées[END_REF] to provide a proof for the so-called Budan-Fourier theorem. Budan's approach was developed further in a recent work of D. Bembe [START_REF] Bembé | An algebraic certificate for Budan?s theorem[END_REF].

As in [START_REF] Galligo | Roots of the Derivatives of some Random Polynomials[END_REF] (see also [START_REF] Bembé | A: Virtual roots of real polynomials and fractional derivatives[END_REF]) we identify the table with a rectangle formed by positive and negative blocks and consider it as a 2D object. The idea developed in this article is to consider successive approximations of the shape of a Budan table, or of portions of this table, defined by their intersection with grids. This will be done in two steps, the first grid is defined as preprocessing taking into account the expected complexity bound of the whole process, then the grid is adapted following a Newton like procedure. We present simple conditions which, when they are satisfied, define what we call a "valid" discretization. However in our subdivision algorithm, the grid is not statically defined but is the result of a dynamic divide and conquer process, based on an improved Budan-Fourier count.

We adopt a 2D point of view and re-interpret the classical Budan-Fourier bisection method, in such a way that in the "valid" intervals, we get an exact count of real roots (like with Sturm sequences). The presentation of our approach is conceptually simpler when we restrict to the generic case where the roots of all derivatives are two by two distinct, but can form clusters. We are able to analyze and compute isolating intervals for these complicated situations. Then, a slight generalization of our techniques allow to deal with the most general case where either the input polynomial or any of its derivative can have multiple roots; the key tool is a priori given separation bounds.

The paper is organized as follows. Section 2 gives the definitions and the properties of Budan tables and (augmented) virtual roots of a (P) polynomial. and illustrate them with some examples, Section 3 provides conditions to be satisfied for a good discretization of a Budan table (possibly truncated) by a grid; i.e such that an improved Budan-Fourier count gives the exact number of real roots of f in an interval. Section 4 presents the Newton approximation schemes to refine the grid in order to satisfy the previous criteria. Section 5 presents our root finding algorithms of all virtual roots of f and addresses complexity issues. Section 6 reports and comments some experiments.

Notations, separation bounds and arithmetics

R denotes the field of real numbers, R[x] denotes the ring of real univariate polynomial, and f a monic polynomial of degree n. E n denotes the set of all monic polynomial of degree n, which is identified to R n . We will also use the notation f (0) for f .

In most part (but not all part) of this paper, f and all its derivatives are assumed square-free. We need two separation bounds to be able to distinguish, after subdivisions, by evaluations on the border of an interval where f has a single root between three elementary situations: "no root of f ", "a double root" or "two distinct roots". The first separation bound, denoted by s := 2 -N , minors the distance between any two roots of g(x) when g denotes f or any of its derivatives; it allows to certify that an interval with length smaller than s does not contain two distinct roots of such a g. With these notations, it remains to be able to certify that there is no tangent point or equivalently to be able to minor the distance between the graph of g in this interval and the x axis. Assuming that the graph of g (or -g) is convex and taking the intersection of the two tangents to the graph at the border of the interval, the graph is contained in a triangle. We want to certify that the x axis does not touch this triangle. The second separation bound, denoted by t := 2 -N , minors the distance between two points of the interval, such that the triangle touch the x axis, in case when there is no intersection and no tangency with the x axis.

To simplify the analysis and its presentation. in the case of integer arithmetic, we will assume that the length of the coefficients is bounded by τ = Õ(n) and that N = Õ(n 2 ). We will also assume that N = Õ(n 2 ). These assump-tions are pessimistic and will be discussed in the conclusion. A process with a quadratic convergence stopping at the separation bound will have a number of steps bounded by Õ(log(n)), use long integer of length bounded by Õ(n 2 ) for representing the abscissas, and long integer of length bounded by Õ(n 3 ) to represent their evaluation by the polynomials, hence will have a bit cost of Õ(n 4 ). Since we aim to locate by Newton-like processes the n virtual roots of f , the arithmetic worst cost cannot be lower than Õ(n 4 ). This is usually considered as the target bound for the root isolation problem(see previous references). Therefore we will make free use of pre-processing having a lower cost. Of course in an optimized implementation one should be more careful.

Finally let us notice that our approach (by evaluation) is also well adapted for bitstream or approximate computations with big-floats. So one can use successively both approximate and exact computations: approximate computations with big-floats is often very efficient providing an approximate result, but since the set of virtual roots depend continuously on the input coefficients, this can be later refined and certified with exact computations. Therefore this work is also a contribution to SNC (Symbolic Numerical Computation).

Definitions and results

In this section, we first recall classical facts (see e.g. [START_REF] Rahman | Analytic theory of polynomials[END_REF]) or in [START_REF] Gonzales-Vega | Virtual roots of real polynomials[END_REF] and [START_REF] Coste | Generalized Budan-Fourier theorem and virtual roots[END_REF], then in the second subsection we report results from [START_REF] Galligo | Budan Tables of Real Univariate Polynomials[END_REF], while in the third subsection we present a useful ingredient for the developments in the next sections.

Facts on Budan tables and virtual roots

Definition 2.1 Let f be a monic univariate polynomial of degree n. The Budan table of f is the union of n + 1 infinite rectangles of height one L i := R × [i -1/2, i + 1/2[ for i from 0 to n, called rows. For i from 0 to n, each row L i is the union of a set of open rectangles (possibly infinite), separated by vertical segments. We color in black the rectangles corresponding to negative values of the (n -i)-th derivative f (n-i) of f , and we color in gray the rectangles corresponding to positive values of f (n-i) .

Remark 2.2

1. Once we know the coefficients of f , the real roots of all its derivatives are contained in an interval [-2 M , 2 M ] for some integer M . So the table is in fact finite, and when we say ∞ we mean 2 M . 2. Since f is assumed monic, every infinite right rectangle of each row is gray.

3. Since f (n) is a positive constant, the row L 0 is a gray infinite rectangle.

4. The first (infinite) rectangle of each row L i is alternatively gray or black, depending on the parity of i: it is gray if n -i is even.

5. We are interested by the connected components of the union of the closures of the gray rectangles; and respectively for the black rectangles.

It is clear that there is a gray connected component containing the infinite right rectangles of all rows. The other connected components (gray or black) are said bounded on the right.

A classical descriptor attached to a Budan table is the function V f (x) of the real indeterminate x with values in the set of integers N, it counts the number of sign changes in the sequence formed by f and its derivatives evaluated at x. Definition 2.3 For a sequence (a 0 , . . . , a n ) ∈ (R \ {0}) n+1 the number of sign changes V(a 0 , . . . , a n ) is defined inductively in the following way:

V(a 0 ) := 0; V(a 0 , . . . , a i ) := V(a 0 , . . . , a i-1 ) if a i-1 a i > 0, V(a 0 , . . . , a i-1 ) + 1 if a i-1 a i < 0.
To determine the number of sign changes of a sequence (a 0 , . . . , a n ) ∈ R n+1 , delete the zeros in (a 0 , . . . , a n ) and apply the previous rule. (V of the empty sequence equals 0).

We notice that the function V f is computed from the Budan table of f , but two different tables (of two polynomials f and g) may have the same function V f = V g . Therefore the Budan table is a finer invariant than V f attached to the polynomial f . Proposition 2.4 (Budan-Fourier theorem) Let f ∈ R[X] be monic of degree n. Then,

• V f (-∞) = n, V f (∞) = 0.
• Near a real root c of multiplicity k of f , which is not a root of another derivative of f , V f decreases by k when x moves from c -h to c + h, for sufficiently small positive h .

• Near a real root c of multiplicity k of f (m) , which is not a root of another non successive derivative of f , the following happens: If k is even, V f decreases by k.

If k is odd, V f decreases by the even integer k + s 1 s 2 , where s 1 and s 2 are the signs at c of f (m-1) and f (m+k) .

• Near c, a real root of several non successive derivative of f , V f decreases by the sum of the quantities corresponding to each of them.

• Near the other points of R, V f is constant. The function V f is decreasing (but not strictly) on R.

• For a, b ∈ R with a < b, the number of real roots of f in the interval ]a, b] counted with multiplicities is at most V f (b) -V f (a). Moreover the defect is an even integer.

Definition 2.5

Let f be a monic real polynomial of degree n. The x value of the rightmost upper segment of a connected component (either gray or black) of the Budan table of f is called a virtual root of f . Any real root (in the usual sense) of f is a virtual root of f . Any multiple real root (in the usual sense) of any derivative of f is also a virtual root of f . The virtual multiplicities are counted as follows:

• the multiplicities of events appearing along a same x-value are added,

• the multiplicity of a simple root of f counts 1,

• the multiplicity of a simple virtual non real root (i.e. it is not a multiple root of a derivative of f ) counts 2,

• the multiplicity of a multiple root of f of order k counts k,

• the multiplicity of a multiple virtual non real root which is a multiple root of order k of a derivative of f counts k if k is even, and otherwise k + s 1 s 2 where s 1 and s 2 are the signs at c of f (m-1) and f (m+k) .

Budan-Fourier theorem (Proposition 2.4) implies that f admits n virtual roots counted with multiplicities. Moreover the following result holds.

Proposition 2.6 Let f ∈ R[X] be monic of degree n. Let y 1 ≤ • • • ≤ y n the

ordered virtual roots of f , repeated according to their multiplicities, and

y 0 = -∞, y n+1 = ∞. Then we have for 1 ≤ r ≤ n + 1, with y r-1 = y r , x ∈ [y r-1 , y r [⇐⇒ V(f (x), f (x), . . . , f (n) (x)) = n + 1 -k (resp. for r = 1 the interval x ∈] -∞, y 1 [).

Theorem 2.7 ([17], see also [8]) The ordered sequence of virtual roots of a monic polynomial f depend continuously on the coefficients of f .

The virtual roots of f and f satisfy the "classical" interlacing property.

Generic case

In this subsection we assume a condition (P), generically satisfied. This will ease the presentation of our new tools. Then we will consider the general case. We introduce and study several data, attached to f and its Budan table.

Definition 2.8 A polynomial g in R[x] satisfies condition (P) if and only if:

each derivative of g has simple roots, and all these roots are two by two distinct. A monic polynomial satisfying this condition will be called a (P)-polynomial.

Obviously, the parities of the number of real roots and degree of a (P)polynomial are equal. The property (P) is generically satisfied. Moreover, the set of monic polynomials in E n , identified with R n , satisfying (P) form a semialgebraic set of E n .

Generic Budan tables

Now, we determine the features of the Budan table B of a generic monic polynomial f with degree n. Definition 2.9 We say that a table B with (n+1) rows L i formed by rectangles of alternating colors, gray and black, separated by vertical segments, is a GB table of degree n if it satisfies the following properties.

• The row L 0 is a gray infinite rectangle. The infinite rightmost rectangle of each row is gray. The first (infinite) rectangle of each row L i is alternatively gray or black, depending on the parity of i: it is gray if n -i is even.

• If i is even (resp. odd) the number of rectangles on the row L i is even (resp. odd).

• Let (l + 1) be the number of rectangles of the top row L n , then l ≤ n and n -l is an even number 2p. There are l + p + 1 same-color-connected components of B. Each non first rectangle of L i , i > 0 is connected on the left to a rectangle of the same color of the row L i-1 .

• The previous item is true, replacing n by any m, 0 < m < n, and B by the table formed by the lower m + 1 rows.

Theorem 2.10 Let f be a (P)-polynomial of degree n, and let m ≤ n be the number or real (simple) roots of f . Then m and n have the same parity, n = m + 2p and the Budan table of f is a GB table of degree n.

Example 1

The following Figure 1 shows a GB table of degree 10 with m = 4, p = 3.

Lemma 2.11

Let f be a (P)-polynomial of degree n. The x value of the rightmost upper segment of a connected component (either gray or black) of the Budan table of f is a virtual root of f . Any real root (in the usual sense) of f is a virtual root. Let m ≤ n be the number or real (simple) roots of f , and let n -m = 2p. There are p virtual non real roots of f , they are of multiplicity two; each of them is a root of some derivative of f of positive order.

We recover that f admits n virtual roots, counted with multiplicities.

Definition 2. [START_REF] Fourier | Analyse des équations déterminées[END_REF] We call augmented virtual root of f the pair (y, k) formed by a virtual root of f and the order of the derivative of f which vanishes at y, i.e.

f (k) (y) = 0.
The augmented virtual roots of f only depend on the Budan table B of f . Virtual roots of a GB table are well defined. Proposition 2.13 Let f be a (P)-polynomial of degree n. By Rolle's theorem between two successive roots a < b of some derivative f (m) with 0 ≤ m ≤ n -2, (or in R if f (m) has no root), there is an odd number 2r + 1 of roots (X 1 < ... < X 2r+1 ) of the next derivative f (m+1) . Then the r roots with an even index (X 2 , ...X 2r ) are virtual non real roots of f . Similarly if a is the smallest (resp. the largest) root of f (m) , in the infinite interval ] -∞, a[ (resp. ]a, ∞[) there is an even number of roots 2r of roots (X 1 < ... < X 2r ) of the next derivative f (m+1) . Then the r roots with an odd index (X 1 , ...X 2r-1 ), (resp. the k roots with an even index (X 2 , ...X 2r )) are virtual non real roots of f . For each augmented virtual non real root (y, k) of f , we have

f (k-1) (y)f (k+1) (y) > 0.

General case

In the general case where the condition (P) is not necessarily satisfied, the definitions of (augmented) virtual roots can be generalized and moreover a continuity result holds. Definition 2. [START_REF] Galligo | Roots of the Derivatives of some Random Polynomials[END_REF] We call augmented virtual root of f a triple (y, k, r) formed by a virtual root of f, an order of derivation k 0 and a multiplicity r 1 such that we have:f (k) (y) = 0, ..., f (k+r-1) (y) = 0 but f (k+r) (y) = 0; and if k > 0,f (k-1) (y) = 0. Different augmented virtual root of f can correspond to the same virtual root of f, moreover they only depend on the Budan table of f . Note that if f is a (P) polynomial then the last coordinates of each triple is 1, so it can be identified with a pair, as we did in the previous subsection.

Each augmented virtual root (y i , k i , r i ) may have a multiplicity ≥ 3. See Figure 2 the Budan table of a degree 7 polynomial f with a simple real root, a Figure 3: a TWO-truncated Table real root of multiplicity 4, and a virtual non real root of multiplicity 2 according to the rule of Proposition 2.4 (but is a root of multiplicity 3 of f ).

Theorem 2.10 which was proved when condition (P) holds, generalizes to the case when this conditions is violated. The only difference is that the end points of some same-color-connected components have the same first coordinate.

Notice that Figure 2 resembles Figure 1, except that the end points of some same-color-connected components have the same first coordinate.

Truncated Budan table

Let f be a monic polynomial of degree n, we analyze the properties of a sub table P := P (f, a, b, u, v) of its Budan table BT . P is delimited on the x axis by two real numbers a and b which are not root of a derivative of f , a < b, and on the second coordinate by two integers u and v such that 0 ≤ u < v ≤ n.

Let us denote by W (x) := W (f, u, v)(x) the function giving the number of sign changes in the sequence formed by the derivatives

f (n-k) , with u ≤ k ≤ v evaluated at x. Let m 2 := W (b) -W (a)
Among the p 1 + l 1 real roots of f (n-u) between a and b, p 1 are virtual roots of f and l 1 are not.

Theorem 2.15

With the previous notations, let m := l 1 + m 2 . Then the sub table P has l + p same-color-connected components, with m = l + 2p; the top row of P has l + 1 rectangles (their l right segments indicate the l roots of f (n-v) between a and b); and the p other ends of the same-color-connected components indicate the virtual non real roots of f (n-v) in P (hence virtual non real roots of f ).

The proof of this theorem uses the same reasonning and is a very similar to the proof of Theorem 2.10. Note that if a rectangle corner of the lower row of P (which corresponds to a root of f (n-v) between a and b) is surrounded above and on the right by 2 rectangles of opposite color, then in the Budan table of f it is also surrounded below by a rectangle of opposite color; hence it corresponds to a virtual non real root of f .

For an illustration, consider any portion of Figures 1 and2, or Figure 3 which represents a truncated Budan table analyzed in subsection 5.1.

Discretized Budan table

In this section we consider a (P) polynomial f and approximations of its truncated Budan table via intersections with grids. Then sketch on an example a strategy of computation of augmented virtual roots, hence of real roots of f .

We will consider two kinds of grids, the first one is either formed by an arithmetic progression on the x-axis times an interval of integer (degrees of derivatives of f ) or by a simple bisection method; while the second one is obtained via a Newton like refining process.

In the illustrations, to be clearer, we present only small portions of the grid; we replace the sign + by a grey solid box and the sign -by a black empty box.

A simple example

We consider a polynomial of degree n = 256 given by f = n i=0 rand() n i x i ; rand() is a function which returns a signed random integer with a uniform distribution between -N and N , where N is a large integer, such a polynomial is often called a random SO(2) polynomial. It is almost surely a generic polynomial which satisfies condition (P). We choose 50 points forming an arithmetic progression between -0.1 and 0.1 and consider f and its first 15 derivatives. So we get a discretized 2D picture with 800 pixels shown in Figure 4. Then we proceed similarly replacing -0.1 and 0.1 by -1 and 1 to get Figure 5.

Let us notice that the Budan-Fourier count for these two intervals gives:

V f (0.1) -V f (-0.1) = 220 , V f (1) -V f (-1) = 246.
Therefore for this example a great part of the virtual roots are between -0.1 and 0.1.

It turns out that the first discretization gives a faithful representation of the truncated Budan table (see below the definition of a "valid" discretization), in contrast with the coarser second one. We assume that we have already certified (by performing inductively on the lower degrees the same construction) that f (15) has only 1 root between -0.1 and 0.1 with the location shown on Figure 4,i With the help of some technical tool (Newton approximations until the separation bound or via interval arithmetic), we certify the signs above and below, hence that they are indeed virtual roots of f . Finally Theorem 2.15 asserts that there are at most 16-2*7=2 real roots of f between -0.1 and 0.1. We see on Figure 4 that there are exactly 2 such roots located in two intervals: ] -0.1 + 5 * 0.2/50, -0.1 + 6 * 0.2/50[ and ] -0.1 + 25 * 0.2/50, -0.1 + 26 * 0.2/50[. Now, let us comment the information given by the grid on Figure 5. It is too coarse and should be refined to become a faithful representation of the corresponding Budan table. Indeed we expect that all derivatives of f are square-free, but this property is not well represented on Figure 5 since in some places the same small interval serves for delimiting the root of some f (k) and of its derivative. E.g. the interval starting at the integer grid coordinates: [38, 240], [38,241]. A refinement or a partition will solve this point, but partition are less expansive from a computational point of view. In this case, first assume (by induction) that the last row is correctly discretized. A partial Budan-Fourier count, the W (f, u, v) function of the previous section, shows that there are no non real virtual roots between x 27 and x 50 , similarly between x 0 and x 23 . Then there are 7 candidates virtual roots in [x 23 , x 27 ]. This interval is included in [-0.1, 0.1] and it now requires the refinement performed in Figure 4 to separate the augmented virtual roots.

A Mignotte polynomial

Now let us describes what happens with the very different case of a Mignotte polynomial: n = 100, f := x n + 2(5x -1). It is well known that this polynomial has only 2 very near real roots. It has also clearly a virtual non real root of multiplicity n -2 at x = 0, since its second derivative is n(n -1)x n-2 . An arithmetic grid with 300 elements between -0.1 and 0.3 will guess the virtual non real root of multiplicity n -2 but will not distinguish between two near by real roots and virtual root near 0.2 as shown on Figure 6.

Valid discretization 4.1 Conditions

In order to express conditions for a good discretization of the Budan table of a P polynomial f , we consider truncated tables of height 1 (resp. 2), i.e. P (f, a, b, n, n -1) (resp. P (f, a, b, n, n -2)) and their 2 points discretizations, they give a 2 by 2 grid (resp. 2 by 3). We first assume that we know by induction that if the signs f (a) and f (b) are equal (resp. different) there are no root (resp. 1 root) between a and b. Then we classify the configurations of In the upper row in Figure 7, the 4 cases correspond to valid configurations of signs. In the lower row, the first case that we call K 1 may correspond to a virtual root (possibly of higher multiplicity) or two roots, the second and third cases that we call K 2 should be refined since they locate a root of a polynomial and of its derivative in a same interval, while the last case is impossible.

[i, k], B[i + 1, k], B[i, k + 1], B[i + 1, k + 1]
], the last 3 cases in the lower row of Figure 7 (or the cases obtained by inverting the signs + and -), never happen.

In this section we consider a first computational strategy which amounts to check inductively by increasing degrees when either of the cases K 1 or K 2 appear, and refine the grid near these points until they are replaced by a case of the upper row of Figure 7 until it is certified that the considered case K 1 corresponds to a virtual root. We are aware that this strategy performs unnecessary computations, but it clearly describe our approach; in the next section we will present a fast algorithm relying on exclusion tests.

The logical condition for case K 2 is

B[i, k] * B[i + 1, k] = -1 and B[i, k + 1] * B[i + 1, k + 1] = -1.
The logical condition for case K 1 is for e = 1 or e = -1,

B[i, k] = -e, B[i + 1, k] = e and B[i, k + 1] = B[i + 1, k + 1] = e.
Moreover, we have the following simple but useful result.

Lemma 4.2

For each of the 3 cases of type K 1 or K 2 , only one configuration of a 2 by 3 rectangular grid is allowed, if we assume (by induction) that the lower part of the discretization is valid. These configurations have equal signs on the lower third row, they are shown on Figure 8.

Approximation theorem

Let f be a monic (P)-polynomial of degree n, consider a sub table P := P (f, a, b, u, v) of its Budan table BT . P is delimited on the x axis by two real numbers a and b which are not root of a derivative of f , a < b, and on the second coordinate by two integers u and v such that 0 ≤ u < v ≤ n. Denote by W (x) := W (f, u, v)(x) the function giving the number of sign changes in the sequence formed by the derivatives

f (n-k) , with u ≤ k ≤ v evaluated at x. Let m 2 := W (a) -W (b).
Assume that f (n-u) has r 1 real roots between a and b.

Theorem 4.3 With the previous notations and definitions, consider the discretization of the table

P := P (f, a, b, u, v
) by a grid of points. Assume that the cases K 2 never happen and that all the squares of P verifying K 1 are certified to correspond to virtual roots. Assume that the discretization shows r 1 changes of signs on the lower row (which corresponds to the real roots of f (n-u) ) and p 1 of them verify

K 1 . Let m 1 := r 1 -p 1 and m := m 1 + m 2 .
Then the sub table P has l+p same-color-connected components, with m = l+2p; the top row of P has l change of signs (they locate the l roots of f (n-v) between a and b); and there are p same-color-connected components ends indicting the virtual non real roots of f (n-v) in P (hence virtual non real roots of f ).

Proof: It is a direct consequence of the definitions and of the previous results.

Refinements

Here, we present for each of the two cases K 1 and K 2 described in the previous section, a refinement scheme relying on Newton-Raphson algorithm.

For K 2

The two sub cases corresponding to K 2 have symmetric shapes, so it suffices to consider the one with positive convexity. Take a derivative g of f corresponding to case K 2 on an interval The hypothesis K 2 says: g(a) > 0, g (a) < 0, g(b) < 0, g (b) > 0, and g [a,b] > 0. We look for γ such that a < γ < b and g(γ) ≤ 0, g (γ) < 0.

We apply the classical Newton-Raphson algorithm to g(x) -g(b) starting from a, and stop when the root of g is isolated. More precisely:

c := a -g(a)-g(b) g (a)
if g(c) ≤ 0, then RETURN(γ := c), else update a := c end if. Note that by positive convexity we have always g (c) > 0. Since the algorithm has quadratic convergence and the distance between a root of g and a root of g is assumed greater than the separation bound s. We want to determine if the graph of g cuts the x axis between a and b. This will not be the case if the second coordinate of the intersection point C is positive. While if the graph cuts the x axis, we detect it via a γ such that a < γ < b and g(γ) ≤ 0, and more likely when the decision is tough g (γ) will be "near" 0.

For

We will test the two possibilities with an iterative algorithm which alternates two computations: a Newton step for g starting at a (or b) which updates a (or b) and a step computing the coordinates [c, L] of the intersection C of the two tangents and proceeds as follows.

If g(c) < 0 then we return that the interval corresponds to two distinct roots. If L > 0 then we return that the interval corresponds to a virtual root, moreover without tangency; else if g (c) ≤ 0 then update a := c, else update b := c end if end if.

If the size of the interval becomes lower than the two separation bounds then there is a virtual root with tangency.

An alternative test that g remains positive on a small interval [a, b] when we do not know that g > 0 but when b -a < 2 -l is much smaller than 1 is to test if

g(a) -(b -a)g (a) > k=n k=2 g (k) (a) k! 2 -lk .
The last possibility to consider is when we already had a virtual root which is a multiple root of f or of a derivative of f . Then we rely on the separation bounds.

Algorithms for a (P) polynomial

In this section, we describe a fast algorithm for virtual roots isolation for a (P) polynomial, it already contains the main difficulties since we also consider the possibility of clusters. The general case addressed in the next section will follow. We present a (pessimistic) estimation of the worst case bit complexity of our algorithm for a polynomial f with integer coefficients of length τ = Õ(n), then it is known that the roots of f and all its derivatives are in [2 -M , 2 M ] with τ = Õ(n). Moreover, we assume given a separation bound between each pair of roots of each derivative of f , equal to s = 2 -N with N = Õ(n 2 ). (That is one of the point that we find pessimistic and would like to improve in a future work).

Subroutine TWO

We first consider the basic case where we have an interval I = [a, b], such that the Budan-Fourier count BF (f, I), (the difference between the two signs variations) is equal to 2. We want to determine if there are 2 distinct real roots of f , else compute the augmented virtual root i.e determine a degree k > 0 such that for an interval I = [a , b ] included in I, f (n-k+1) keeps a constant sign on I and f (n-k) has one simple root in I with BF (f (n-k) , I) = 1. Assume that neither a nor b is a root of a derivative of f . Let k 1 be the greatest integer such that BF (f (n-k) , I) < 2. Necessarily k ≥ k 1 . Then for any k > k 1 , f (n-k ) has either 0 or 2 roots in I, hence f (n-k ) (a) and f (n-k ) (b) have the same sign. Moreover by monotony, f (n-k -1) (a) < f (n-k -1) (b). These remarks allow to restrict the interval ]k 1 , k 2 ] where k should be searched. See Figure 3. We proceed by a binary search so in at most log(k 2 -k 1 ) steps. At each step, we test with an integer k , let g := f (n-k ) (x), and apply to g Newton steps starting from b to search a value c < b such that g(c) < 0. If c is found in ]a, b[ then we update k 1 : k and a := c. Else we update k 2 : k . Notice that g decreases from b to c, (or to a if there is no c) and that we have g > 0 on this interval [c, b], except for k = k + 1, detected at the end of the loop and where we proceed as in subsection 4.4.2.

Preprocessing

We perform 3M = Õ(n) steps of the following bisection, construct a binary tree of segments I = At the end of the preprocessing, the set A contains u intervals

I i = [a i , b i ] such that b i -a i < 2 -2M and 1≤i≤u BF (I i ) ≤ n.
Moreover following [START_REF] Pan | Solving a polynomial equation: some history and recent progress[END_REF], we assume that these intervals are separated by Õ(n) times their size and will call them "clusters" of virtual roots.

The intervals in BB will correspond to interval isolating a root of f . We will later apply Subroutine TWO to the intervals in B.

The total bit cost the preprocessing is of order Õ(n 4 ).

Processing

Now instead of performing bisection of I ∈ A by the middle m(I), we perform the 2 following steps which aims bounding the clusters of augmented virtual roots in rectangles of type

I × [k 1 , k 2 ]
. Now in A and B we collect not only the intervals I but indeed the products

I × [k 1 , k 2 ]
. Just after the preprocessing we start with k 1 = 0 and k 2 = n. 1) Cutting the bottom and refining: For a chosen

I × [k 1 , k 2 ]
in A, we first compute the lower degree k + 1 such that the Budan Fourier count BF (f (n-k-1) (I) becomes greater or equal to 2. Therefore f (n-k) admits a simple root on I and all its derivatives have one or zero (simple) root on I. We propose to perform two Newton steps (or as usual in numerical recipes, mix it with some bisection to avoid to encounter a cycle, since we cannot certify convexity, but still achieve quadratic convergence) from a and b to compute a and b with a ≤ a ≤ b ≤ b. Then update the sets A and B as explained above.

2) Cutting the top: If for some I × [k 1 , k 2 ] in A the total multiplicities of the cluster of virtual roots in I, detected by the changes in the signs variations, is greater than k 2 -k 1 , it means that in the cluster should be divided at least in two parts. Starting from the top, a probable "weakest link" is the the value k 3 where the partial difference of signs variations W (f, u, v) on I (see section 3) pauses. So we perform one of the same sign test K 1 on I. If it succeeds, we delete

I × [k 1 , k 2 ] from A, then add I × [k 1 , k 3 ] and I × [k 3 + 1, k 2 ] in A. Note that f (n-k3-1
) admits a simple root on I.

We stop either when A is empty or if the sizes of all remaining intervals I are smaller than the separation bound.

Illustration with Figure 4

We consider the example shown in Figure 4, assume that at the end of the preprocessing I = [-0.1, 0.1] and I × [240, 256] ∈ A; and I × [0, 240] has already been processed, in particular the real root shown in the last row is simple. We start with a cluster of S(I, 240, 256) = 16 virtual roots counted with multiplicities, and cut it into smaller clusters. In this illustration the values obtained by Newton steps are denoted using the letter N .

Step 1: we perform two Newton steps, N (f (15) , x 0 ) → x 20 , N (f (15) , x 50 ) → x 40 . so we easily decompose I × [240, 256] in 3 parts, the left one and the right one are added to B so let us concentrate on the new cluster. We update A adding I 1 := [x 20 , x 40 ] and S([x 20 , x 40 ], 240, 256) = 13.

Step 2: we perform two Newton steps, N (f (15) , x 20 ) → x 26 , N (f (15) Finally check that all the augmented virtual roots of f has been well separated and certified.

Output the list of the augmented virtual roots of f . Output the list of the number of real roots (and their multiplicities) of each derivative of f .

Worst case complexity, bit costs

The depth of a Newton iteration tree is O(log(n)), i.e. Õ [START_REF] Abbott | Quadratic interval refinement for real roots[END_REF]. Hence the size of the total subdivision tree is in Õ(n). However the last evaluations involve smaller intervals so are more costly, assuming pessimistically that all process go till the separation bound and that each separation bound is in 2 -N with N = Õ(n 2 ); we arrive at a total bit cost of Õ(n 5 ) bits. This lags by a factor n behind the fastest (complex) root finding algorithms of Schoenague and Pan [START_REF] Pan | Solving a polynomial equation: some history and recent progress[END_REF].

Remark 5.1 Notice that for the same order of computational bit cost, i.e Õ(n 5 ), one can get a set of isolating intervals of the roots for all derivative of f , and the shape of the Budan table of f . Our bit costs also lags by a factor n behind the new subdivision algorithm for square-free polynomials presented by Sagraloff in his very recent arXiv preprint [START_REF] Sagraloff | When Newton meets Descartes: A simple and fast algorithm to isolate the real roots of a polynomial[END_REF]. However this paper relies on a more subtle separation bound which is the product of all the separation bounds of all the complex roots of f ; he uses it for his fine interpretation of Obrechkoff theorems.

We do not have yet in our setting any estimation of a similar smart separation bound for the augmented virtual roots of f . However, if we restrict to the problem of separating simple real roots of a square-free polynomial, we can implement an early detection subroutine explained below, and discard most branches of the subdivision tree.

In that case we believe (but we did not prove yet) that our approach could also take advantage of a smart separation bound, drop the extra factor n and also meets quasi optimal complexity bounds.

Early detection

Assuming that f is square-free, after the pre-processing we can concentrate on finding only the real roots of f . We proceed as in the previous subsection but we discard from the set A all rectangles

[a, b] × [k 1 , k 2 ] such that k 2 < n.
With the assumptions and b -a < 2 (-L) with L = Õ(n), in "many" cases the same-sign test will work and certify in the early stage of the subdivision process that k 2 < n. Remark 5.2 If we suspect that k < n an alternate procedure could be to apply a Moebius transform (i.e. a translation composed with an inversion). A generic inversion will not lower the multiplicity of a multiple root (or of a compact cluster) of f , however it will spread out and separate the multiple root (or of a compact cluster) of a derivative of f which are not root of f . One Moebius transform by cluster is enough. Notice that other differentiable bijections, e.g polynomial transform, will do the same effect but they will increase the degree of f .

General case

We follow essentially the same algorithm and analysis as in the previous subsection.

The only difference is that the multiplicity of a virtual root can be greater than 2. In the Budan table of f this means as in Figure 3, that some ends of same-sign-connected components instead of being surrounded by points with the opposite sign can have just above them a zero. In other words the corresponding graph admits a tangency.

This can give rise to singular Newton like approximations, but in the previous section we already took it into account with the clusters, and proposed a safe bottom-up process.

Therefore the only difference will appear at the very last steps. at this point we rely on the second separation bound (we assumed a lower bound 2 -N with N = Õ(n 2 )) which allow to distinguish the tangency hence the equality of the x values of the augmented virtual roots. At the same order bit complexity cost: Õ(n 5 ).

The algorithm similarly output all virtual roots of f but also their multiplicities.

Experiments

We have implemented a prototype of our algorithm which needs to be tuned and optimized, nevertheless it produces interesting informations.

The usual benchmarks where f is either a classic random polynomial or Laguerre or Wilkinson or Mignotte polynomials are somehow rough since either the separation bounds are not small or one of their first derivative is a polynomial with many well separated real roots, or if they admit clusters there are only one or two of them. Therefore the efficiency of our approach reduces to the efficiency of the used Taylor shifts which computes the f (x + a) in the preprocessing. We are not aware of other benchmarks, but it would be interesting to develop a great variety of them.

A composed example

Here, to present and comment our algorithm on a complete example, we consider the following polynomial (composed with the previous ones with the following notations: f := W ilkison(n) := 0≤i≤n (x -i) Kac( 96)) := 0≤i≤n a i x i where a i are random real numbers following a standard normal distribution.

f := W ilkison(32) * M ignotte(64) + round(10 10 * Kac(96)).

It has degree n = 96 with integer coefficients of of about 30 digits, hence τ about n. We chose a rather small degree to ease the description. We denote by V f (a) the Budan Fourier count at the value a (i.e. the number of sign changes in the sequence of derivatives evaluated at a).

After few (less than log(n)) checks. We see that V f (-1/2) = 96, V f (7/2) = 1, and V f (4) = 0. So f admits a simple real root between 3.5 and 4 and all the subdivisions will happen in [-1/2, 7/2] an interval of size 4. In order to illustrate the potential of our approach we perform a "short" preprocessing.

Preprocessing

We construct a bisection tree of depth 7 i.e. log(128), we collect the intervals of size about 1/32 in 3 sets: BB contains the intervals [a, b] such that V (a) -V (b) = 1 isolating a simple root of f , B contains the intervals [a, b] such that V (a) -V (b) = 2 isolating a virtual non real root of f , or two simple roots of f , (they will be solved later in a post processing), A contains the intervals [a, b] such that V (a) -V (b) ≥ 3, that we call "clusters". So it remains 96-30-4 = 62 virtual roots counted with multiplicities. They are included in the 7 following intervals of A, for each of them we indicate the

Early detection

If we are only interested by computing the real roots, there is a test to get rid of a cluster of virtual roots which is not a multiple root of f . In this example consider the "compact" cluster ([-1/32, 0], 36).

We have f (0) 5.10 35 , but then the signs of the last 30 coefficients of f alternate, therefore on [-1/32, 0] their contribution can only increase positively the approximation of f (x) given by the Taylor expansion. To estimate a lower bound of f (x) we only need to find an upper bound of the "remainder". Since |x| ≤ 2 -5 a bound is i=n i=31 |coeff(f, x, i)|2 -5i , but this sum is expected to be very small, it is indeed about 10 -35 . Therefore after the preprocessing we can forget this cluster and subtract 36 from the cumulated multiplicities to be controlled. More generally, we can take advantage of our knowledge of the signs to bound the remainder in a Taylor expansion.

Conclusion

Although there have been many scientific works (and implementations) on real root finding algorithm via subdivision methods, the Budan-Fourier count which historically initiated the subject was not considered as an efficient tool. The reason is that, in contrast with Sturm count or Descartes rule of signs associated to Moebius transforms, it did not provide a termination criterium. In this article we presented a new subdivision based on an improved Budan Fourier count, extended to the derivatives of the input polynomial f , which now provides such a termination criterium.

Instead of just finding the real root of f , our approach allows to also find all the "augmented virtual roots" of f , and eventually the Budan tree and the Budan table of f . Those are new concepts that we introduced in a previous work and that we consider important abstract data associated to f . We view the proposed algorithm as successive approximations of these data. More precisely we have considered a quad-tree like approximations of truncated rectangles of the Budan table of f , in the spirit of the classical and improved Weyl algorithms as explained in [START_REF] Pan | Solving a polynomial equation: some history and recent progress[END_REF] for complex root finding, and a bottom-up (with increasing degrees) processing.

The semantic of our approach is geometric and quite classical in singular Newton processes: the multiplicities or compact clusters of roots of g are controlled by a derivative of g having simple well separated roots. The depth and size of the subdivision tree are controlled by separations bounds.

The subject of virtual roots is still new and we do not have yet the fine estimates separation bounds known for the complex roots of a polynomial. We obtained a satisfactory bit complexity cost of Õ(n 5 ) for our all procedure.

But if we compare it with the best real root finding algorithms, it lags behind by a factor n, due to the lack of a good global estimate of all the separation bounds.

This will be a subject for future researches. Another direction of research is to study finely the effect of bijective transforms such as Moebius or Graeffe transforms on clusters of augmented virtual roots.

From another point of view, following the philosophy of [START_REF] Labhalla | Espace métrique rationnellement présentés et complexité[END_REF] if the coeficients of the input polynomial f are approximate real numbers known with some precision (or given by some oracle), one can only expect to compute within some precision the virtual roots (or clusters of virtual roots). Our approach allows to achieve this goal.

In a joint paper in preparation, with Mariemi Alonso Garcia, we are applying the approach presented in this article to the important case of fewnomials.

We hope that our approach will be adopted and further developed, even in higher dimensions, by other researchers.

Figure 1 :Figure 2 :

 12 Figure 1: A Budan table of degree 10

Figure 4 :Figure 5 :

 45 Figure 4: A valid discretized table

Figure 6 :

 6 Figure 6: a Mignotte polynomial

Figure 7 :

 7 Figure 7: Distributions of signs

Figure 8 :

 8 Figure 8: Signs for K1 or K2

  [a, b], and assume (by induction) that the lower grid is valid. Consider the points A and B of coordinates A := [a, g(a)], B := [b, g(b)] and the intersection point C of the two tangent lines to the graph of g at A and B. By positive convexity the graph of g is contained in the triangle ACB.

K 1

 1 Take a derivative g of f corresponding to case K 2 on an interval [a, b], and assume (by induction) that the lower grid is valid. Consider the points A and B of coordinates A := [a, g(a)], B := [b, g(b)] and the intersection point C of the two tangent lines to the graph of g at A and B. By positive convexity the graph of g is contained in the triangle ACB. The hypothesis K 1 says: g(a) > 0, g (a) < 0, g(b) > 0, g (b) > 0, and g [a,b] > 0.

  [a, b], and update three sets A, B and BB. When we start Band BB are empty and A contains [2 -M , 2 M ]. Pick I = [a, b] from A, and delete it from A, by fast Taylor shift expand f (x + a) and f (x + b) then compute the Budan-Fourier count S := BF (I), i.e the difference between the two signs variations. If S = 0 then discard I. If S = 1 then put I in BB. If S = 2 then put I in B. Else divide I by its middle m(I) then add [a, m(I)] and [m(I), b] to A.

  .e with the integer grid coordinates [34, 240]. Then we apply Theorem 2.15 to this truncated Budan table.

	Now we see on the Figure 4 that we have only 7 candidate virtual non real
	roots, namely with the integer grid coordinates:
	[24, 245], [24, 249], [27, 254], [30, 247], [32, 251], [44, 244], [30, 241].

We count the signs changes and get W (f, 240, 256)(-0.1) = 15 , W (f, 240, 256)(0.1) = 0 hence with the notations of Theorem 2.15, m:=15+1=16.

  Definition 4.1 A discretization of a truncated Budan table P by a grid [x 1 , ..., x N ] × [u, ..., v] is called valid, if and only if for each 2 by 2 square formed by the signs [B

  , x 40 ) → x 35 . Then the rectangle [x 20 , x 40 ] × [240, 256] is discarded; we set I 2 := [x 20 , x 26 ], I 3 := [x 26 , x 35 ], add I 2 × [245, 256] and I 3 × [240, 256] to A, with S([x 26 , x 35 ], 240, 256) = 8 and S([x 20 , x 26 ], 245, 256) = 5. Now since 256 -245 = 11 > 5 we try to cut I 2 × [245, 256] from the top: since the jumps of W are 5, 4, 4, 4, 4, 4, 3, 2, 2, 1, we test if f keeps the same sign (here positive) on I 2 . Since it is so, we delete I 2 × [245, 256] from A, then add I 2 × [256, 256] to B and I 2 × [245, 255] to A. etc ...

	Post processing:
	We consider all elements I × [k 1 , k 2 ] in B, which have multiplicity 2, and
	apply Subroutine TWO.

  For this example, BB contains 4 simple roots: [[1, 33/32], [63/32, 2], [3, 97/32], [3.5, 4]]. At this stage we cannot certify that they are the only ones. B contains 15 intervals: B := [[3/32, 1/8], [1/8, 5/32], [1/4, 9/32], [3/8, 13/32], [9/16, 19/32], [5/8, 21/32], [21/32, 11/16], [11/16, 23/32], [3/4, 25/32], [25/32, 13/16], [29/32, 15/16], [15/16, 31/32], [31/32, 1], [39/32, 5/4], [29/16, 59/32]].

Acknowledgments

The author thanks Henri Lombardi and Mariemi Alonso Garcia for valuable discussions. The paper was written the author visited the department of Algebra of the University Complutense in Madrid (Spain). This work was partially supported by the contract MathAmSud (11MATH-04-Complexity-Deterministic and probabilistic complexity of algorithms for solving equations) and by the European Marie Curie network SAGA.

value V (a) -V (b).

A := [([-1/16, -1/32], 4), ([-1/32, 0], 36), ([1/16, 3/32], 4), ([3/16, 7/32], 6), ([5/16, 11/32], 4), ([7/16, 15/32], 4), ([17/32, 9/16], 4)].

Processing

For this example, the small clusters can be solved and the new intervals put in the set B. So let us concentrate our description on the more compact cluster ([-1/32, 0], 36).

We compare the two lists of signs at a = -1/32 and b = 0, at their ends we read: [..., +, -, +] and [..., +, +, +] the last derivatives. Hence the degree 1 polynomial f (95) vanishes on that interval, we compute a rough decimal approximation of its root,-0.0098. Then evaluate V (-0.01) = 80, V (-0.009) = 78, and recall that V (-1/32) = 92, V (0) = 56 were computed previously. Therefore after this Newton step, the first compact cluster is replaced in A by 2 new clusters ([-1/32, -0.01], 12) and ([-0.009, 0], 22); [-0.01, -0.009] is put in B.

Again let us e.g. concentrate our description on ([-0.009, 0], 22).

We compare the two lists of signs at a = -0.009 and b = 0. They are rather similar and differ only between the degrees 32 and 53: at b = 0 there are only 21 signs + while at a = -0.009 the signs -and + alternate (as described in section 2 for a multiple virtual root of multiplicity 22). We apply two Newton steps to f (64) which is of degree 32, and compute two approximations, of the simple root of that polynomial in the interval, with a doubled precision of 10 -4 , we found a = -10 -4 and b = 0. Moreover V (-10 -4 ) = 64. Hence we get two new smaller clusters: ([-0.009, -10 -4 ], 14) and ([-10 -4 , 0], 8).

Again we perform Newton steps, the first one with a precision of 10 -4 and the second one with a doubled precision of 10 -8 . etc... When the interval is small enough, here 10 -8 , we check that we found the correct augmented virtual root. It is not a real root of f . And so on, until A is empty and B contains 46 elements. For this cluster we computed the Newton steps with a precision up to 2 -50 , i.e. Õ(n).

Post processing

Here we consider each element

the greatest integer such that f (n-k) has only one root in [a, b]. We assume wlog that there is a safe value k such that W (f, k , k 1 , a) = W (f, k , k 1 , b) + 1 and W (f, k , k 1 +1, a) = W (f, k , k 1 +1, b)+2. Then we apply our subroutine TWO. Its cost depends on the separation distance between the roots of the derivatives.