Victor Y Pan
email: victor.pan@lehman.cuny.edu

Elias P Tsigaridas
email: elias.tsigaridas@inria.fr

Nearly Optimal Computations with Structured Matrices

We estimate the Boolean complexity of multiplication of structured matrices by a vector and the solution of nonsingular linear systems of equations with these matrices. We study four basic most popular classes, that is, Toeplitz, Hankel, Cauchy and Vandermonde matrices, for which the cited computational problems are equivalent to the task of polynomial multiplication and division and polynomial and rational multipoint evaluation and interpolation. The Boolean cost estimates for the latter problems have been obtained by Kirrinnis in [11], except for rational interpolation, which we supply now. All known Boolean cost estimates for these problems rely on using Kronecker product. This implies the d-fold precision increase for the d-th degree output, but we avoid such an increase by relying on distinct techniques based on employing FFT. Furthermore we simplify the analysis and make it more transparent by combining the representation of our tasks and algorithms in terms of both structured matrices and polynomials and rational functions. This also enables further extensions of our estimates to cover Trummer's important problem and computations with the popular classes of structured matrices that generalize the four cited basic matrix classes.

INTRODUCTION

Table 1 displays four classes of most popular structured matrices, which are omnipresent in modern computations for Sciences, Engineering, and Signal and Image Processing. These basic classes have been naturally extended to the four larger classes of matrices, T , H, V, and C, that have structures of Toeplitz, Hankel, Vandermonde and Cauchy types, respectively. They include many other important classes of structured matrices such as the products and inverses of the matrices of these four basic classes, as well as the companion, Sylvester, subresultant, Loewner, and Pick matrices. All these matrices can be readily expressed via their displacements of small ranks [START_REF] Pan | Structured Matrices and Polynomials: Unified Superfast Algorithms[END_REF]Chapter 4], which implies their • Solution of nonsingular linear systems of equations with these matrices in quadratic or nearly linear arithmetic time. These properties enable efficient computations, closely linked and frequently equivalent to fundamental computations with polynomials and rational polynomial functions, in particular to the multiplication, division, multipoint evaluation and interpolation [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF]. Low arithmetic cost is surely attractive, but substantial growth of the computational precision quite frequently affects the known algorithms having low arithmetic cost (see, e.g., [START_REF] Bunch | Stability of methods for solving toeplitz systems of equations[END_REF]). So the estimation of the complexity under the Boolean model is more informative, although technically more demanding.

To the best of our knowledge, the first Boolean complexity bounds for multipoint evaluation are due to Ritzmann [START_REF] Ritzmann | A fast numerical algorithm for the composition of power series with complex coefficients[END_REF]. We also wish to cite the papers [START_REF] Van Der Hoeven | Fast composition of numeric power series[END_REF] and [START_REF] Kobel | Fast approximate polynomial multipoint evaluation and applications[END_REF], although their results have been superceded in the advanced work of 1998 by Kirrinnis, [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF], apparently still not sufficiently well known. Namely in the process of studying approximate partial fraction decomposition he has estimated the Boolean complexity of the multipoint evaluation, interpolation, and the summation of rational functions. He required the input polynomials to be normalized, but actually this was not restrictive at all. We generalize his estimates. For simplicity we assume the evaluation at the points of small magnitude, but our estimates can be rather easily extended to the case of general input. Kirrinnis' study as well as all previous estimates of the Boolean complexity of these computational problems rely on multiplying polynomials as integers, by using Kronecker's product, aka binary segmentation, as proposed in [START_REF] Fischer | String-matching and other products[END_REF]. This implies the d-fold increase of the computational precision for the d-th degree output. The results that we present rely on FFT algorithms for multiplying univariate polynomials and avoid this precision growth.

We represent our FFT-based estimates and algorithms in terms of operations with both structured matrices and polynomial and rational functions. In both representations the computational tasks and the solution algorithms are equivalent, and so the results of [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF] for partial fraction decomposition can be extended to most although not all of these tasks. By using both representations, however, we make our analysis more transparent. Furthermore in Section 7 we extend Table 1. Four classes of structured matrices .

Toeplitz matrices T = (ti-j ) n-1 i,j=0

Hankel matrices H = (hi+j ) n-1 i,j=0

       t0 t-1 • • • t1-n t1 t0 . . . . . . . . . . . . . . . t-1 tn-1 • • • t1 t0              h0 h1 • • • hn-1 h1 h2 . . . hn . . . . . . . . . . . . hn-1 hn • • • h2n-2       Vandermonde matrices V = Vs = s j i n-1 i,j=0
Cauchy matrices C = Cs,t =

1 s i -t j n-1 i,j=0       1 s1 • • • s n-1 1 1 s2 • • • s n-1 2 . . . . . . . . . 1 sn • • • s n-1 n             1 s 1 -t 1 • • • 1 s 1 -tn 1 s 2 -t 1 • • • 1 s 2 -tn . . . . . . 1 sn-t 1 • • • 1 sn-tn      
Kirrinnis' results to the solution of a Cauchy linear system of equations (which unlike [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF] covers rational interpolation) and in Section 7.2 to the solution of Trummer's celebrated problem [START_REF] Gerasoulis | A fast algorithm for trummer's problem[END_REF], [START_REF] Greengard | A fast algorithm for particle simulations[END_REF], [START_REF] Carrier | A fast adaptive multipole algorithm for particle simulations[END_REF], having important applications to mechanics (e.g., to particle simulation) and representing the secular equation, which is the basis for the MPSolve, the most efficient package of subroutines for polynomial rootfinding [START_REF] Bini | Solving secular and polynomial equations: A multiprecision algorithm[END_REF].

Our estimates cover multiplication of the matrices of the four basic classes of Table 1 by a vector and solving Vandermonde and Cauchy linear systems of equations. (These tasks are equivalent to the listed tasks of the multiplication, division, multipoint evaluation and interpolation of polynomials and rational functions.) Expressing the solution of these problems in terms of matrices has a major advantage: it can be extended to matrices from the four larger matrix classes T , H, V, and C. Actually the algorithms for multiplication by vector can be extended quite readily, as we explain in Section 7. There we also briefly discuss the solution of linear systems of equations with the matrices of the cited classes, which can be a natural subject of our further study.

Notation.

In what follows OB, resp. O, means bit, resp. arithmetic, complexity and OB, resp. O, means that we are ignoring logarithmic factors. "Ops" stands for "arithmetic operations". For a polynomial A = d i=0 ai x i ∈ Z Z[x], deg(A) = d denotes its degree and L (A) = τ the maximum bitsize of its coefficients, including a bit for the sign. For a ∈ ( Q, L (a) ≥ 1 is the maximum bitsize of the numerator and the denominator. µ(λ) denotes the bit complexity of multiplying two integers of size λ; we have µ(λ) = OB(λ). 2 Γ is an upper bound on the magnitude of the roots of A. We write ∆α(A) or just ∆α to denote the minimum distance between a root α of a polynomial A and any other root. We call this quantity local separation bound. We also write ∆i instead of ∆α i . ∆(A) = minα ∆α(A) or just ∆ denotes the separation bound, that is the minimum distance between all the roots of A. The Mahler bound (or measure) of A is M (A) = a d |α|≥1 |α|, where α runs through the complex roots of A, e.g. [START_REF] Mignotte | Mathematics for Computer Algebra[END_REF][START_REF] Yap | Fundamental Problems of Algorithmic Algebra[END_REF].

If A ∈ Z Z[x] and L (A) = τ , then M (A) ≤ A 2 ≤ √ d + 1 A ∞ = 2 τ √ d + 1.
If we evaluate a function F (e.g. F = A) at a number c using interval arith-metic, then we denote the resulting interval by [F (c)], provided that we fix the evaluation algorithm and the precision of computing. We write

D(c, r) = {x : |x-c| ≤ r}. f ∈ ( C[x] denotes a λ-approximation to a polynomial f ∈ ( C[x], such that f -f ∞ ≤ 2 -λ . In particular a ∈ (
C denotes a λapproximation to a constant a ∈ ( C such that |a -a| ≤ 2 -λ . lg stands for log.

PRELIMINARIES

Univariate Separation Bounds

The following proposition provides upper and aggregate bounds for the roots of a univariate polynomial. There are various version of these bounds. We use the one presented in [START_REF] Strzeboński | Univariate real root isolation in an extension field[END_REF], to which we also refer the reader for further details and a discussion of the literature. For multivariate separation bounds we refer the reader to [START_REF] Emiris | The DMM bound: Multivariate (aggregate) separation bounds[END_REF].

Proposition 1. Let f = d i=0 aix i ∈ ( C[
x] be a square-free univariate polynomial of a degree d such that a d a0 = 0. Let Ω be any set of k pairs of indices (i, j) such that 1 ≤ i < j ≤ d, let the complex roots of A be 0

< |γ1| ≤ |γ2| ≤ • • • ≤ |γ d |, and let disc(f ) be the discriminant of f . Then |a0| f 2 ≤ |γi| ≤ f 2 |a d | , (1)
(i,j)∈Ω |γi -γj| ≥ 2 k-d-d(d-1) 2 |a0| k M(f ) 1-d-k |disc(f )| ≥ 2 k-d-d(d-1) 2 |a0| k f 1-d-k 2 |disc(f )|. (2) If f ∈ Z Z[x] and the maximum coefficient bitsize is τ then 2 -τ -1 ≤ |γ| ≤ 2 τ +1 , (3)
-lg

(i,j)∈Ω |γi -γj| ≤ 3d 2 + 3dτ + 4d lg d. (4)
The following lemma from [START_REF] Strzeboński | Univariate real root isolation in presence of logarithms[END_REF] provides a lower bound on the evaluation of a polynomial that depends on the closest root and on aggregate separation bounds. Lemma 2. Suppose L ∈ ( C, f is a square-free polynomial, and its root γ1 is closest to L. Then

|f (L)| ≥ |a d | 7 |L -γ1| 6 M (f ) -6 2 lg i ∆ i -6 .

Complex Interval arithmetic

We also need the following bounds for the width of complex intervals when we perform computations with interval arithmetic. We will use them to bound the error when we perform basic computation with complex (floating point) numbers. We refer the reader to [START_REF] Rokne | Complex interval arithmetic[END_REF] for further details.

Proposition 3 (Complex intervals). Given complex intervals I and J, where |I|, resp. |J|, denotes the modulus of any complex number in the complex interval I, resp. J. If 2 -ν ≤ |I| ≤ 2 τ and |J| ≤ 2 σ , then wid(I + J) ≤ 2 wid(I) + 2 wid(J), wid(I J) ≤ 2 τ +1 wid(J) + 2 σ+1 wid(I), and wid(1/I) ≤ 2 4ν+2τ +3 wid(I).

Approximate multiplication of two polynomials

We need the following two lemmas from [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF] on the evaluation of a polynomial at the powers of a root of unity and on polynomial multiplication. A result similar to the first lemma appeared in [21, Section 3] where Bluestein's technique from [START_REF] Bluestein | A linear filtering approach to the computation of discrete fourier transform. Audio and Electroacoustics[END_REF] is applied (see also [START_REF] Knuth | The art of computer programming[END_REF]Chapter 4.3.3,Exercise 16]). We use that lemma to provide a bound on the Boolean complexity of multiplying two univariate polynomials when their coefficients are known up to a fixed precision. An algorithm for this problem appeared in [START_REF] Schönhage | Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coeficients[END_REF]Theorem 2.2] based on employing Kronecker's product, but instead we rely on FFT and the estimates of Corollary 4.1 from [2, Chapter 3].

Lemma 4. Suppose A ∈ ( C[x] of a degree at most d such that A ∞ ≤ 2 τ . Let K = 2 k ≥ d for a positive integer k.
Assume that we know the coefficients of A up to the precision -ℓ -τ -lg K -3; that is the input is assumed to be a polynomial A such that A -A ∞ ≤ 2 -ℓ-τ -lg K-3 ≥ 10. Let ω = exp( 2π K √ -1) denote a K-th root of unity. Then we can evaluate the polynomial A at 1, ω, . . . , ω K-1 in OB(K lg K µ(ℓ + τ + lg K)) such that max 0≤i≤K-1 |A(ω i ) -

A(ω i )| ≤ 2 -ℓ . Moreover, |A(ω i )| ≤ K A ∞ ≤ 2 τ +lg K , for all 0 ≤ i ≤ K -1. Lemma 5. Let A, B ∈ ( C[x] of degree at most d, such that A ∞ ≤ 2 τ 1 and B ∞ ≤ 2 τ 2 .
Let C denote the product AB and let K = 2 k ≥ 2d + 1 for a positive integer k. Write λ = ℓ + 2τ1 + 2τ2 + 5.1 lg K + 4. Assume that we know the coefficients of A and B up to the precision λ, that is that the input includes two polynomials A and B such that

A -A ∞ ≤ 2 -λ and B -B ∞ ≤ 2 -λ . Then we can compute in OB(d lg d µ(ℓ + τ1 + τ2 + lg d)) a polynomial C such that C -C ∞ ≤ 2 -ℓ . Moreover, C ∞ ≤ 2 τ 1 +τ 2 +2 lg K for all i.
Remark 6. In the sequel, for simplicity we occasionally replace the value λ = ℓ + 2τ1 + 2τ2 + 5.1 lg(2d + 1) + 4 by its simple upper bound ℓ + 2τ1 + 2τ2 + 6 lg d + 15 .

APPROXIMATE FFT-BASED POLYNOMIAL DIVISION

In this section we present an efficient algorithm and its complexity analysis for dividing univariate polynomials approximately. This result is the main ingredient of the fast algorithms for multipoint evaluation and interpolation. The evaluation is involved into our record fast real root-refinement, but all these results are also interesting on their own right because, unlike the previous papers such as [START_REF] Schönhage | Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coeficients[END_REF], [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF] and [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF], we keep the Boolean cost bounds of these computations at the record level by employing FFT rather than the Kronecker product and thus decreasing the precision of computing dramatically.

Assume two polynomials s(x) = m i=0 six i and t(x) = n i=0 tix i such that smtn = 0, m ≥ n, and seek the quotient q(x) = m-n i=0 qix i and the remainder r(x) = n-1 i=0 rix i of their division such that s(x) = t(x) q(x) + r(x) and deg(r) < deg(t). Further assume that tn = 1. This is no loss of generality because we can divide the polynomial t by its nonzero leading coefficient. We narrow our task to computing the quotient q(x) because as soon as the quotient is available, we can compute the remainder r(x) = s(x) -t(x) q(x) at the dominated cost by multiplying t(x) by q(x) and subtracting the result from s(x).

The complexity analysis that we present relies on root bounds of t(x), contrary to [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF] where it relies on bounds on the infinity norm of t(x). To keep the presentation selfcontained we copy from [START_REF] Pan | Nearly Optimal Refinement of Real Roots of a Univariate Polynomial[END_REF] the matrix representation of the algorithm, which occupies the next two pages, up to to Lemma 9.

We begin with an algorithm for the exact evaluation of the quotient. Represent division with a remainder by the vector equation

               1 tn-1 1 . . . . . . t1 t0 t1 • • • 1 t0 t1 . . . t1 t0                       qm-n qm-n-1 . . . q1 q0        +              rn-1 rn-2 . . . r0              =                sm sm-1 . . . sn sn-1 . . . s0                .
The first m -n + 1 equations form the following vector equation,

       1 tn-1 1 . . . . . . t1 t0 t1 • • • 1               qm-n qm-n-1 . . . q1 q0        =              sm sm-1 . . . sn+1 sn . . . sm-n-1              ⇔ T q = s, (5)
where q = (qi) m-n i=0 , s = (si) m i=m-n+1 , and T is the nonsingular lower triangular Toeplitz matrix, defined by its first column vector t = (ti) n i=0 ), tn = 1. Write T = Z(t) and Z = Z(e2) where e2 = (0, 1, 0 . . . , 0) T is the second coordinate vector, and express the matrix T as a polynomial in a generator matrix Z = Zn+1 of size (n + 1) × (n + 1) as follows,

Z =            0 . . . 0 1 . . . . . . . . . . . . . . . . . . 0 0 . . . 1 0            , T = Z(t) = t(Z) = n i=0 tiZ i , Z n+1 = O.
The matrix T is nonsingular because tn = 0, and the latter equations imply that the inverse matrix T -1 = t(Z) -1 mod Z n+1 is again a polynomial in Z, that is again a lower triangular Toeplitz matrix defined by its first column. We compute this column by applying a divide and conquer algorithm. Assume that n + 1 = γ = 2 k is a power of two, for a positive integer k. If this is not the case, embed the matrix T into a lower triangular Toeplitz γ × γ matrix t(Zγ ) for γ = 2 k and k = ⌈lg(n + 1)⌉ with the leading (that is northwestern) block T = t(Zγ ), such that t(Zγ) = t(Zγ ) mod Z n+1 γ , compute the inverse matrix and output its leading (n + 1)

× (n + 1) block T -1 . Now represent T as the 2×2 block matrix, T =   T0 O T1 T0  
where T0 and T1 are γ 2 × γ 2 Toeplitz submatrices of the Toeplitz matrix T , T0 is invertible, and observe that

T -1 =   T0 O T1 T0   -1 =   T -1 0 O -T -1 0 T1 T -1 0 T -1 0   . ( 6
)
We only seek the first column of the matrix T -1 . Its computation amounts to solving the same problem for the half-size triangular Toeplitz matrix T0 and to multiplication of each of the γ 2 × γ 2 Toeplitz matrices T1 and T -1 0 by a vector. Let T T I(s) and T M (s) denote the arithmetic cost of s × s triangular Toeplitz matrix inversion and multiplying an s × s Toeplitz matrix by a vector, respectively. Then the above analysis implies that T T I(γ) ≤ T T I(γ/2) + 2T M (γ/2). Recursively apply this bound to T T I(γ/2 g ) for g = 1, 2, . . . , and deduce that T T I(γ) ≤ h g=1 T M (γ/2 g ). The following simple lemma (cf. [16, equations (2.4.3) and (2.4.4)]) reduce Toeplitz-by-vector multiplication to polynomial multiplication and the extraction of a subvector of the coefficient vector of the product, thus implying that T M (s) ≤ cs lg s for a constant c and consequently T T I(γ) < 2cγ lg γ.

Lemma 7. The vector equation

               u0 O . . . . . . . . . . . . u0 um . . . . . . . . . . . . O um                    v0 . . . vn     =              p0 . . . . . . pm . . . pm+n              (7) is equivalent to the polynomial equation m i=0 uix i n i=0 vix i = m+n i=0 pix i . ( 8
)
We wish to estimate the Boolean (rather than arithmetic) cost of inverting a triangular Toeplitz matrix T and then extend this to the Boolean cost bound of computing the vector T -1 s and of polynomial division. So next we assume that the input polynomials are known up to some precision 2 -λ and employ the above reduction of the problem to recursive (approximate) polynomial multiplications.

To study the Boolean complexity of this procedure, we need the following corollary, which is a direct consequence of Lemma 5 and the inequality lg(2d + 1) ≤ 2 + lg d.

Corollary 8 (Bounds for the product P 2 0 P1). Let a polynomial P0 ∈ ( C[x] have a degree d, let its coefficients be known up to a precision 2 -ν , and let P0 ∞ ≤ 2 τ 0 . Similarly, let P1 ∈ ( C[x] have the degree 2d, let its coefficients be known up to a precision 2 -ν , and let P1 ∞ ≤ 2 τ 1 .

Then the polynomial P = P 2 0 P1 has degree 4d, its coefficients are known up to the precision 2 -ν+8τ 0 +2τ 1 +15 lg d+40 , and P 2 0 P1 ∞ ≤ 2 2τ 0 +τ 1 +6 lg d+8 . The following lemma is a normalized version of Lemma 4.4 in [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF].

Lemma 9. Let F, G ∈ ( C[x] such that deg(F ) = m ≥ n = deg(G) ≥ 1, let 2 ρ
be an upper bound on the magnitude of roots of G, and let

F = GQ + R with deg(Q) = m -n and deg(R) = n -1. Then Q ∞ ≤ 2 m+lg m+mρ F ∞ and R ∞ ≤ 2 m+n+lg m+mρ F ∞ .
Proof: To bring the roots inside the unit circle, transform the polynomials by scaling the variable x as follows,

f (x) = F (x 2 ρ ), g(x) = G(x 2 ρ ), q(x) = Q(x 2 ρ ), and r(x) = R(x 2 ρ ). Now apply [11, Lemma 4.4] to the equation f = gq + r to obtain q ∞ ≤ q 1 ≤ 2 m-1 f 1 ≤ 2 m+lg m f ∞ and r ∞ ≤ r 1 ≤ 3 4 2 m+n f 1 ≤ 2 m+n+lg m f ∞ . Combine these inequalities with the equation f ∞ = 2 mρ F ∞ and the inequalities Q ∞ ≤ q ∞ and R ∞ ≤
r ∞ to deduce the claimed bounds.

We will estimate by induction the cost of inverting the matrix T , by using Eq. ( 6) recursively. The proof of the following lemma could be found in the Appendix.

Lemma 10. Let n + 1 = 2 k for a positive integer k and let T be a lower triangular Toeplitz (n + 1) × (n + 1) matrix of. Eq. ( 5), having ones on the diagonal. Let its subdiagonal entries be complex numbers of magnitude at most 2 τ known up to a precision 2 -λ . Let 2 ρ be an upper bounds on the magnitude of the roots of the univariate polynomial t(x)

associated with T . Write T -1 = (T -1 i,j ) n i,j=0 . Then max i,j |T -1 i,j | ≤ 2 (ρ+1)n+lg(n)+1 .
Furthermore, to compute the entries of T -1 up to the precision of ℓ bits, that is to compute a matrix

T -1 = ( T -1 i,j ) n i,j=0 such that maxi,j |T -1 i,j -T -1 i,j | ≤ 2 -ℓ , it is sufficient to know the entries of T up to the precision of ℓ + 10τ lg n + 70 lg 2 n + 8(ρ + 1)n lg n or O(ℓ + (τ + lg n + nρ) lg n) = O(ℓ + τ + nρ) bits. The computation of T -1 costs OB(n lg 2 (n) µ(ℓ+(τ +lg n+ nρ) lg n)) or OB(nℓ + nτ + n 2 ρ).
As usual in estimating the complexity of approximate division we assume that m = 2n to simplify our presentation. Recall that s(x) = t(x) q(x) + r(x).

Theorem 11. Assume s, t ∈ ( C[x] of degree at most 2n and n, such that s ∞ ≤ 2 τ 1 , t ∞ ≤ 2 τ 2 ,
and 2 ρ is an upper bound on the magnitude of the coefficients of t(x). Assume that we know the coefficients of s and t up to a precision λ, that is that the input includes two polynomials s and t such that s -s ∞ ≤ 2 -λ and t -t ∞ ≤ 2 -λ , where λ = ℓ + τ1 + 12τ2 lg n + 80 lg 2 n + 10(ρ + 1)n lg n + 30 or λ = O(ℓ + τ1 + τ2 lg n + n ρ lg n). Let q denote the quotient and let r denote remainder of the division of the polynomials s by t, that is s = t • q + r where deg r < deg t.

Then we can compute in OB(n lg 2 (n) µ(ℓ

+ τ1 + (τ2 + nρ) lg n)) or OB(nℓ + nτ1 + nτ2 + n 2 ρ) two polynomials q and r such that q -q ∞ ≤ 2 -ℓ and r-r ∞ ≤ 2 -ℓ , q ∞ ≤ 2 n+lg n+1+nρ+τ 1 and r ∞ ≤ 2 3n+lg n+1+nρ+τ 1 .
Proof: We compute the coefficients of q(x) using Eq. ( 5), ie q = T -1 s. Each coefficient of the polynomial q comes as the inner product of two vectors, ie qi = n j=0 T -1 i,j sj. From Lemma 10 we know that lg|T

-1 i,j | ≤ n(ρ + 1) + lg n + 1 = N and lg |T -1 i,j -T -1 i,j | ≤ -λ + l2 for l2 = 10τ2 lg n + 70 lg 2 n + 8(ρ + 1)n lg n.
For the coefficients of the polynomials s = 2n j=0 sjx j and t = n j=0 tj x j , we have assumed the following bounds, lg|sj| ≤ τ1, lg |sj -sj | ≤ -λ, lg|tj | ≤ τ2, lg |tj -tj| ≤ -λ, lg|T -1 i,j sj| ≤ τ1 + N , and lg|T -1 i,j sj -T -1 i,j sj | ≤ -λ + ℓ2 + τ1 for all i and j. Therefore

lg q -q ∞ ≤ lg| j T -1 i,j sj - j T -1 i,j sj | ≤ -λ + ℓ2 + τ1 + lg n ≤ -λ + 10τ2 lg n + 70 lg 2 n + 8(ρ + 1)n lg n + τ1 + lg n .
To compute the remainder we apply the formula r(x) = s(x)-t(x)q(x). It involves an approximate polynomial multiplication and a subtraction. For the former we use Lemma 5 and obtain the inequality lg t q -t q ∞ ≤ -λ + 2τ2 + 6 lg n + 26 + ℓ2 + 2N .

Let us also cover the impact of the subtraction. After some calculations and simplifications that make the bounds less scary (albeit less accurate wrt the constant involved), we obtain

lg r -r ∞ ≤ -λ + τ1 + 2τ2 + 6 lg n + 26 + ℓ2 + 2N ≤ -λ + τ1 + 2τ2 + 6 lg n + 26 + 10τ2 lg n + 70 lg 2 n + 8(ρ + 1)n lg n + 2(n(ρ + 1) + lg n + 1)
≤ -λ + τ1 + 12τ2 lg n + 80 lg 2 n + 10(ρ + 1)n lg n + 30 .

By using Lemma 9 we bound the norms of the quotient and the remainder as follows: lg r ∞ ≤ 3n + lg n + 1 + nρ + τ1 and lg q ∞ ≤ n + lg n + 1 + nρ + τ1 . The maximum number of bits that we need to compute with is ℓ + τ1 + 12τ2 lg n + 80 lg 2 n + 10(ρ + 1)n lg n + 30 or

O(ℓ + τ1 + τ2 lg n + lg 2 n + n ρ lg n).
The complexity of computing

T -1 i,j is OB(n lg 2 (n) µ(ℓ + τ1 + τ2 lg n + lg 2 n + n ρ lg n)) or OB(nℓ + nτ1 + nτ2 + n 2 ρ).
According to Lemma 5 the complexity of computing the product t q is OB(n lg(n) µ(ℓ

+ τ1 + τ2 lg n + lg 2 n + n ρ lg n)) or OB(nℓ + nτ1 + nτ2 + n 2 ρ).
Remark 12. We can eliminate the dependence of the bounds of Theorem 11 on τ2 by applying Vieta's formulae and the following inequality,

|t k | ≤ n k (2 ρ ) k ≤ 2 2n+nρ
, where t k is the k-th coefficient of t(x). In this way, after some further simplifications, the required precision is ℓ + τ1 + 150(ρ + 1)n lg n and the complexity bound becomes OB(n lg 2 (n) µ(ℓ+ τ1 + ρ lg n)) or OB(nℓ + nτ1 + n 2 ρ).

MULTIPOINT POLYNOMIAL EVALUA-TION

Problem 1. Multipoint polynomial evaluation. Given the coefficients of a polynomial p(x) = n-1 i=0 pix i and a set of knots t0, . . . , tn-1, compute the values r0 = p(t0), . . . , rn-1 = p(tn-1) or equivalently compute the vector r = V p where r = (ri) n-1 i=0 , p = (pi) n-1 i=0 , and V = (x j i ) n-1 i,j=0 .

In the case where the knots ti = ω i are the n-th roots of 1 for all i, ω = exp(2π √ -1)/n, and V = Ω = (ω ij ) n-1 i,j=0 , Problem 1 turns into the problem of the DFT(v) computation.

Solution:. The Moenck-Borodin algorithm of [START_REF] Moenck | Fast modular transforms via division[END_REF] solves Problem 1 in O(M (n) log n) ops for M (n) in (2.4.1), (2.4.2) based on the two following simple observations. Fact 13. p(a) = p(x) mod (x -a) for any polynomial p(x) and any scalar a. Fact 14. w(x) mod p(x) = (w(x) mod (u(x)p(x))) mod p(x) for any triple of polynomials u(x), p(x), and w(x).

Algorithm 4: the Moenck-Borodin algorithm for multipoint polynomial evaluation.

initialization:. Write k = ⌈log 2 n⌉, m (0) j = x -xj, j = 0, 1, . . . , n -1; m (0) j = 1 for j = n, . . . , 2 k -1 (that is, pad the set of the moduli m (0) j = x -xj with ones, to make up a total of 2 k moduli). Write r (k) 0 = p(x).
Computation:

1. Fan-in process (see Figure 1). Compute recursively the "supermoduli" m

(h+1) j = m (h) 2j m (h) 2j+1 , j = 0, 1, . . . , 2 k-h - 1; h = 0, 1, . . . , k -2.
2. Fan-out process (see Figure 2). Compute recursively the remainders r

(h) j = r (h+1) ⌊j/2⌋ mod m (h) j , j = 0, 1, . . . , min{n, ⌈n/2 h ⌉-1}; h = k -1, k -2, . . . , 0. Output:. p(xi) = r (0) i , i = 0, 1, . . . , n -1.
Let us include a brief outline of the analysis of the algorithm (cf. [START_REF] Moenck | Fast modular transforms via division[END_REF]). To prove its correctness, first apply Fact 14 recursively to obtain that r To estimate the computational cost of the algorithm, represent its two stages by the same binary tree (see Figures 3.1 At each level h of the tree, the algorithm computes 2 k-h products of pairs of polynomials of degree 2 h at stage 1 and 2 k-h remainders of the division of polynomials of degree of at most 2 h+1 by "supermoduli" of degree 2 h . Each time multiplication/division uses

O(M (2 h )) ops for M (n) in (2.4.1), (2.4.2). So we use O(2 k-h M (2 h )) ops at the h-th level and O( k-1 h=0 2 k-h M (2 h )) = O(M (2 k )k) ops at all levels. Re- call that n ≤ 2 k < 2n and obtain the claimed bound of O(M (n) log n) ops.
Remark 15. The fan-in computation at stage 1 depends only on the set {t0, . . . , tn-1} and can be viewed as (costfree) preprocessing if the knot set is fixed and only the polynomial p(x) varies. Similar observations hold for the solution of many other problems in this chapter.

Remark 16. Problem 1 and its solution algorithms are immediately extended to the case where we have m points t0, . . . , tm-1 for m > n or m < n. The solution requires O(E(l)r/l) ops provided l = min{m, n}, r = max{m, n}, and E(l) ops are sufficient for the solution where n = l. E(l) = O(M (l) log l) for a general set {ti} but decreases to O(M (l)), where ti = at 2i + bt i + c for fixed scalars a, b, c, and t and for all i. This also leads to a similar improvement of the estimates for the Boolean complexity [START_REF] Aho | Evaluating polynomials at fixed sets of points[END_REF].

Boolean complexity estimates

In the following two lemmata we present the bit complexity of the fan-in and the fan-out process. These results are of independent interest. We do not estimate the accuracy needed and the bit complexity bound of the algorithm for multipoint evaluation because in Lemma 21 we cover a more general algorithm. Multipoint evaluation is its special case.

Lemma 17 (Complexity of Fan-in process). Assume that we are given n complex numbers xi known up to a pre- Now assume that the claimed bounds hold for k-1, that is deg(m

cision λ = ℓ + (4n -4)τ + 32n -(lg n + 5) 2 -7, that is |xi -xi| ≤ 2 -λ ,
(k-1) i ) = 2 k-1 , lg m (k-1) i ∞ ≤ 2 k-1 τ + 2 k+2 -2k -6, and lg m (k-1) i -m (k-1) i ∞ ≤ -λ + (2 k+1 -4)τ + 2 k+4 -(k + 4) 2 -7 for i ∈ {0, 1}. Since m (k) 0 = m (k-1) 0 m (k-1) 1
, it follows that deg(m

(k) 0 ) = 2 k . By applying Lemma 5 we deduce that lg m (k) i ∞ = lg m (k-1) 0 m (k-1) 1 ∞ ≤ 2 k τ + 2 k+3 -2k -8 and lg m (k) 0 -m (k) 0 ∞ = lg m (k-1) 0 m (k-1) 1 -m (k-1) 0 m (k-1) 1 ∞ ≤ -λ + (2 k+2 -4)τ + 2 k+5 -(k + 5) 2 -7
as claimed.

To estimate the overall complexity, note that at the hth level of the tree we perform n/2 h multiplications of polynomials of degrees at most 2 h-1 for h = 2, . . . , k -1. We can assume that we perform all the computation with precision O(ℓ + nτ + lg n), and so the overall cost of the algorithm is We can compute an ℓ-approximation of the fan-out process in OB(n lg 2 n µ(ℓ + τ1 lg n + ρ n lg n)) provided that λ = ℓ + 2τ1 lg n + 300(ρ + 1)n lg n.

k n 2 k OB(2 k-1 lg 2 k-1 µ(ℓ + nτ + lg n)) = OB(n lg 2 n µ(ℓ + nτ + lg n)).
Proof: We keep assuming for simplicity that n = 2 k and proceed as in the proof of Lemma 17. Recall that |xi| ≤ 2 ρ for all the subscripts i, and so 2 ρ bounds the roots of all polynomials m (k) j . We can prove by induction, by using the bounds of Theorem 11 and the simplifications of Remark 12, that the precision of λ = ℓ + 2τ1 lg n + 300(ρ + 1)n lg n bits is sufficient.

At the hth step of the algorithm, for each h, we perform 2 h approximate polynomial divisions of polynomials of degree n 2 h using Theorem 11. We assume performing all the operations with the maximum precision, and bound the overall complexity by

lg n h=0 2 h OB( n 2 k (lg n 2 h ) 2 µ(ℓ + τ1 lg n + ρ n lg n))
= OB(n lg 2 n µ(ℓ + τ1 lg n + ρ n lg n))

BOUNDS ON THE COMPLEXITY OF BA-SIC ALGORITHMS

Lemma 19 (Multiplication of m polynomials). Let Pj ∈ ( C[x] of degree n and Pj ∞ ≤ 2 τ . Let Pj be λ-approximation of Pj with λ = ℓ + (4m -4)τ + (4m + 2 lg m -4) lg n + 32m, where 1 ≤ j ≤ m. We can compute j Pj such that

j Pj -j Pj ∞ ≤ 2 -ℓ in OB(m n lg m lg(m n) µ(λ)) or OB(m n (ℓ + mτ )). Moreover, lg j Pj ∞ ≤ mτ + (m - 1) lg n + 4m -lg m -4.
Proof: The algorithm is similar to the Fan-in process of Moenck-Borodin algorithm. Let p (0) j = Pj and compute recursively the polynomials p

(h+1) j = p (h) 2j p (h)
2j+1 , for 0 ≤ j ≤ 2 k-h , h = 0, . . . , k -2. Let m = 2 h . We prove the bounds on the infinite norm and the approximation using induction on h.

For h = 1, we compute the polynomials p Assume the claimed bounds for h -1, that is lg p

(h-1) j ∞ ≤ 2 h-1 τ + (2 h-1 -1) lg n + 4 • 2 h-1 -lg 2 h-1 -4 and lg p (h-1) j -p (h-1) j ∞ ≤ -λ + (4 • 2 h-1 -4)τ + (4 • 2 h-1 + 2 lg 2 h-1 -4) lg n + 32 • 2 h-1 for j = 0, 1. By ap- plying Lemma 5 for 2lg(K) ≤ lg(n), deduce the following bounds, lg p (h) j ∞ ≤ 2 h τ + (2 h -1) lg n + 4 • 2 h -lg 2 h -4,
which agrees with the claimed norm bound, and lg p

(h) j - p (h) j ∞ ≤ -λ + (4 • 2 h -4)τ + (4 • 2 h + 2 lg 2 h -4) lg n + 24 • 2 h + 2h -12
which is smaller than the claimed bound on the precision.

To estimate the overall complexity note that at each level, h, of the tree we have to perform m/2 h multiplications of polynomials of degrees at most 2 h-1 n. We can assume that we perform all the computations with the precision λ + (4m -4)τ + (4m + 2 lg m -4) lg n + 32m, or O(ℓ + mτ + m lg n), and so the overall Boolean cost of performing the algorithm is h m

2 h OB(2 h-1 n lg (2 h-1 n)µ(ℓ + nτ + m lg n)) = OB(m n lg m lg(m n) µ(ℓ + nτ + m lg n)) , which concludes the proof.
If the degrees of Pj vary as j varies, then we can apply a more pedantic analysis based on Huffman trees, see [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF].

The problem of computing (approximately) the sum of rational functions reduces to the problem on multiplying polynomials, which admits the same asymptotic complexity bounds. To estimate the overhead constants, we should also take into account the polynomial additions involved.

We have the following lemma.

Lemma 20 (Sum of rational functions). Suppose Pj ∈

( C[x] has degree n, Qj ∈ ( C[x] has a smaller degree, Pj ∞ ≤ 2 τ 2 , and Qj ∞ ≤ 2 τ 1 .
Assume λ-approximations of Pj by Pj and of Qj by Qj where λ = ℓ + τ1 + (4m -4)τ2 + (5m + 2 lg m -4) lg n + 32m and 1 ≤ j ≤ m.

Let Q P = j Q j P j . We can compute an ℓ-approximation of

Q/P , in OB(m n lg m lg(m n) µ(λ)) or OB(m n (ℓ + +τ1 + mτ2)). Moreover, lg Q ∞ ≤ τ2 + (m -1)(τ2 + lg n) + 5m - lg m -4 and lg P ∞ ≤ m τ2 + (m -1) lg n + 4m -lg m -4. Lemma 21 (Modular representation). Let F ∈ ( C[x] of degree 2mn and F ∞ ≤ 2 τ 1 . Let Pj ∈ ( C[x]
of degree n, for 1 ≤ j ≤ m. Moreover, 2 ρ be an upper bound on the magnitude of the roots of all Pj , for all j. Assume λapproximations of F by F and of Pj by Pj such that F -

F ∞ ≤ 2 -λ and Pj -Pj ∞ ≤ 2 -λ .
Furthermore assume that λ = ℓ + τ1 lg m + 60 n m(ρ + 3) lg(m n) + 60 lg m lg 2 (m + n) or λ = ℓ + O(τ1 lg m + m n ρ).

Then we can compute an ℓ-approximation Fj of

Fj = F mod Pj such that Fj -Fj ∞ ≤ 2 -ℓ in OB(m n lg n lg 2 (m + n) µ(ℓ + τ1 lg m + m n ρ))
or OB(m n (ℓ + τ1 + m n ρ)).

Moreover, lg Fj ∞ ≤ τ1 + (ρ + 1) m n + n + lg(m n).

Proof: First we perform the Fan-in process with polynomials Pj using the algorithm of Lemma 19. Assume that Pj ∞ ≤ 2 τ 2 and that we are given λ1-approximations. Following Lemma 19 we compute all the supermoduli, P (i) j so that lg P (i)

j -P (i) j ∞ ≤ -λ1+(4m-4)τ2+(4m+2 lg m-4) lg n+32m .
Remark 12 implies that τ2 ≤ 2n + nρ, and so lg

P (i) j - P (i) j ∞ ≤ -λ + (4m -4)(2n + lg n + nρ) + 2 lg m lg n + 32m = -λ + O(m n ρ) .
For computing ℓ-approximations of Fj = F mod Pj we mimic the procedure of the Fan-out process. This means that we apply repeatedly Theorem 11, which we can refine by following Remark 12. The bounds accumulate at each step, and so

lg Fj ∞ ≤ τ1 + h 3 n 2 h + n + h + 2 h nρ + 1 ≤ (mn -n)(ρ + 3) + m(m + n).
We assume that we are given λ2-approximation of F and all the supermoduli. For the required precision we have

lg Fj -Fj ∞ ≤ -λ2 + h 25(ρ + 2)(h + lg n)n2 h + 80(h + lg n) 2 + τ1 + 30 ≤ -λ2 + τ1 lg m + 25 n m(ρ + 2) lg(m n) + 40 lg m lg 2 (m + n).
To ensure an ℓ-approximation for Fj we require λ = ℓ + τ1 lg m + 60 n m(ρ + 3) lg(m n) + 60 lg m lg 2 (m + n) = ℓ + O(τ1 lg m + m n ρ) approximations of the input to ensure the validity of both the Fan-in and Fan-out process.

We assume that we perform all the computations with maximum accuracy. The complexity of computing the supermoduli is OB(m n lg m lg(m n) µ(ℓ+τ1 lg m+m n ρ)) or OB(m n (ℓ+ τ1 + m n ρ)).

For the complexity of the Fan-out process we proceed as follows. At each step, h, of the algorithm we perform 2 h approximate polynomial divisions of polynomials of degree m n 2 h using Theorem 11. The overall complexity is

lg m h=0 2 h OB( mn 2 h (lg mn 2 h ) 2 µ(ℓ + τ1 lg m + m n ρ)) = OB(m n lg n lg 2 (m+n) µ(ℓ+τ1 lg m+m n ρ)) or OB(m n (ℓ+ τ1 + m n ρ)).

LAGRANGE INTERPOLATION

Problem 2. Lagrange polynomial interpolation. Given the knot set (or vector) {xi} n-1 i=0 of n distinct points x0, . . . , xn-1 and the set (or vector) of values {yi} n-1 i=0 , compute a set (or vector) {aj} n-1 j=0 such that n-1 j=0 aj x j i = yi, i = 0, 1, . . . , n -1, that is, recover the coefficients of a polynomial A(X) = n-1 j=0 aj X j from its values at n distinct points x0, . . . , xn-1.

We follow the approach presented in [START_REF] Pan | Structured Matrices and Polynomials: Unified Superfast Algorithms[END_REF]Section 3.3], to which we also refer for a detailed presentation.

Lemma 22. Let |xi| ≤ 2 τ 1 , |yi| ≤ 2 τ 2 , and ∆i(x) = minj|xi-xj|, for all 0 ≤ i ≤ n-1. Assume λ-approximations of xi and yi, where λ = ℓ+68n(τ1+3) lg n+4nτ2-6 lg i ∆i(x)+50n+ 60 lg 3 n+20 or λ = ℓ+O(nτ1 lg n+nτ2 -lg i ∆i(x)+lg 3 n). Then we can compute an ℓ-approximation of the Lagrange polynomial interpolation in OB(n lg 2 n µ(λ)) or OB(n 2 τ1 + n 2 τ2 -n lg i ∆i(x))).

Proof: The input is given as λ-approximations, where λ is to be specified in the sequel. We track the loss of accuracy at each step of the algorithm.

1. Compute B(X) = i (X -xi) = X n + n-1 k=0 b k X k .
For this task we apply Lemma 17. In this case the infinite norm of B is bounded as follows, lg B ∞ ≤ nτ1 + 8n -2 lg n -4, and the computed approximation, B, is such that lg B -B ∞ ≤ -λ0 + (4n -4)τ1 + 16n + 20. As by-product we can compute the "supermoduli" j (x -xj) and then reuse them at stage L5.

Compute B ′ (X).

This operation increases the norm and the precision bounds by a factor of n in the worst case. That is lg B ′ ∞ ≤ nτ1 + 8n -lg n -4 and lg B ′ -B ′ ∞ ≤ -λ0 + (4n -4)τ1 + 16n + lg n + 20 = -λ1.

Evaluate B ′ at all points xi.

Perform this task using Lemma 21. This lemma implies that lg|B ′ (xi)| ≤ (n -1)(τ1 + 3) + n(n + 1). However, in this special case we can decrease the bound as follows,

|B ′ (xi)| ≤ j B ′ ∞ |xi| n-1 ≤ j 2 nτ 1 +8n-lg n-4 2 (n-1)τ 1 ≤ 2 (2n-1)τ 1 +8n-4 .
We achieve the accuracy lg|B ′ (xi)-B ′ ( xi)| ≤ -λ1 +(nτ1 + n -lg n -1) lg n + 60n(τ1 + 3) lg n + 60 lg 3 n = -λ2 .

Consider the rational functions Ai(X)

= A i,0 (X) A i,1 (X) = y i /B ′ (x i ) (X-x i )
. Deduce that lg Ai,1 ∞ ≤ τ1, and so the approximation bound matches that of xi.

To compute the relevant quantities of the numerator(s) we need a lower bound for B ′ (xi), for all i. We notice that

B ′ (X) = n i=1 j =i (X -xj). Thus B ′ (xi) = j =i (xi - xj) and so |B ′ (xi)| ≥ j =i ∆j (x), and lg Ai,0 ∞ ≤ τ2 - lg j =i ∆j(x) ≤ τ2 -lg j ∆j(x).
For computing an approximation of the denominator we rely on (complex) interval arithmetic. That is |Ai,0 -

Ai,0| ≤ wid([Ai,0]) = wid([yi/B ′ (xi)]).
We compute wid([yi/B ′ (xi)]) based on Prop. 3, and so lg wid

([1/B ′ (xi)]) ≤ -λ2 -4 lg j ∆j (x) + 2(2n -1)τ2 + 2n -8 + 3 = -λ3. Finally wid([yi/B ′ (xi)]) ≤ 2 τ 2 wid([1/B ′ (xi)]) + 2 -lg j ∆ j (x) wid([yi]) ≤ 2 -λ 3 +τ 2 -lg j ∆ j (x) ≤ 2 -λ 4
5. Compute the sum of the rational functions ie A 0 (X)

A 1 (X) = i A i,0 (X)
A i,1 (X) . Using Lemma 20 we get lg A0 ∞ ≤ τ2 -lg j ∆j (x)+(n-1)τ1 + 4n -lg n -4 and lg A1 ∞ ≤ nτ1 + 4n -lg n -4.

For the approximation we have that lg A0 -A0 ∞ ≤ -λ4 + τ2 -lg j ∆j (x) + (4n -4)τ1 + 32n. If we substitute the various values for λi we have lg A0 -A0 ∞ ≤ -λ + 68n(τ1 + 3) lg n + 4nτ2 -6 lg j ∆j (x) + 50n + 60 lg 3 n + 20.

The numerator, A0, is the required polynomial A(X). To achieve an ℓ-approximation of A(x) we assume that we perform all the computations using the maximum precision, that is ℓ + 68n(τ1 + 3) lg n + 4nτ2 -6 lg j ∆j(x) + 50n + 60 lg 3 n+20 or λ = ℓ+O(nτ1 lg n+nτ2 -lg j ∆j(x)+lg 3 n).

The overall complexity is OB(n lg 2 n µ(λ)) or OB(n 2 τ1 + n 2 τ2 -n lg j ∆j(x))).
Remark 23 (The hidden costs). In the previous lemma we have assumed bounds on the minimum distance between the xi's, which we denote by ∆i(x). The complexity results depend on this quantity, as it is very important in the computation of the number of bits that we need to certify the result to a desired accuracy. It is reasonable to assume that such bounds are part of the input.

However, how do we handle the case where such bounds are not known? As the precision required for the computations depends on these bounds we should be able to compute them, given the points xi.

We consider the numbers xi ∈ ( C as points on IR 2 and we compute their Voronoi diagram. This costs O(n lg n) operations, e.g. [START_REF] Preparata | Computational geometry. texts and monographs in computer science[END_REF]. Then for each point xi we find its closest in O(lg n) operations. Therefore, we can compute the quantities ∆i(x) in O(n lg n) operations. But what about the required precision? What are the required primitive operations for these computations? We only need to evaluate the signs of 3 × 3 determinants. The precision of Lemma 22 is sufficient for these operations.

SOLUTION OF A CAUCHY LINEAR SYS-

TEM OF EQUATIONS

Multiplication of a Cauchy matrix by a vector

We consider the problem of computing the matrix vector product Cv, where

C = C(s, t) = ( 1 s i -t j ) n-1 i,j
is a Cauchy matrix and v = (vi) n-1 i=0 . We refer the reader to [16, Problem 3.6.1] for further details of the algorithm.

Let |si| ≤ 2 τ 1 , |ti| ≤ 2 τ 2 , |vi| ≤ 2 τ 3 , and ∆i(t) = minj|ti -tj|, for all i.

The following quantities are also useful ∆j(s, t) = mini|sj -ti| and ∆(s, t) = minj ∆j(s, t).

Lemma 24. If the input is given as a λ-approximation, where λ = ℓ + 90n(τ1 + 3) lg n + 32(n -1)τ2 lg n + 30τ3 lg n -35 -24 lg ∆(s, t) -4 lg k ∆ k (t), then we can compute an ℓ-approximation of the vector Cv in OB(n

lg 2 n µ(λ)) or OB(n 2 τ1 + n 2 τ2 + nτ3 -n lg ∆(s, t) -n lg k ∆ k (t)).
Proof: The input is given as λ-approximations, where λ is to be specified in the sequel. We track the loss of accuracy at each step of the algorithm.

Consider the rational functions

v k x-t k = P k Q k .
For the numerators it holds deg(P k ) = 0, P k ∞ ≤ 2 τ 3 , and lg P k -P k ∞ ≤ -λ. For the denominators it holds

deg(Q k ) = 1, Q k ∞ ≤ 2 τ 2 , and lg Q k -P k ∞ ≤ -λ.

Compute the sum

P Q = k P k Q k .
For this computation we rely on Lemma 20.

For the numerator of the result we have deg(P ) ≤ n -1, lg P ∞ ≤ τ3 + (n -1)τ2 + 5n -lg n -4, and lg P -P ∞ ≤ -λ + τ3 + (4n -4)τ2 + 32n = -λ1.

For the denominator of the result we have deg

(Q) ≤ n, lg Q ∞ ≤ nτ2 + 4n -lg n -4, and lg Q -Q ∞ ≤ -λ + τ3 + (4n -4)τ2 + 32n = -λ1.
3. Compute P (si) and Q(si) for all i.

For this multipoint evaluation we use Lemma 21 (with m = n, n = 1, τ1 = lg P ∞, ρ = τ1).

We have lg|P (si)| ≤ (τ1 + 1)n + (n -1)τ2 + τ3 + 4n -4 and lg|P (si) -P ( si)| ≤ -λ1 + (τ3 + (n -1)τ2 + 4n -lg n -4) lg n + 60n(τ1 + 3) lg n + 60 lg 3 n = -λ2.

Similar bounds hold for Q(si).

Compute the fractions P (s i )

Q(s i ) . These are the elements of the result of the matrix-vector multiplication Cv.

For this task we need to perform n (complex) divisions. We use complex interval arithmetic to compute the loss of precision, as we did for deriving the bounds for Lagrange interpolation. To compute a lower bound for |Q(si)| we use Lemma 2, and so

|Q(si)| ≥ (∆i(s, t)) 6 2 -6 lg Q ∞-6 lg n 2 lg k ∆ k (t)-6 ≥ (∆i(s, t)) 6 2 lg k ∆ k (t) 2 -6nτ 2 -24n+20 ≥ 2 -ν Let |Q(si)| ≤ 2 T , where T = nτ1 + nτ + 2 + 4n -4. Using Prop. 3, wid[1/Q(si)] ≤ 2 4ν+2T +3 2 -λ 2 ≤ 2 -λ 3 and wid[P (si)/Q(si)] ≤ 2 |P (si)| 2 -λ 3 +2 |1/Q(si)| 2 -λ 2 ≤ 2 -λ 3 +T +2
which is also the accuracy of the result.

Putting together the various values of λi, ν, and T , we achieve an ℓ approximation by choosing λ = ℓ + 90n(τ1 + 3) lg n + 32(n -1)τ2 lg n + 30τ3 lg n -35 -24 lg ∆(s, t) -

4 lg k ∆ k (t).
We perform all the computations with maximum required accuracy, and so the overall complexity is OB(n lg 2 n µ(λ)) or OB(n 2 τ1 + n 2 τ2 + nτ3 -n lg ∆(s, t) -n lg k ∆ k (t)).

Trummer's problem

This is the important special case where s = t and the diagonal entries of the Cauchy matrix are set to zero. In this case we compute the matrix vector product by using the following formula,

(C v) n-1 i=0 = 2 P ′ (si) -vi Q ′′ (si) 2 Q ′ (si) n-1 i=0 = A0,i A1,i n-1 i=0 (9)
We refer the reader to [START_REF] Pan | Structured Matrices and Polynomials: Unified Superfast Algorithms[END_REF]Problem 3.6.3] for further details. Proof: First we compute bounds for the numerator of Eq. ( 9). Following the proof of Lemma 24 we have lg|P ′ (si)| ≤ (τ1 + 1)n + (n -1)τ2 + τ3 + 4n -4 + lg n and lg|P (si) -P ( si)| ≤ -λ2 +log n. Similarly for lg|Q ′′ (si)| ≤ (τ1 +1)n+(n-1)τ2 + τ3 + 4n -4 + 2 lg n and lg|Q ′′ (si) -Q ′′ ( si)| ≤ -λ2 + 2 log n.

For the first derivative we add to the logarithm of the norm a term lg n and for the second a term 2 lg n. Moreover, τ1 = τ2, as s = t.

Taking into account that |vi| ≤ 2 τ 3 we deduce that lg|A0,i| ≤ (τ1 + 1)n + (n -1)τ2 + 2τ3 + 4n -4 + 2 lg n and lg|A0,i -A0,i| ≤ -λ2 + τ3 + 2 log n + 2 .

Regarding the denominator we have that lg|Q ′ (si)| ≤ (τ1+ 1)n+(n-1)τ2 +τ3 +4n-4+lg n = T . In addition |Q ′ (si)| = j =i |si -sj| ≥ j =i ∆j (s) and so |1/A1,i| = |1/2Q ′ (si)| ≤ j =i (∆j(s)) -1 . Prop. 3 leads to lg wid([1/A1,i]) = lg wid[1/Q ′ (si)] ≤ -λ2 + lg n -4 lg j =i ∆j (s) + 2T + 3 where λ2 is defined at the (C3) step of the proof of Lemma 24. Putting all the pieces together, we have wid([

A 0,i A 1,i ]) ≤ 2 |A0,i| 2 -λ 2 +lg n-4 j =i ∆ j (s)+2T +3 +2 |1/A1,i| 2 -λ 2 +τ 3 +2 log n+2
and after many simplifications and overestimations lg

A 0,i A 1,i - A 0,i A 1,i ≤ wid([ A 0,i A 1,i ]) ≤ -λ + 70(τ1 + 3)n lg n + 4τ3 lg n -4 lg j ∆(s).
To achieve an ℓ-approximation we need the input to be a λ-approximation, where λ = ℓ + 70(τ1 + 3)n lg n + 4τ3 lg n -4 lg j ∆(s), and we perform all the computations with this number of bits. The overall complexity is OB(n lg 2 n µ(λ)).

Solving a (Cauchy) linear system

We consider the following problem.

Problem 3 (Cauchy linear system of equations). Solve a non-singular Cauchy linear system of n equations, C(s, t) v = r for an unknown vector v and 3 given vectors r, s, and t. Theorem 26. Let |si| ≤ 2 τ 1 , |ti| ≤ 2 τ 2 , |ri| ≤ 2 τ 3 , and ∆i(s) = minj |si -sj|, ∆i(t) = minj |ti -tj|, for all i. Let also ∆(s, t) = mini,j |si -tj|.

If the input is given λ-approximations, where λ = ℓ + 630(τ1+τ2)n lg n+32τ3 lg n-35-35 lg n lg ∆j (s)-5 lg ∆j (t)-25 lg ∆(s, t), then an ℓ-approximation solution to Problem 3 could be obtained in OB(n lg 2 n µ(λ)), or OB(nℓ + n 2 (τ1 + τ2) + nτ3 -n lg ∆j (s) -n lg ∆j (t) -n lg ∆(s, t)).

Proof: Following [16, Eq. 3.6.10] then inverse of C(s, t) could be obtained as follows

C -1 (s, t) = diag(ps(ti)/p ′ t (ti)) n-1 i=0 C(t, s) diag(pt(si)/p ′ s (si)) n-1 i=0 = D1 • C(t, s) • D2 , (10)
where pt(X) = i (X -ti) and ps(X) = i (X -si).

Using this formula we can solve the linear system as

v = C -1 (s, t) r = D1 • C(t, s) • D2 r
We apply Lemma 17 to analyze the computation of the polynomials ps(x) and pt(x). Then we compute the two diagonal matrices D1 and D2 by applying the Moenck-Borodin algorithm for multipoint evaluation (see Section 4). At first we perform n ops to compute the vector r1 = D2 r. Then we use Lemma 24 to obtain the vector r2 = C(t, s)r1. Finally we multiply the matrix D1 by r2 to obtain the vector v. We track the loss of precision at each of these three steps.

1. Computation of the matrices D1 and D2.

For this task we need to compute ps(ti), p ′ s (si), pt(si), and p ′ t (ti). It holds lg ps ∞ ≤ nτ1 + 4n -lg n -4 and lg ps -ps ∞ ≤ -λ + (4n -4)τ1 + 32n, using Lemma 19. By applying Lemma 21 we get lg ps(ti) ∞ ≤ nτ1 + nτ2 + 5n -3. and lg ps(ti)-ps( ti) ∞ ≤ -λ+(n lg n+4 n-4)τ1+60 nτ2 lg n-4 lg n + 184 n lg n + 32 n -(lg n) 2 . However, to simplify the calculations we consider the inferior bounds lg ps(ti) ∞ ≤ 7n(τ1 + τ2). and lg ps(ti)ps( ti) ∞ ≤ -λ + 300(τ1 + τ2)n lg n for all the involved quantities. We also need lower bounds for |p ′ t (ti)| and

|p ′ s (si)|. It holds |p ′ t (ti)| ≥ j =i |ti -tj| ≥ j =i ∆j (t). Similarly |p ′ s (si)| ≥ j =i |si -sj | ≥ j =i ∆j (s).
This leads to the bounds: ps(t i ) s) . By combining the previous bounds with the complex interval arithmetic of Prop. 1 we obtain the following estimation for the approximation:

p ′ t (t i ) ≤ 2 7n(τ 1 +τ 2 )-lg j =i ∆ j (t) and p t (s i ) p ′ s (s i ) ≤ 2 7n(τ 1 +τ 2 )-lg j =i ∆ j (
lg pt(si) p ′ s (si) - pt(si) p ′ s (si) ≤ -λ + 315(τ1 + τ2)n lg n -4 lg j =i ∆j(s) lg ps(ti) p ′ t (ti) - ps(ti) p ′ t (ti) ≤ -λ + 315(τ1 + τ2)n lg n -4 lg j =i ∆j (t) 2. r1 = D2 r .
This computation increases the bounds by a factor of τ3.

To be more specific, for the elements of r1, r1,i we have |r1,i| ≤ 2 7n(τ 1 +τ 2 )+τ 3 -lg j =i ∆ j (s) and lg|r1,i -ri,1| ≤ -λ + 316(τ1 + τ2)n lg n + τ3 -4 lg j =i ∆j (s)

3. r2 = C(t, s) r1.

For this computation we need to apply Lemma 24. We obtain To achieve an ℓ-approximation of the output we should require a λ-approximation of the input, where λ = ℓ+630(τ1 + τ2)n lg n + 32τ3 lg n -35 -35 lg n lg ∆j(s) -5 lg ∆j (t) -25 lg ∆(s, t).

As in all the previous sections we perform all the computations using the maximum precision. The overall complexity of the algorithm is OB(n lg 2 n µ(λ)), or OB(nℓ + n 2 (τ1 + τ2) + nτ3 -n lg ∆j(s) -n lg ∆j(t) -n lg ∆(s, t)).

|T -1 i,k | ≤ |T -1 i,1 |+|T -1 i,k | = |qi| ≤ 2 n+lg n+nρ φ ν,k ∞ ≤ 2 n+lg n+nρ . (11
) Let P0 and P1 denote the two polynomials defined earlier in the proof such that maxi,j{|T -1 i,j |} ≤ ||P 2 0 P1||∞. Then Eq. ( 11) implies the inequality

lg P 2 0 P1 ∞ ≤ n + lg n + nρ = N . (12)
Recall that n+1 = 2 k by assumption and that the polynomial P 2 0 P1 has degree 2n -2. Write h = ⌈lg(2n -2)⌉. Under the assumption that the input elements are known up to the precision of λ bits, we compute a polynomial P 2 0 P1 such that

lg P 2 0 P1 -P 2 0 P1 ∞ ≤ ≤ -λ + (2 lg n + 8)τ + 2 lg n(4 lg n + 27) + 8(lg n -1)N ≤ -λ + (2k + 8)τ + 2k(4k + 27) + 8(k -1)N . (13)
We proceed by induction. Let the base case be n + 1 = 4; then k = 2. We need to invert the following matrix

T =        1 0 0 0 t2 1 0 0 t1 t2 1 0 t0 t1 t2 1        = T0 0 T1 T0
, where

T0 = 1 0 t2 1 , T -1 0 = 1 0 -t2 1 , T1 = t1 t2 t0 t1 ,
and |ti| ≤ 2 τ . The associated polynomials are P0(x) = 1-t2x, for T -1 0 , and P1(x) = t2+t1x+t0x 2 , for T1. Therefore

P1P0 = (1 -t2x)(t2 + t1x + t0x 2 ) = t2 + (t1 -t 2
2 )x + (t0 -t1t2)x 2 -t0t2x 3 . The subvector v = (t1 -t 2 2 , t0 -t1t2) T of the coefficient vector of the polynomial product P1P0 is the first column of the matrix product T1T -1 0 . Furthermore the vector -T -1 0 v = (t 2 2 -t1, 2t1t2 -t 3 2 -t0) T is a subvector of the coefficient vector of the polynomial product -P0v = -(1 -t2x)(t1 -t 2 2 + (t0 -t1t2)x) where v denotes the polynomial t1 -t 2 2 + (t0 -t1t2)x with the coefficient vector v. As we proved earlier, the subvector is the first column of the matrix

-T -1 0 T1 T -1 0 = t2 2 -t1 -t2 -t2 3 + 2 t1t2 -t0 t2 2 -t1 .
(This is not a subvector of the coefficient vector of the polynomial

P 2 0 P1 = (1 -t2x) 2 (t2 + t1x + t0x 2 ) = t2 + (t1 -2t 2 2 )x + (t0 -2t1t2 + t 3 2 )x 2 + (t 2 2 t1 -2t2t0)x 3 + t 2 2 t0x
4 because multiplication of polynomials and the truncation of their coefficient vectors are not commutative operations, but as we observed, we still can reduce our study of the product T -1 0 T1T -1 0 to the study of p 2 0 p1.) We perform the multiplications by using the algorithm of Lemma 5 and the bounds from Cor. 8, to obtain a polynomial P 2 0 P1 such that lg P 2 0 P1-P 2 0 P1 ∞ ≤ -λ+10τ +15 lg 1+40 ≤ -λ+12τ +140+8N , where the last inequality is obtained by substituting k = 2 in Eq. [START_REF] Kobel | Fast approximate polynomial multipoint evaluation and applications[END_REF].

It remains to prove the induction. Assume that the claimed bounds are true for n + 1 = 2 k and extend them to n + 1 = 2 k+1 . In our case P0 is a polynomial of degree 2 k-1 -1. By induction hypothesis we know the coefficient of P0 within 2 -ℓ for ℓ defined by ( 13), and we can apply Eq. ( 12) to obtain P0 ∞ ≤ 2 N .

The polynomial P1 has degree 2 k -2, P1 ∞ ≤ 2 τ . Its coefficients are the entries of the matrix T and we know them within 2 -λ .

We apply Cor. 8 for d = 2 k-1 -1, τ0 = N , τ1 = τ , and These inequalities imply that all our computations require a precision bound of at most λ ′ = ℓ + (2k + 8)τ + 2k(4k + 27) + 8(k -1)N . To simplify the formula, we substitute k = lg(n) rather than k = lg(n + 1) and then rewrite λ ′ as follows, ℓ + (2 lg n + 8)τ + 2 lg n(4 lg n + 27) + 8(lg n -1)(n + lg n + nρ + 1) bits, which we simplify to the bound of ℓ+10τ lg n+70 lg 2 n+8(ρ+1)n lg n or O(ℓ+(τ +lg n+nρ) lg n) bits.

-ν = -λ+(2(k-1)+8)τ +2(k-1)(4(k-1)+27)+8(k-2)N , substitute k ≥ 3,
To estimate the overall complexity of approximating the first column of the inverse matrix T -1 , we notice that we perform k steps overall, for k = lg(n + 1), but we will keep writing k = lg(n) to simplify the notation.

At each step we perform two multiplications of polynomials of degrees at most 2 k with the coefficients having absolute values less than 2 O(nτ ) , and we use the precision of ℓ+10τ lg n+70 lg 2 n+8(ρ+1)n lg n = O(ℓ+(τ +lg n+nρ) lg n) bits. All these bounds together imply that

A NORMALIZATION OF KIRRINNIS' RE-SULTS

The following results express the bounds of Kirrinnis [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF] for polynomials of arbitrary norms; for this an appropriate scaling is applied. The complexity bounds depends on polynomial multiplication algorithms that rely on Kronecker substitution and not on FFT as our results that we presented earlier. The bounds are asymptotically the same as ours, up to logarithmic factors.

The following is from [11, Theorem 3. We can compute j Pj such that j Pjj Pj ∞ ≤ 2 -ℓ in OB(µ(n • lg n • (ℓ + nτ + j lg nj ))) or OB(n(ℓ + nτ )).

Proof: Let pj = 2 -τ -lg n j +n j Pj . Then Pj ∞ ≤ 2 τ ⇒ Pj 1 ≤ 2 τ +lg n j ⇒ 2 -τ -lg n j +n j Pj 1 ≤ 2 n j ⇒ pj 1 ≤ 2 n j . Moreover, Pj -Pj ∞ ≤ 2 -λ ⇒ Pj -Pj 1 ≤ 2 -λ+lg n j ⇒ pj -pj 1 ≤ 2 -τ -lg n j +n j 2 -λ+lg n j = 2 -λ-τ +n j .

The algorithm of Kirrinnis [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF]Theorem 3.7 and Algorithm 5.1] computes and approximation of j pj such that j pjj pj 1 ≤ 2 -s+τ +3n , which implies j 2 -τ -lg n j +n j Pj - If we want the error to be less than 2 -ℓ , then we need to choose initial precision λ ≥ ℓn(τ + 2) + n lg n. The total cost is OB(µ(n • lg n • (ℓ + nτ + j lg nj))) or OB(n(ℓ + nτ )).

The algorithm for computing the sum of rational functions [11, Theorem 3.8 and Algorithm 5.2] relies on Lemma 27 and the bounds are similar, so we do not elaborate further.

The following lemma relies on [11, Theorem 3.9 and Algorithm 5.3].

Lemma 28 (Modular representation). Let F ∈ ( C[x] of degree m and F ∞ ≤ 2 τ 1 . Let Pj ∈ ( C[x] of degree nj, for 1 ≤ j ≤ ν, such that j nj = n and m ≥ n. Moreover, 2 ρ be an upper bound on the magnitude of the roots of all Pj . Let F , resp. Pj, be an approximation of F , resp. Pj , such that F -F ∞ ≤ 2 -λ , resp. Pj -Pj ∞ ≤ 2 -λ .

If λ = ℓ + τ1 + 2(n + m)ρ then we compute approximations Fj of Fj = F mod Pj , such that Fj -Fj ∞ ≤ 2 -ℓ in OB(µ((m + n lg n)(ℓ + τ1 + (m + n)ρ))).

Proof: The polynomials Pj (2 ρ x) have all their roots inside the unit disc. We make them monic and then denote them pj. It holds:

Pj -Pj ∞ ≤ 2 -λ ⇒ Pj (2 ρ x)-Pj(2 ρ x) 1 ≤ 2 -λ+n j ρ+lg n j ⇒ pj -pj 1 ≤ 2 -λ+n j ρ+lg n j +1 .
We apply the same transformation to F . Let f = 2 -τ 1 -mρ-lg m F (2 ρ x), then [START_REF] Kirrinnis | Partial fraction decomposition in C z and simultaneous Newton iteration for factorization in C[z][END_REF]Theorem 3.9] choosing λ = ℓ + τ1 + 2(n+m)ρ to guarantee f mod pj -fj 1 = fj -fj 1 ≤ 2 -ℓ . The complexity of the procedure is OB(µ((m + n lg n)(ℓ + τ1 + (m + n)ρ))).

F ∞ ≤ 2 τ 1 ⇒ F 1 ≤ 2 τ 1 +lg m ⇒ F (2 ρ x) 1 ≤ 2 τ 1 +mρ+lg m ⇒ 2 -τ 1 -mρ+lg m F (2 ρ x) 1 ≤ 1 ⇒ f 1 ≤ 1 .

Now we can use

Remark 29. In the case where m = n the latter bound becomes OB(n(ℓ + τ1 + nρ)).

✟ ✟ ✟ ✟ ✟ ✟ ✟ ✟ ❍ ❍ ❍ ❍ ❍ ❍ ❍ ❍ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆ ✁ ✁ ✁ ✁
= v(x) mod m (h) j ✟ ✟ ✟ ✟ ✟ ✟ ✟ ✟ ❍ ❍ ❍ ❍ ❍ ❍ ❍ ❍ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ❅ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆ ✁ ✁ ✁ ✁ ❆ ❆ ❆ ❆

  x) mod m (h) j for all j and h. Now, correctness of the output p(xi) = r

  and 3.2), whose nodes are the "supermoduli" m (h) j at the fanin stage 1, but turn into the remainders r (h) j at the fan-out stage 2.

∞ ≤ 2 ∞

 2 and that |xi| ≤ 2 τ for a positive integer τ . At the cost OB(n lg 2 n µ(ℓ + nτ + lg n)) the Fan-in process of the Moenck-Borodin algorithm approximates the "supermoduli" m -ℓ for all i and j. Moreover, lg m (i) j ∞ ≤ nτ + 8n -2 lg n -8 for all i and j.Proof: Assume that n = 2 k . The proof is by induction on k. Write m(0) i = x -xi and m (0) i = x -xi.Wlog we provide the estimates just in the case where j = 0.Consider the case where k = 1. Apply Lemma 5 for A = m ≤ -λ + 4τ + 14.2. This proves the induction basis.

Lemma 18 (

 18 Complexity of Fan-out process). Let v(x) ∈ ( C[x] of degree n -1 and v ∞ ≤ 2 τ 1 , and let v be a λapproximation. Let m (k) j be the supermoduli of the fan-in process and m (k) j their λ-approximations.

∞∞

  ≤ -λ. Apply Lemma 5 (and Remark 6) for K = n and τ1 = τ2 = τ to deduce that lg p (1) 0 ∞ ≤ 2τ + lg n + 3 , which agrees with our formula, and lg p ≤ -λ+4τ +5.1 lg n+4 ≤ -λ+4τ +6 lg n+64 , where the right hand-side represents the claimed bound for k = 1.

Corollary 25 .

 25 Using the notation of Lemma 24, we can solve Trummer's problem in OB(n lg 2 n µ(λ)), where λ = ℓ + 70(τ1 + 3)n lg n + 4τ3 lg n -4 lg j ∆(s).

  + 8τ0 + 2τ1 + 14k + 42 ≤ -λ + (2k + 8)τ + 8k 2 -2k + 104 + 8(k -1)N .

  lg (n+1) k=1 k • 2 • OB(2 k • lg 2 k • µ(ℓ + (τ + lg n + nρ) lg n)) = = OB(n lg 2 (n) µ(ℓ + (τ + lg n + nρ) lg n)) ,which concludes the proof.

  7 and Algorithm 5.1]. Lemma 27 (Multiplication of polynomials). Let Pj ∈ ( C[x] of degree nj and Pj ∞ ≤ 2 τ and Pj be an approximation of Pj such that Pj -Pj ∞ ≤ 2 -λ , with λ = ℓn(τ + 2) + n lg n., where 1 ≤ j ≤ m, n1 ≤ • • • ≤ nm, and n j = n.

j 2 - 1 ≤ 2 1 ≤ 2

 21212 τ -lg n j +n j Pj -λ+τ +3n ⇒ 2 -nτ -j lg n j +n j -λ+n(τ +2)+n lg n .

Figure 1 .

 1 Figure 1. Fan-in process, m (h+1) j = m (h) 2j m (h) 2j+1

Figure 2 .

 2 Figure 2. Fan-out process, r (h) j = r (h+1) ⌊j/2⌋ mod m (h) j

Acknowledgments. VP is supported by NSF Grant CCF-1116736 and PSC CUNY Awards 64512-0042 and 65792-0043. ET is partially supported by GeoLMI (ANR 2011 BS03 011 06), HPAC (ANR ANR-11-BS02-013) and an FP7 Marie Curie Career Integration Grant.

Appendix

Lemma 10. Let n + 1 = 2 k for a positive integer k and let T be a lower triangular Toeplitz (n + 1) × (n + 1) matrix of. Eq. ( 5), having ones on the diagonal. Let its subdiagonal entries be complex numbers of magnitude at most 2 τ known up to a precision 2 -λ . Let 2 ρ be an upper bounds on the magnitude of the roots of the univariate polynomial t(x)

Furthermore, to compute the entries of T -1 up to the precision of ℓ bits, that is to compute a matrix

Proof (of Lemma 10): We will prove the claimed estimates by reducing the inversion to recursive multiplication of polynomials defined by equation ( 6) and Lemma 7.

Consider the n+1 2 × n+1 2 Toeplitz matrices T -1 0 and T1 of equation ( 6). The Toeplitz matrix T -1 0 is triangular, and so its first column, p = (pi)

, with p0 = 1, defines this matrix and the polynomial p(x) = (n-1)/2 i=0 pix i of degree (n -1)/2. Likewise the vector t = (tn-i) n i=1 (made up of two overlapping vectors, that is, the reversed first row, (tn-1, tn-2, . . . , t (n-1)/2 ) of the matrix T1 and its first column, (t (n-1)/2 , t (n-3)/2 , . . . , t0) T ) defines this matrix and the polynomial t(x) = n i=1 tn-ix i-1 of degree n -1. The first column of the n+1 2 × n+1 2 Toeplitz matrix T1T -1 0 is a subvector, v, of dimension (n + 1)/2 of the coefficient vector of the polynomial product t(x)p(x), having degree 3(n -1)/2. Likewise the first column of the n+1

is the vector q = -T -1 0 v, which is a subvector of dimension (n + 1)/2 of the coefficient vector of the polynomial product p(x)v(x) of degree n -1, where the polynomial v(x) of degree (n -1)/2 is defined by its coefficient vector v. In sum the vector q is the coefficient vector of a polynomial q(x) obtained by two successive multiplications of polynomials, each followed by the truncation of the coefficient vectors. Namely we first compute the polynomial t(x)p(x), then truncate it to obtain the polynomial v(x), then compute the polynomial -p(x)v(x), and finally truncate it to obtain the polynomial q(x).

The truncation can only decrease the degree of a polynomial and the maximum length of its coefficients, and so we can bound the precision and the cost of computing the matrix product -T -1 0 T1T -1 0 by the bounds on the precision and the cost of computing the polynomial product P 2 0 P1, estimated in Corollary 8 for P0 = p(x), P1 = t(x), and d = (n -1)/2.

First we prove the upper bound on the elements of T -1 . Consider the division φ i,k (x) = t(x)q(x) + r(x), where t(x) is the univariate polynomial of degree n associated with

By abusing notation we also write φ i,k and q to denote the coefficient vectors of these polynomials. Using Eq. ( 5) we can compute the elements of q from the equation q = T -1 φ i,k . In this way |qi| = |sign(T -1 i,1 )T -1 i,1 + sign(T -1 i,k )T -1 i,k | = |T -1 i,1 | + |T -1 i,k |. Using Lemma 9 we obtain the inequality |qi| ≤ q ∞ ≤ 2 n+lg n+nρ φ i,k ∞, which in turn implies


A Amon

R Delannay

A Valance

Shear Instabilities in Bidisperse Granular Flows

Keywords: granular flow ; shear instabilities PACS: 45, 70, Mg

We study experimentally the stability of the interface between two shearing granular flows. Both flows consist of monosized spherical glass beads but they differ in particle size. When running on an inclined rough plane, the two flows reach different stationary velocities. Hence, if their lateral sides get into contact, a shear zone appears at the interface. We show that in a limited region of the phase space (inclination angle/flow rate), the interface is unstable and develops a wavy pattern.

INTRODUCTION

Instabilities in granular flows have not been explored deeply yet [START_REF] Aranson | Patterns and collective behavior in granular media : Theoretical concepts[END_REF]. In the case of shear instability (analogous to the Kelvin-Helmholtz instability in hydrodynamics), first experimental observations were reported by Goldfarb et al. [START_REF] Goldfarb | Shear instabilities in granular flows[END_REF]. However, such kind of instability is not fully characterized and the understanding of its mechanism is still open.

Goldfarb et al. [START_REF] Goldfarb | Shear instabilities in granular flows[END_REF] studied shear flows on an inclined chute. Two hoppers with different apertures allowed to obtain two flows with different rates. They used for both flows the same type of irregular sand, sieved or not. The rigid plane at the bottom of the channel was smooth and consequently the flows were accelerated along the chute. They observed instabilities at the interface between the two shearing flows only at small angles of inclination. When using spherical glass beads instead of sand, they did not observe any instability. The authors proposed a dimensionless number to characterize the competition between elongation due to acceleration and shear. They argued from their experiments that instabilities are favored when there is no elongation.

These observations motivated us to design a new experiment allowing to obtain a shear between two steady granular flows. In order to achieve such a configuration, we took advantage of the properties of monodisperse granular flows on bumpy inclined plane [START_REF] Gdr Midi | On dense granular flows[END_REF][START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF]. We created two granular streams of same height but differing in grain sizes. Experiencing a differential basal friction, they flow at different velocity producing a shear zone at the surface of contact. We investigated the behavior of the interface between these two streams according to the flow height and chute inclination. We showed in particular that the interface becomes unstable in a limited region of the phase space.

EXPERIMENTAL SETUP

Figure 1 shows a sketch of our experimental setup. The side walls of the chute are made of plexiglass. Two different chute configurations were used: a wide chute of 28.5 cm width and a narrow one of 13.2 cm width. The chute is 130 cm long and the base consists of 500 µm spherical beads glued on a rigid plane. At the upper end of the chute, two independent reservoirs with variable outlets allow to trigger two granular streams. For our flow experiment, we used two kinds of spherical glass beads : white fine beads of 250 µm diameter and black large ones of 500 µm diameter. Control parameters are the inclination angle of the plane with horizontal, θ , controlled with a precision of 0.1 • , and the thickness of the flowing layer, h, tuned with the outlet of the reservoirs. Two laser sheets hit the flow surface with small incident angles at two different positions along the chute. The displacement of their traces on the flow surface gives a measurement of the thickness of the flowing layer. A camera is placed above the flow to record the temporal evolution of the interface between the two granular streams.

FLOW FEATURES

We first recall briefly the different regimes observed with monodisperse granular assemblies flowing on rough inclines [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF]. For a given angle θ of inclination, granular layers thiner than a minimum thickness, h stop (θ ), comes into arrest. For thicknesses greater than h stop , flows can be sustained: they are steady and fully developed for moderate inclination angle and thickness, but become accelerated for larger values of these parameters [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF].

The ratio between the diameter of the flowing beads and that of the glued beads determines the effective friction experienced by the flow and consequently controls the flow velocity when steady uniform flow is obtained. For a given bumpiness of the plane, there exists a critical bead diameter for which the flow velocity is minimum, that is when the diameter of the flowing beads is approximatively half of the diameter of the glued beads [START_REF] Goujon | Monodisperse dry granular flows on inclined planes: Role of roughness[END_REF]. We determined the curves h stop (θ ) for monodisperse flows and mixtures using the two grain species. At a given fixed angle, we released the material placed in the reservoirs and waited for the flow to stop. The thickness of the remaining layer is h stop . Figure 2 shows the data obtained for the wide chute. Note that h stop has been normalized by the diameter of the large beads, d l =500 µm. Data corresponding to monodisperse flows can be well approximated by the following phenomenological law [START_REF] Gdr Midi | On dense granular flows[END_REF][START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF] :

h stop (θ ) d = γ tan θ 2 -tan θ tan θ -tan θ 1 ( 1
)
where d is the diameter of the flowing beads, and γ, θ 1 , θ 2 are fit parameters whose values are given in table 1 for both bead diameters and chute widths. For a given grain size, the values obtained are approximately the same for both chute widths. Nevertheless, we noted a variability along the plane, with a systematic tendency for h stop to decrease downstream along the slope. For a given inclination angle, h stop is higher for the fine beads than for the large ones. This is in agreement with the fact that 250 µm beads experience a larger effective friction on a bumpy plane made with 500 µm glued beads than the 500 µm particles [START_REF] Goujon | Monodisperse dry granular flows on inclined planes: Role of roughness[END_REF]. Data corresponding to mixtures collapse on the h stop curve obtained for the fine beads. We observed that the deposit obtained with mixtures is almost exclusively composed of fine grains. This is in agreement with the known behavior of polydisperse flows [START_REF] Goujon | Bidisperse granular avalanches on inclined planes : A rich varity of behaviors[END_REF] : segregation causes the large grains to move towards the front and the surface of the flow, and when the flow stops, the remaining layer is mostly composed of fine grains. For steady flows of velocity v, the Froude number v/ √ gh is expected to vary as [START_REF] Gdr Midi | On dense granular flows[END_REF] :

v √ gh = α + β h h stop , (2)
where α and β are parameters which depend on material properties of the beads but not on the grain diameter. For glass beads flowing over a rough plane made of glass beads, α = 0 and β = 0.136 [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF], whatever the bead diameter. Figure 3 shows data obtained with the wide chute. The slopes deduced from a linear fit are shown in the last column of table 1. Taking into account the dispersion of the measurements, the slope values obtained for the 500 µm beads fairly agree with the expected value. For the fine beads, there is a deviation from the expected value. Besides, the values of the slope obtained for both chute widths are different. We checked the granulometry of our beads and found that the fine beads have a much larger size distribution than the large ones. This could explain the discrepancy observed for large and small glass beads. By using the difference of friction experienced by two types of beads on a bumpy plane, we can obtain for a given chute inclination steady flows with same height but different velocity. It is then possible to produce a shear between two granular streams running side by side on the same bumpy plane. In the next section, we present succinctly the behavior of the contact interface between both streams.

RESULTS

Figure 4 shows a snapshot of unstable shear flows : the interface between the black fine grains and the white large ones develops a wavy pattern. Analysis of successive images of the flow allows to extract the interface position and to characterize its temporal evolution. The wavelength of the pattern typically ranges from 10 to 25 cm and the frequency from 5 to 10 Hz. The instability of the interface is not observed in the whole parameter space. Figure 5 shows the parameter values for which flows were carried out. We display data for large (fig. 5(a)) and narrow (fig. 5(b)) chute widths. Filled symbols correspond to parameters for which instability of the interface was observed, while open symbols are associated with stable flows. The curves h stop (θ ) for fine and large grains are presented as well.

One can note that for small values of h, interface instability is rarely observed. A natural question arises: what is the criterion for the instability? As the flow rule is different for both streams (see eq. two flows can vanish and change sign according to the values of the control parameters h and θ . In figure 5, we colored in grey the region for which |v lv s | > √ gh. Experimental data suggest that a criterion based on the velocity difference may be plausible.

CONCLUSION

In this article, we present a simple granular experiment where two steady streams of grains experience a shear at their surface of contact. Each stream is composed of monosized spherical glass beads. Particles from the two streams differ only in size. We show that the interface between the two flows can become unstable and develop a wavy pattern. Contrary to experiments of Goldfarb et al., the instability is observed with regular spherical particles. In addition, we provided a phase diagram in the parameter space (h, θ ) identifying the stable and unstable regions of the flow, and proposed a criterion for the instability onset based on the difference in flow velocity.

FIGURE 1 .

 1 FIGURE 1. Sketch of the experimental set-up. θ is the inclination of the chute with respect to the horizontal and g is the acceleration of gravity.

FIGURE 2 .

 2 FIGURE 2. h stop /d l versus θ for different mixtures of different proportions of fine and large grains. The solid and dashed lines correspond to best fits using eq. (1).

FIGURE 3 .

 3 FIGURE 3. Froude number v/ √ gh as a function of h/h stop for monodisperse flows. Filled symbols : 500 µm beads; open symbols : 250 µm beads. Solid and dotted lines are linear fits.

FIGURE 4 .

 4 FIGURE 4. Top view of the two streams.

  2 and table 1), it follows that the velocity difference v lv s between the

FIGURE 5 .

 5 FIGURE 5. Phase diagram in the parameter space (h/d l , θ ) for both chute configurations: wide (a) and narrow (b) chute. Filled symbols: unstable flow; open symbols: stables flows. Solid and dotted lines stand for the curves h stop (θ ) respectively for fine and large grains. The domains in grey correspond to regions where the velocity difference |v lv s | is greater than √ gh.

TABLE 1 .

 1 Fit parameters.

	d	channel γ	θ 1	θ 2	β
	250 µm	wide narrow	.8 .8	23 • 22 •	42 • 39 •	.20 .31
	500 µm	wide narrow	2.6 2.5	21 • 19 •	30 • 28 •	.13 .11

ACKNOWLEDGMENTS

We are grateful to several students who contributed to this work : R.M Guillermic, I. Hussein and S. Pasumarthy. We also thank A. Faisant for his advices concerning the design of the experimental set-up, S. Bourles and P. Chasle for technical assistance.


Matthieu Denoual

Patrick Attia

Standard VHDL Modeling and Top-Simulation for the Development of an Integrated Smart-Bolometer