Matías R Bender

Jean-Charles Faugère

Ludovic Perret

Elias Tsigaridas

A nearly optimal algorithm to decompose binary forms

Keywords: Decomposition of binary forms, Tensor decomposition, Symmetric tensor, Symmetric tensor rank, Polynomial Waring's problem, Waring rank, Hankel matrix, Algebraic degree, Canonical form

Symmetric tensor decomposition is an important problem with applications in several areas, for example signal processing, statistics, data analysis and computational neuroscience. It is equivalent to Waring's problem for homogeneous polynomials, that is to write a homogeneous polynomial in n variables of degree D as a sum of D-th powers of linear forms, using the minimal number of summands. This minimal number is called the rank of the polynomial/tensor. We focus on decomposing binary forms, a problem that corresponds to the decomposition of symmetric tensors of dimension 2 and order D, that is, symmetric tensors of order D over the vector space K 2 . Under this formulation, the problem finds its roots in invariant theory where the decompositions are related to canonical forms.

We introduce a superfast algorithm that exploits results from structured linear algebra. It achieves a softly linear arithmetic complexity bound. To the best of our knowledge, the previously known algorithms have at least quadratic complexity bounds. Our algorithm computes a symbolic decomposition in O(M(D) log(D)) arithmetic operations, where M(D) is the complexity of multiplying two polynomials of degree D. It is deterministic when the decomposition is unique. When the decomposition is not unique, it is randomized. We also present a Monte Carlo 1 variant as well as a modification to make it a Las Vegas one.

 to bound the size of the representation of the coefficients involved in the decomposition and we bound the algebraic degree of the problem by min(rank, D-rank +1).

We show that this bound can be tight. When the input polynomial has integer coefficients, our algorithm performs, up to poly-logarithmic factors, O B (D + D 4 + D 3 τ) bit operations, where τ is the maximum bitsize of the coefficients and 2 -is the relative error of the terms in the decomposition.

Introduction

The problem of decomposing a symmetric tensor consists in writing it as the sum of rank-1 symmetric tensors, using the minimal number of summands. This minimal number is known as the rank of the symmetric tensor 1 . The symmetric tensors of rank-1 correspond to, roughly speaking, the D-th outer-product of a vector. The decomposition of symmetric tensor is a common problem which appears in divers areas such as signal processing, statistics, data mining, computational neuroscience, computer vision, psychometrics, chemometrics, among others. For a modern introduction to the theory of tensor, their decompositions and applications we refer to e.g., [START_REF] Comon | Tensors: a brief introduction[END_REF]; [START_REF] Landsberg | Tensors: geometry and applications[END_REF].

There is an equivalence between decomposing symmetric tensors and solving Waring's problem for homogeneous polynomials, e.g., [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF]; [START_REF] Helmke | Waring's problem for binary forms[END_REF]. Given a symmetric tensor of dimension n and order D, that is a symmetric tensor of order D over the vector space Email addresses: matias.bender@inria.fr (Matías R. Bender), jean-charles.faugere@inria.fr (Jean-Charles Faugère), ludovic.perret@lip6.fr (Ludovic Perret), elias.tsigaridas@inria.fr (Elias Tsigaridas)

URL: http://www-polsys.lip6.fr/~bender/ (Matías R. Bender), http://www-polsys.lip6.fr/~jcf/ (Jean-Charles Faugère), http://www-polsys.lip6.fr/~perret/ (Ludovic Perret), http://www-polsys.lip6.fr/~elias/ (Elias Tsigaridas)

1 Some authors, e.g., [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF], refer to this number as the symmetric rank of the tensor.

2

K n , we can construct a homogeneous polynomial in n variables of degree D. We can identify the symmetric tensors of rank-1 with the D-th power of linear forms. Hence, to decompose a symmetric tensor of order D is equivalent to write the corresponding polynomial as a sum of D-th powers of linear forms using the minimal numbers of summands. This minimal number is the rank of the polynomial/tensor.

Under this formulation, symmetric tensor decomposition dates back to the origin of modern (linear) algebra as a part of Invariant Theory. In this setting, the decomposition of generic symmetric tensors corresponds to canonical forms [START_REF] Sylvester | An essay on canonical forms, supplement to a sketch of a memoir on elimination, transformation and canonical forms[END_REF][START_REF] Sylvester | On a remarkable discovery in the theory of canonical forms and of hyperdeterminants[END_REF][START_REF] Gundelfinger | Zur theorie der binären formen[END_REF].

Together with the theory of apolarity, this problem was of great importance because the decompositions provide information about the behavior of the polynomials under linear change of variables [START_REF] Kung | The invariant theory of binary forms[END_REF].

Binary Form Decomposition. We study the decomposition of symmetric tensors of order D and dimension 2. In terms of homogeneous polynomials, we consider a binary form

f (x, y) := ∑ D i=0 D i a i x i y D-i , (1)
where a i ∈ K ⊂ C and K is some field of characteristic zero. We want to compute a decomposition f (x, y) = ∑ r j=1 (α j x + β j y) D ,

where α 1 , . . . , α r , β 1 , . . . , β r ∈ K, with K being the algebraic closure of K, and r is minimal. We say that a decomposition unique if, for all the decompositions, the set of points {(α j , β j ) : 1 ≤ j ≤ r} ⊂ P 1 (K) is unique, where P 1 (K) is the projective space of K (Reznick, 2013a).

Previous work. The decomposition of binary forms, Equation (2), has been studied extensively for K = C. More than one century ago [START_REF] Sylvester | On a remarkable discovery in the theory of canonical forms and of hyperdeterminants[END_REF][START_REF] Sylvester | An essay on canonical forms, supplement to a sketch of a memoir on elimination, transformation and canonical forms[END_REF] described the necessary and sufficient conditions for a decomposition to exist, see Section 2.1. He related the decompositions to the kernel of Hankel matrices. For a modern approach of this topic, we refer to [START_REF] Kung | The invariant theory of binary forms[END_REF]; [START_REF] Kung | Canonical forms of binary forms: variations on a theme of Sylvester[END_REF]; Reznick (2013a); [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF]. Sylvester's work was extended to a more general kind of polynomial decompositions that we do not consider in this work, e.g., [START_REF] Gundelfinger | Zur theorie der binären formen[END_REF]; [START_REF] Reznick | Homogeneous polynomial solutions to constant coefficient pde's[END_REF]; [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF].

Sylvester's results lead to an algorithm (Algorithm 1) to decompose binary forms (see Comon and Mourrain, 1996, Sec. 3.4.3). In the case where the binary form is of odd degree, then we can compute the decompositions using Berlekamp-Massey algorithm (see [START_REF] Dür | On computing the canonical form for a binary form of odd degree[END_REF]. When the decomposition is unique, the Catalecticant algorithm, which also works for symmetric tensors of bigger dimension [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF][START_REF] Oeding | Eigenvectors of tensors and algorithms for waring decomposition[END_REF], improves

Sylvester's work. For an arbitrary binary form, [START_REF] Helmke | Waring's problem for binary forms[END_REF] presented a randomized algorithm based on Padé approximants and continued fractions, in which he also characterized the different possible decompositions. Unfortunately, all these algorithms have complexity at least quadratic in the degree of the binary form.

Besides the problem of computing the decomposition(s) many authors considered the subproblems of computing the rank and deciding whether there exists a unique decomposition, e.g., [START_REF] Sylvester | On a remarkable discovery in the theory of canonical forms and of hyperdeterminants[END_REF][START_REF] Sylvester | An essay on canonical forms, supplement to a sketch of a memoir on elimination, transformation and canonical forms[END_REF]; [START_REF] Helmke | Waring's problem for binary forms[END_REF]; [START_REF] Comas | On the rank of a binary form[END_REF]; [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF]. For example, [START_REF] Sylvester | On a remarkable discovery in the theory of canonical forms and of hyperdeterminants[END_REF][START_REF] Sylvester | An essay on canonical forms, supplement to a sketch of a memoir on elimination, transformation and canonical forms[END_REF] considered generic binary forms, that is binary forms with coefficients belonging to a dense algebraic open subset of K D+1 [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF] Section 3), and proved that when the degree is 2k or 2k + 1, for k ∈ N, the rank is k + 1 and that the minimal decomposition is unique only when the degree is odd. In the non-generic case, [START_REF] Helmke | Waring's problem for binary forms[END_REF]; [START_REF] Comas | On the rank of a binary form[END_REF]; [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF], among others, proved that the rank is related to the kernel of a Hankel matrix and that the decomposition of a binary form of degree 2k or 2k -1 and rank r, is unique if and only if r ≤ k. With respect to the problem of computing the rank there are different variants of algorithms, e.g., [START_REF] Comas | On the rank of a binary form[END_REF]; [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF]; [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF]. Even though there are not explicit complexity estimates, by exploiting recent superfast algorithms for Hankel matrices [START_REF] Pan | Structured matrices and polynomials: unified superfast algorithms[END_REF], we can deduce a nearly-optimal arithmetic complexity bound for computing the rank using the approach of [START_REF] Comas | On the rank of a binary form[END_REF].

For the general problem of symmetric tensor decomposition, Sylvester's work was successfully extended to cases in which the decomposition is unique, e.g., [START_REF] Brachat | Symmetric tensor decomposition[END_REF]; [START_REF] Oeding | Eigenvectors of tensors and algorithms for waring decomposition[END_REF]. There are also homotopy techniques to solve the general problem, e.g., to decompose generic symmetric tensors [START_REF] Bernardi | Tensor decomposition and homotopy continuation[END_REF] or, when there is a finite number of possible decompositions and we know at least one of them, to compute all the other decompositions [START_REF] Hauenstein | Homotopy techniques for tensor decomposition and perfect identifiability[END_REF]. There are no complexity estimations for these methods. Besides tensor decomposition, there are other related decompositions for binary forms and univariate polynomials that we do not consider in this work, e.g., [START_REF] Reznick | Homogeneous polynomial solutions to constant coefficient pde's[END_REF]Reznick ( , 2013b)); [START_REF] Giesbrecht | Algorithms for computing sparsest shifts of polynomials in power, chebyshev, and pochhammer bases[END_REF]; [START_REF] Giesbrecht | Interpolation of shifted-lacunary polynomials[END_REF]; [START_REF] García-Marco | Reconstruction algorithms for sums of affine powers[END_REF].

Formulation of the problem. Instead of decomposing the binary form as in Equation ( 2), we compute λ 1 . . . λ r , α 1 . . . α r , β 1 . . . β r ∈ K, where r is minimal, such that,

f (x, y) = ∑ r j=1 λ j (α j x + β j y) D . (3)
Since every λ j belongs to the algebraic closure of the field K, the problems are equivalent.

This approach allows us to control the algebraic degree [START_REF] Bajaj | The algebraic degree of geometric optimization problems[END_REF][START_REF] Nie | The algebraic degree of semidefinite programming[END_REF] of the parameters λ j , α j , and β j in the decomposition (Section 4.1).

Note that if the field is not algebraically closed and we force the parameters to belong to the base field, that is λ j , α j , β j ∈ K, the decompositions induced by Equation (2) and Equation ( 3) are not equivalent. We do not consider the latter case and we refer to [START_REF] Helmke | Waring's problem for binary forms[END_REF]; [START_REF] Reznick | Sums of even powers of real linear forms[END_REF]; [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF]; [START_REF] Boij | Monomials as sums of powers: the real binary case[END_REF]; [START_REF] Blekherman | Typical real ranks of binary forms[END_REF] for K = R, and to [START_REF] Reznick | Homogeneous polynomial solutions to constant coefficient pde's[END_REF]Reznick ( , 2013a)); [START_REF] Reznick | Binary forms with three different relative ranks[END_REF] for K ⊂ C.

Main results. We extend Sylvester's algorithm to achieve a nearly-optimal complexity bound in the degree of the binary form. By considering structural properties of the Hankel matrices, we restrict the possible values for the rank of the decompositions and we identify when the decomposition is unique. We build upon [START_REF] Helmke | Waring's problem for binary forms[END_REF] and we use the Extended Euclidean Algorithm to deduce a better complexity estimate than what was previously known. Similarly to Sylvester's algorithm, our algorithm decomposes successfully any binary form, without making any assumptions on the input.

First, we focus on symbolic decompositions, that is a representation of the decomposition as a sum of a rational function evaluated at the roots of a univariate polynomial (Definition 36).

We introduce an algorithm to compute a symbolic decomposition of a binary form of degree

D in O(M(D) log(D))
, where M(D) is the arithmetic complexity of polynomial multiplication (Theorem 43). When the decomposition is unique, the algorithm is deterministic and this is a worst case bound. When the decomposition is not unique, our algorithm makes some random choices to fulfill certain genericity assumptions; thus the algorithm is a Monte Carlo one. However, we can verify if the genericity assumptions hold within the same complexity bound, that is

O(M(D) log(D))
, and hence we can also deduce a Las Vegas variant of the algorithm.

Following the standard terminology used in structured matrices [START_REF] Pan | Structured matrices and polynomials: unified superfast algorithms[END_REF], our algorithm is superfast as its arithmetic complexity matches the size of the input up to poly-logarithmic factors. The symbolic decomposition allow us to approximate the terms in a decomposition, with a relative error of 2 -ε , in O D log 2 (D) log 2 (D) + log(ε) arithmetic operations [START_REF] Pan | Univariate polynomials: Nearly optimal algorithms for numerical factorization and root-finding[END_REF][START_REF] Mcnamee | Numerical methods for roots of polynomials[END_REF]. Moreover, we can deduce for free the rank and the border rank of the tensor, see [START_REF] Comas | On the rank of a binary form[END_REF], Section 1).

Using results from [START_REF] Kaltofen | Improved sparse multivariate polynomial interpolation algorithms[END_REF], we bound the algebraic degree of the decompositions by min(rank, Drank + 1) (Theorem 28). Moreover, we prove lower bounds for the algebraic degree of the decomposition and we show that in certain cases the bound is tight (Sec- This Boolean worst case bound holds independently of whether the decomposition is unique or not.

This work is an extension of the conference paper [START_REF] Bender | A superfast randomized algorithm to decompose binary forms[END_REF]. With respect to the conference version, our main algorithm (Algorithm 3) omits an initial linear change of coordinates as we now rely on fewer genericity assumptions. In contrast with our previous algorithm, we present an algorithm which is deterministic when the decomposition is unique (Theorem 43).

When the decomposition is not unique, our algorithm is still randomized but we present bounds for the number of bad choices that it could make (Proposition 29). With respect to the algebraic degree of the problem, we study the tightness of the bounds that we proposed in the conference paper (Theorem 27). We introduce explicit lower bounds showing that our bounds can be tight (Section 4.1.3).

Organization of the paper. First, we introduce the notation. In Section 2, we present the preliminaries that we need for introducing our algorithm. We present Sylvester's algorithm (Section 2.1), we recall some properties of the structure of the kernel of the Hankel matrices (Section 2.2), we analyze its relation to rational reconstructions of series/polynomials (Section 2.3), and we present the Extended Euclidean Algorithm (Section 2.4). Later, in Section 3, we present our main algorithm to decompose binary forms (Algorithm 3) and its proof of correctness (Section 3.3).

This algorithm uses Algorithm 4 to compute the kernel of a family of Hankel matrices, which we consider in Section 3.1. Finally, in Section 4, we study the algebraic degree of the problem (Section 4.1), we present tight bounds for it (Section 4.1.

3), and we analyze the arithmetic (Section 4.2) and bit complexity of Algorithm 3 (Section 4.3).

Notation. We denote by O, respectively O B , the arithmetic, respectively bit, complexity and we use O, respectively O B , to ignore (poly-)logarithmic factors. M(n) is the arithmetic complexity of multiplying two polynomial of degree n. Let K be a zero characteristic subfield of C, and K its algebraic closure. If v = (v 0 , . . . , v n ) T then P v = P (v 0 ,...,v n ) := ∑ n i=0 v i x i y n-i . Given a binary form f (x, y), we denote by f (x) the univariate polynomial f (x) := f (x, 1). By f (x) we denote the derivative of f (x) with respect to x. For a matrix M, rk(M) is its rank and Ker(M) its kernel.

Preliminaries

An algorithm based on Sylvester's theorem

Sylvester's theorem (Theorem 2) relates the minimal decomposition of a binary form to the kernel of a Hankel matrix. Moreover, it implies an (incremental) algorithm for computing the minimal decomposition. The version that we present in Algorithm 1 comes from Comon and Mourrain (1996, Section 3.2).

Definition 1. Given a vector a = (a 0 , . . . , a D ) T , we denote by {H k a } 1≤k≤D the family of Hankel matrices indexed by k, where H k a ∈ K (D-k+1)×(k+1) and

H k a :=            a 0 a 1 • • • a k-1 a k a 1 a 2 • • • a k a k+1 . . . . . . . . . . . . . . . a D-k-1 a D-k • • • a D-2 a D-1 a D-k a D-k+1 • • • a D-1 a D            . ( 4
)
We may omit the index a in H k a when it is clear from the context.

Theorem 2 [START_REF] Sylvester | On a remarkable discovery in the theory of canonical forms and of hyperdeterminants[END_REF]. Let f (x, y) = ∑ D i=0 D i a i x i y D-i with a i ∈ K ⊆ C. Also, consider a non-zero vector c = (c 0 , . . . , c r ) T ∈ K r+1 , such that the polynomial

P c = ∑ r i=0 c i x i y r-i = ∏ r j=1 (β j x -α j y)
is square-free and α j , β j ∈ K, for all 1 ≤ j ≤ r. Then, there are λ 1 , . . . λ r ∈ K such that we can write f (x, y) as

f (x, y) = r ∑ j=1 λ j (α j x + β j y) D ,
if and only if (c 0 , . . . , c r ) T ∈ Ker(H r a ).

Algorithm 1 INCRDECOMP [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF], Figure 1)

1. r := 1 2. Get a random c ∈ Ker(H r )
3. If P c is not square-free, r := r + 1 and GO TO 2 4. Write P c as ∏ r j=1 (β j xα j y)

5. Solve the transposed Vandermonde system:

        β D 1 • • • β D r β D-1 1 α 1 • • • β D-1 r α r . . . . . . . . . α D 1 • • • α D r              λ 1 . . . λ r      =      a 0 . . . a D      . (5) 6. Return ∑ r j=1 λ j (α j x + β j y) D
For a proof of Theorem 2 we refer to Reznick (2013a, Theorem 2.1 & Corollary 2.2). Theorem 2 implies Algorithm 1. This algorithm will execute steps 2 and 3 as many times as the rank.

At the i-th iteration it computes the kernel of H i . The dimension of this kernel is ≤ i and each vector in the kernel has i + 1 coordinates. As the rank of the binary form can be as big as the degree of the binary form, a straightforward bound for the arithmetic complexity of Algorithm 1 is at least cubic in the degree.

We can improve the complexity of Algorithm 1 by a factor of D by noticing that the rank of the binary form is either rk(H [START_REF] Comas | On the rank of a binary form[END_REF], Section 3) (Helmke, 1992, Theorem B). Another way to compute the rank is by using minors (Bernardi et al., 2011, Algorithm 2).

D 2 ) or D -rk(H D 2 ) + 2 (
The bottleneck of the previous approaches is that they have to compute the kernel of a Hankel matrix. However, even if we know that the rank of the binary form is r, the dimension of the kernel of H r can still be as big as O(D); the same bound holds for the length of the vectors in the kernel. Hence, the complexity is lower bounded by O(D 2 ).

Our approach avoids the incremental construction. We exploit the structure of the kernel of the Hankel matrices and we prove that the rank has only two possible values (Lemma 17), see also (Comas and Seiguer, 2011, Section 3), (Helmke, 1992, Theorem B), or [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF]. Moreover, we use a compact representation of the vectors in the kernel. We describe them as a combination of two polynomials of degree O(D).

Kernel of the Hankel matrices

The Hankel matrices are among the most studied structured matrices [START_REF] Pan | Structured matrices and polynomials: unified superfast algorithms[END_REF]. They are related to polynomial multiplication. We present results about the structure of their kernel. For details, we refer to Heinig and Rost (1984, Chapter 5).

Proposition 3. Matrix-vector multiplication of Hankel matrices is equivalent to polynomial multiplication. Given two binary forms A := ∑ D i=0 a i x i y D-i and U := ∑ k i=0 u i x i y k-i , consider R := ∑ D+k i=0 r i x i y D+k-i = A •U. If we choose the monomial basis {y D+k , . . . , x D+k }, then the equality A • U = R is equivalent to Equation (6), where the central submatrix of the left matrix is H k (a 0 ,...,a D ) (Definition 1).

                                 a 0 a 0 a 1 . . . . . . . . . a 0 • • • a k-2 a k-1 a 0 a 1 • • • a k-1 a k a 1 a 2 • • • a k a k+1 . . . . . . . . . . . . . . . a D-k a D-k+1 • • • a D-1 a D a D-k+1 a D-k+2 • • • a D . . . . . .      r 0 r 1 . . . r k-1 r k r k+1 . . . r D r D+1 . . .
Consider a family of Hankel matrices {H k a } 1≤k≤D as in Definition 1. There is a formula for the dimension of the kernel of each matrix in the family {H k a } 1≤k≤D that involves two numbers, N a 1 and N a 2 . To be more specific, the following holds:

Proposition 4. For any family of Hankel matrices {H k a } 1≤k≤D there are two constants, N a 1 and N a 2 , such that the following hold:

1. 0 ≤ N a 1 ≤ N a 2 ≤ D. 2. For all k, 1 ≤ k ≤ D, it holds dim(Ker(H k a )) = max(0; k -N a 1 ) + max(0; k -N a 2 ). 3. N a 1 + N a 2 = D.
We may omit the index a in N a 1 and N a 2 when it is clear from the context.

rk(H k ) The graphs of the functions consist of three line segments, as we can see in the Figures 1a and1b. The dimension of the kernel is an increasing function of k. For k from 1 to N 1 , the kernel of the matrix is trivial, so the rank increases with the number of columns. That is, the slope of the graph of the rank (Figure 1a) is 1, while the slope of the graph of the dimension of the kernel (Figure 1b) is 0. For k from N 1 + 1 to N 2 , the rank remains constant as for each column that we add, the dimension of the kernel increases by one. Hence, the slope of the graph of the rank is 0 and the slope of the graph of the dimension of the kernel is 1. For k from N 2 + 1 to D, the rank decreases because the dimension of the kernel increases by 2, and so the slope of the graph of the rank is -1, while the slope of the graph of the dimension of the kernel is 2.

N 1 + 1 1 N 1 N 2 D 2 D k 1 (a) Rank of the Hankel matrices dim(Ker(H k )) N 1 N 2 D D N 2 -N 1 (b) Dimension of the kernel
If N 1 = N 2 , then the graph of both functions degenerates to two line segments. Regarding the graph of the rank, the first segment has slope 1 for k from 1 to N 1 + 1 and the second segment has slope -1 for k from N 1 + 1 to D. For the graph of the dimension of the kernel, the first segment has slope 0 from 1 to N 1 + 1, and the second one has slope 2 from N 1 to D.

The elements of the kernel of the matrices in {H k } are related. To express this relation from a linear algebra point of view, we need to introduce the U-chains.

Definition 5 [START_REF] Heinig | Algebraic methods for Toeplitz-like matrices and operators[END_REF]

, Definition 5.1)). A U-chain of length k of a vector v = (v 0 , . . . , v n ) T ∈ K n+1 is a set of vectors {U 0 k v, U 1 k v, . . . , U k-1 k v} ⊂ K n+k . The i-th element, for 0 ≤ i ≤ k -1, is U i k v = (0, . . . , 0 i , n+1 v 0 , . . . , v n , 0, . . . , 0 k-1-i ),
where (Heinig and Rost, 1984, page 11).

U i k is an i-shifting matrix of dimension (n + k) × (n + 1)
If v is not zero, then all the elements in a U-chain of v are linearly independent. The following theorem uses U-chains to relate the vectors in the kernels of a family of Hankel matrices.

Proposition 6 (Vectors v and w). Given a family of Hankel matrices {H k } 1≤k≤D , let N 1 and N 2 be the constants of Proposition 4. There are two vectors, v ∈ K N 1 +2 and w ∈ K N 2 +2 , such that

• If 0 ≤ k ≤ N 1 , then Ker(H k ) = {0}. • If N 1 < k ≤ N 2 , then the U-chain of v of length (k -N 1 ) is a basis of Ker(H k ), that is Ker(H k ) = U 0 k-N 1 v, . . . , U k-N 1 -1 k-N 1 v . • If N 2 < k ≤ D, then the U-chain of v of length k -N 1 together with the U-chain of w of length k -N 2 is a basis of Ker(H k ), that is Ker(H k ) = U 0 k-N 1 v, . . . , U k-N 1 -1 k-N 1 v, U 0 k-N 2 w, . . . , U k-N 2 -1 k-N 2 w .
The vectors v and w of Proposition 6 are not unique. The vector v could be any vector in Ker(H N 1 +1 ). The vector w could be any vector in Ker(H N 2 +1 ) that does not belong to the vector space generated by the U-chain of v of length N 2 -N 1 + 1. From now on, given a family of

Hankel matrices, we refer to v and w as the vectors of Proposition 6.

Let u be a vector in the kernel of H k and P u its corresponding polynomial (see Notation).

We say that P u is a kernel polynomial. As P U j k v = x j y k-1-j P v , we can write any kernel polynomial of a family of Hankel matrices as a combination of P v and P w (Heinig and Rost, 1984, Propositions 5.1 & 5.5). Moreover, P v and P w are relatively prime.

Proposition 7. Consider any family of Hankel matrices {H k } 1≤k≤D . The kernel polynomials P v and P w are relative prime. Moreover, for each k, the set of kernel polynomials of the matrix H k is as follows:

• If 0 < k ≤ N 1 , then it is {0}. • If N 1 < k ≤ N 2 , then it is {P µ • P v : µ ∈ K k-N 1 }. • If N 2 < k ≤ D, then it is {P µ • P v + P ρ • P w : µ ∈ K k-N 1 , ρ ∈ K k-N 2 }. Corollary 8. Let ω ∈ Ker(H N 2 +1 ) such that P ω {P µ • P v : µ ∈ K N 2 -N 1 +1
}, then we can consider ω as the vector w from Proposition 6.

Rational Reconstructions

A rational reconstruction for a series, respectively a polynomial, consists in approximating the series, respectively the polynomial, as the quotient of two polynomials. Rational reconstructions are the backbone of many problems e.g., Padé approximants, Cauchy Approximations, Linear Recurrent Sequences, Hermite Interpolation. They are related to Hankel matrices. For an introduction to rational reconstructions, we refer to Bostan et al. (2017, Chapter 7) and references therein.

Definition 9. Consider a := (a 0 , . . . , a D ) T ∈ K D+1 and a polynomial A :

= ∑ D i=0 a i x i ∈ K[x]
. Given a pair of univariate polynomials (U, R), we say that they are a rational reconstruction of A

modulo x D+1 if A •U ≡ R mod x D+1 .
Such a reconstruction is not necessarily unique. Our interest emanates from the relation between the rational reconstructions of A modulo x D+1 and the kernels of the family of Hankel matrices {H k a } k .

Lemma 10. Consider ω ∈ Ker(H k a ) and (r 0 , . . . ,

r k-1 ) ∈ K k such that      0 • • • 0 a 0 . . . . . . . . . . . . 0 a 0 • • • a k-1              ω 0 ω 1 . . . ω k         =      r 0 . . . r k-1      . Then, (P ω (1, x), ∑ k-1 i=0 r i x i ) is a rational reconstruction of A modulo x D+1 . Proof. Following Equation (6), if ω ∈ Ker(H k a ), then                a 0 . . . . . . a 0 • • • a k-1 a 0 a 1 • • • a k . . . . . . . . . . . . a D-k a D-k+1 • • • a D                        ω 0 ω 1 . . . ω k         =                r 0 . . . r k-1 0 . . . 0                . (7) Hence, P ω (1, x) = ∑ k i=0 ω k-i x i and A • P ω (1, x) ≡ ∑ k-1 i=0 r i x i mod x D+1 .
Therefore,

(P ω (1, x), ∑ k-1 i=0 r i x i ) is a rational reconstruction of A modulo x D+1 .
Lemma 11. If (U, R) is a rational reconstruction of A module x D+1 , then there is a vector

ω ∈ Ker(H max(deg(U),deg(R)+1) a
) such that

P ω = U y x x max(deg(U),deg(R)+1) . Proof. Let k = deg(U), q = deg(R), U = ∑ i u i x i and R = ∑ i r i x i . Following Equation (6), AU ≡ R mod x D+1 is equivalent to,               a 0 . . . . . . a 0 • • • a k-1 a 0 a 1 • • • a k . . . . . . . . . . . . a D-k a D-k+1 • • • a D                       u k . . . u 1 u 0         =                  r 0 r 1 . . . r q 0 . . . 0                  . ( 8
)
If k > q, Equation ( 8) reduces to Equation ( 7), and so ω = (u k , . . . , u 0 ) ∈ Ker(H k a ). Moreover,

U y x x k = k ∑ i=0 u i y i x k-i = ( j↔k-i) k ∑ j=0 u k-j x j y j-k = P ω .
If q ≥ k, we extend the vector (u k , . . . , u 0 ) by adding (q + 1 -k) leading zeros. We rewrite Equation ( 8) as Equation ( 9). The concatenation of the two bottom submatrices form the matrix H q+1 a , and so ω = (0, . . . , 0, u k , . . . , u 0 ) ∈ Ker(H q+1 a ). Also,

P ω = k ∑ j=0 u j x q+1-j y j + q+1 ∑ j=k+1 0 x q+1-j y j = U y x x q+1 .                       a 0 . . . . . . a 0 • • • a k a 0 a 1 • • • a k+1 . . . . . . . . . . . . . . . a 0 . . . a q-k a q+1-k . . . a q+1 . . . . . . . . . . . . . . . . . . a D-q-1 . . . a D-k-1 a D-k • • • a D                                      0 . . . 0 u k . . . u 0                =                      r 0 . . . r k . . . r q 0 . . . 0                      . ( 9
)
Remark 12. If (U, R) is a rational reconstruction, then the degree of the kernel polynomial

P ω (x, y) = U y x x max(deg(U),deg(R)+1) is max(deg(U), deg(R) + 1).
In particular, the maximum power of x that divides the kernel polynomial P ω is x max(0,deg(R)+1-deg(U)) .

Greatest Common Divisor and Bézout identity

The Extended Euclidean algorithm (EGCD) is a variant of the classical Euclidean algorithm that computes the Greatest Common Divisor of two univariate polynomials A and B, gcd(A, B), together with two polynomials U and V , called cofactors, such that U A +V B = gcd(A, B). In the process of computing these cofactors, the algorithm computes a sequence of relations between A and B that are useful to solve various problems, in particular to compute the rational reconstruction of A modulo B. For a detailed exposition of this algorithm, we refer to Bostan et al. (2017, Chapter 6) and Gathen and Gerhard (2013, Chapter 3 and 11).

Algorithm 2 Calculate the EGCD of A and B

(U 0 ,V 0 , R 0 ) ← (0, 1, B) (U 1 ,V 1 , R 1 ) ← (1, 0, A) k ← 1 while R k 0 do k ← k + 1 Q k-1 ← R k-2 quo R k-1 (U k ,V k , R k ) ← (U k-2 ,V k-2 , R k-2 ) -Q k-1 (U k-1 ,V k-1 , R k-1 ) end while Return {(U i ,V i , R i )} i The Extended Euclidean Algorithm (Algorithm 2) computes a sequence of triples {(U i ,V i , R i )} i which form the identities U i A +V i B = R i , for all i. (10)
Following [START_REF] Gathen | Modern computer algebra[END_REF], we refer to these triplets as the rows of the Extended Euclidean algorithms of A and B. Besides Equation ( 10), the rows are related to each other as follows.

Remark 13. The degrees of the polynomials {R i } i form a strictly decreasing sequence, that is

deg(R i ) > deg(R i+1
) for every i.

Lemma 14 (Bostan et al., 2017, Sec 7.1). For each i, U i V i+1 -U i+1 V i = (-1) i , and so the polynomials U i and V i are coprime.

Lemma 15 (Bostan et al., 2017, Lem 7.1). For each i > 0, the degree of U i is the degree of B minus the degree of R i-1 , that is

deg(U i ) = deg(B) -deg(R i-1 ), ∀ i > 0.
Every row of the Extended Euclidean Algorithm leads to rational reconstruction of A modulo B.

Remark 16. For each i ≥ 0, (U i , R i ) is a rational reconstruction of A modulo B.

The Algorithm

One of the drawbacks of Algorithm 1 and its variants is that they rely on the computation of the kernels of many Hankel matrices and they ignore the particular structure that it is present in all of them. Using Lemma 17, we can skip many calculations by computing only two vectors, v and w (Proposition 6). This is the main idea behind Algorithm 3 that leads to a softly-linear arithmetic complexity bound (Section 4.2).

Algorithm 3 performs as follows: First, step 1 computes the kernel polynomials P v and P w which, by Proposition 7, allow us to obtain the kernel polynomials of all the Hankel matrices (see Section 3.1). Then, step 2 computes a square-free kernel polynomial of the minimum degree r (see Section 3.2). Next, step 3 computes the coefficients λ 1 , . . . , λ r (see Section 4.1.2). Finally, step 4 recovers a decomposition for the original binary form.

Let f be a binary form as in Equation ( 1) and let {H k } 1≤k≤D be its corresponding family of Hankel matrices (see Definition 1). The next well-known lemma establishes the rank of f .

Lemma 17. Assume f , {H k } k , N 1 and N 2 of Proposition 4, and v and w of Proposition 6.

If P v (Proposition 7) is square-free then the rank of f is N 1 + 1, else, it is N 2 + 1.
Proof. By Proposition 4, for k < N 1 + 1, the kernel of H k is trivial. Hence, by Sylvester's theorem (Theorem 2), there is no decomposition with a rank smaller than N 1 + 1. Recall that v ∈ Ker(H N 1 +1 ). So, if P v is square-free, by Sylvester's theorem, there is a decomposition of rank N 1 + 1.

Assume P v is not square-free. For N 1 + 1 ≤ k ≤ N 2 , P v divides all the kernel polynomials of the matrices H k (Proposition 7). Therefore, none of them is square-free, and so the rank is at least N 2 + 1.

By Proposition 7, P v and P w are coprime. So, there is a polynomial P µ of degree N 2 -N 1 such that Q µ := P v • P µ + P w is square-free. A formal proof of this appears in Theorem 22.

By Proposition 7, Q µ is a square-free kernel polynomial of degree N 2 + 1. Consequently, by

Sylvester's theorem, there is a decomposition with rank N 2 + 1.

For alternative proofs of Lemma 17 we refer to (Helmke, 1992, Theorem B), (Comas and Seiguer, 2011, Section 3), [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF], or (Carlini et al., 2018, Section 4).

To relate Lemma 17 with the theory of binary form decomposition, we recall that the decompositions are identified with the square-free polynomials in the annihilator of the ideal f [START_REF] Kung | The invariant theory of binary forms[END_REF]; [START_REF] Iarrobino | Power sums, Gorenstein algebras, and determinantal loci[END_REF], Chapter 1). All the kernel polynomials of {H k } k belong to the annihilator of f and they form an ideal. If f is a binary form of degree D = 2k or 2k + 1, then this ideal is generated by two binary forms of degrees rk(H k ) and D + 2rk(H k ),

with no common zeros (Iarrobino and Kanev, 1999, Theorem 1.44). These are the polynomials P v and P w . Using this interpretation Algorithm 1 and its variants compute a (redundant) generating set of the annihilator of f , while Algorithm 3 computes a (minimal) basis.

Algorithm 3 FASTDECOMP

Input: A binary form f (x, y) of degree D

Output: A decomposition for f (x, y) of rank r.

1. Compute P v and P w of {H k a } k

We use Algorithm 4 with (a 0 , . . . , a D ).

IF P

v (x, y) is square-free, Q ←-P v r ←-N 1 + 1 {The rank of the decomposition is the degree of Q}

ELSE

Compute a square-free binary form Q

We compute a vector µ of length (N 2 -N 1 + 1), such that (P µ • P v + P w ) is square-free (Section 4.1.1).

Q ←-P µ • P v + P w r ←-N 2 + 1 {The rank of the decomposition is the degree of Q} 3. Compute the coefficients λ 1 , . . . , λ r

Solve the system of Equation ( 5) where Q(x, y) = ∏ r j=1 (β j xα j y).

For details and the representation of λ j , see Section 4.1.2.

4. Return f (x, y) = ∑ r j=1 λ j (α j x + β j y) D

3.1. Computing the polynomials P v and P w

We use Lemma 10 and Lemma 11 to compute the polynomials P v and P w from Proposition 7 as a rational reconstruction of A := ∑ D i=0 a i x i modulo x D+1 . Our algorithm exploits the Extended Euclidean Algorithm in a similar way as [START_REF] Cabay | Algebraic computations of scaled Padé fractions[END_REF] do to compute scaled Padé fractions.

In the following, let v be the vector of Proposition 6, consider

U v := P v (1, x) and R v ∈ K[x]
as the remainder of the division of (A • P v (1, x)) by x D+1 . Note that the polynomial R v is the unique polynomial of degree strictly smaller to

N 1 + 1 such that A • P v (1, x) ≡ R v mod x D+1 , see Equation (7). Lemma 18. If (U, R) is a rational reconstruction of A modulo x D+1 such that max(deg(U), deg(R) + 1) ≤ N 2 , then there is a polynomial Q ∈ K[x] such that Q • P v (x, 1) = U and Q • R v = R.
Proof. Let k := deg(U) and q := deg(R). By Lemma 11, there is a non-zero vector

ω ∈ Ker(H max(k,q+1) a
) such that the kernel polynomial P ω is equal to U y x x max(k,q+1) . Hence,

Ker(H max(k,q+1) a
) 0 and so, by Proposition 6, N 1 < max(k, q + 1). We assume that max(k, q + 1) ≤ N 2 , hence the degree of P ω is smaller or equal to N 2 and, by Proposition 7, P ω is divisible by P v . Therefore, there is a polynomial Q

∈ K[x, y] such that QP v = P ω . Let Q := Q(1, x). By definition, U v = P v (1, x) and U = P ω (1, x), so U = QU v . Hence, Q R v ≡ R mod x D+1 , because R v ≡ U v A mod x D+1 and Q R v ≡ QU v A ≡ UA ≡ R mod x D+1 . If the degrees of (Q R v ) and R are smaller than D + 1, then Q R v = R,
deg(Q R v ) ≤ deg(U v Q R v ) ≤ N 2 + N 1 = D (Proposition 4).
We can use this lemma to recover the polynomial P v from certain rational reconstructions.

Corollary 19. If (U, R) is a rational reconstructions of A modulo x D+1 such that max(deg(U), deg(R) + 1) ≤ N 2 and for every polynomial Q of degree strictly bigger than zero that divides U and R,

( U Q , R Q ) is not a rational reconstruction of A modulo x D+1 , then there is a non-zero constant c such that P v = c •U y x x max(deg(U),deg(R)+1) (Proposition 7). In particular, N 1 = max(deg(U) -1, deg(R)).
Proof. By Lemma 18, there is a Q

∈ K[x] such that Q • (U v , R v ) = (U, R). By Lemma 10, (U v , R v
) is a rational reconstruction, and so deg (Bostan et al., 2017, Section 7.1). When this Padé approximant exists, it is unique, meaning that for any rational reconstruction with this property the quotient R U is unique (we can invert U mod x D+1 because U(0) = 1). When N 1 < N 2 , we have that D+1 2 ≤ N 2 (Proposition 4) and so, if the the Padé approximant of A of type ( D+1 2 -1, D+1 2 ) exists, by Lemma 18, we can recover P v from it. The existence of this Padé approximant is equivalent to the condition U v (0) = 1, which means v N 1 +1 = 1. In the algorithm proposed in the conference version of this paper (Bender et al., 2016, Algorithm 3), the correctness of our algorithms relied on this condition. In that version, we ensured this property with a generic linear change of coordinates in the original polynomial f . In this paper, we skip this assumption. Following Bostan et al. (2017, Theorem 7.2), when N 1 < N 2 , we can compute v no matter the value of v N 1 +1 . This approach has a softly-linear arithmetic complexity and involves the computation of a row of the EGCD of A and x D+1 . We can compute P w from a consecutive row.

(Q) = 0. Hence, N 1 +1 = deg(P v ) = max(deg(U), deg(R) + 1) and Q • P v (1, y x )x N 1 +1 = U( y x )x N 1 +1 . If (U, R) is a rational reconstruction of A modulo x D+1 such that deg(U) + deg(R) ≤ D and U(0) = 1, then R U is the Padé approximant of A of type (deg(R), deg(U))
Before going into the proof, we study the case N 1 = N 2 . In this case, there are not rational reconstructions with the prerequisites of Lemma 18, and so we treat this case in a different way. Proof. First note that, as D = N 1 + N 2 (Proposition 4), then N 1 = D 2 . Following Equation ( 6), if we write U = ∑ N 1 i=0 u i x i and R = ∑ N 1 i=0 u i x i , then we get the linear system,

H N 1                                  a 0 a 0 a 1 . . . . . . . . . a 0 • • • a N 1 -1 a N 1 a 1 • • • a N 1 a N 1 +1 . . . . . . . . . . . . a D-N 1 • • • a D-1 a D                           u N 1 . . . u 1 u 0         =                   r 0 . . . r N 1 -1 r N 1 0 . . . 0                  
.

The matrix H N 1 ∈ K (D-D rational reconstruction. Hence, we can consider deg(U) = N 1 . If R is monic, then r N 1 = 1 and so we compute the coefficients of U and R as

        u N 1 . . . u 1 u 0         = (H N 1 +1 ) -1         1 0 . . . 0         ,         r 0 . . . r N 1 -1 1         =         a 0 a 0 a 1 . . . . . . . . . a 0 • • • a N 1 -1 a N 1                 u N 1 u N 1 -1 . . . u 0         .
Lemma 21 (Correctness of Algorithm 4). Let {(U j ,V j , R j )} j be the set of triplets obtained from the Extended Euclidean Algorithm for the polynomials A and x D+1 , see Section 2.4. Let i be the index of the first row of the extended Euclidean algorithm such that deg(R i ) < D+1 2 . Then, we can compute the polynomials P v and P w of Proposition 7 as

(A) P v = U i ( x y ) • x max(deg(U i ),deg(R i )+1) . (B) If deg(R i ) > deg(U i ), P w = U i+1 ( x y ) • x deg(U i+1 ) . (C) If deg(R i ) ≤ deg(U i ), P w = U i-1 ( x y ) • x deg(R i-1 +1) .
Proof. (A). First observe that if i is the first index such that the degree of R i is strictly smaller than D+1 2 , then, by Remark 13, the degree of R i-1 has to be bigger or equal to D+1 2 . Hence, the degree of U i is smaller or equal to D+1 2 , because by Lemma 15, deg

(U i ) = D + 1 -deg(R i-1 ) ≤ D + 1 -D+1 2 = D+1 2 .
We can consider R i-1 , that is i is strictly bigger than 0, because the degree of R 0 = x D+1 is strictly bigger than D+1 2 . If N 1 = N 2 , then D is even and 10). Coupling the two equalities together, we conclude that

N 1 = D 2 (Proposition 4). As D+1 2 = D 2 , deg(R i ) ≤ D 2 and deg(U i ) ≤ D 2 . By Lemma 20, max(deg(U i ), deg(R i ) + 1) = N 1 + 1 and we can consider P v as U i ( y x )x N 1 +1 . If N 1 < N 2 , assume that there is a non-zero Q ∈ K[x] such that Q divides U i and R i and ( U i Q , R i Q ) is a rational reconstruction of A modulo x D+1 . Then, U i Q A ≡ R i Q mod x D+1 and so there is a poly- nomial V such that V x D+1 + U i Q A = R i Q . Multiplying by Q, we get the equality Q V x D+1 +U i A = R i . Consider the identity V i x D+1 + U i A = R i from Equation (
V i = Q V . As Q divides U i and V i , which are coprime (Lemma 14), Q is a constant, that is deg(Q) = 0. If N 1 < N 2 , then D < 2N 2 (Proposition 4). Hence, deg(U i ) ≤ D+1 2 ≤ N 2 and deg(R i ) + 1 < D+1 2 + 1 ≤ N 2 + 1, that is, max(deg(U i ), deg(R i ) + 1) ≤ N 2 .

Hence, by 20

Lemma 18, we can consider U i ( y x )x max(deg(U i ),deg(R i )+1) as the kernel polynomial P v from Proposition 7, spanning Ker(H N 1 +1 ).

(B). Assume that the degree of U i is strictly bigger than the one of R i , that is deg

(U i ) > deg(R i ). Then N 1 = deg(U i ) -1, as deg(U i ) = deg(P v ) = N 1 + 1 (Remark 12). Note that in this case i > 1 because U 1 = 1, R 1 = A 0, and so deg(U 1 ) ≤ deg(R 1 ). The degree of R i-1 is N 2 because, by Lemma 15, deg(R i-1 ) = D + 1 -deg(U i ) = D + 1 -N 1 -1 = N 2 (Proposition 4). Consider the degree of U i-1 . By Lemma 15, deg(U i-1 ) = D + 1 -deg(R i-2 ). As deg(R i-2 ) > deg(R i-1 ) (Remark 13), then deg(R i-2 ) > N 2 . Therefore, the degree of U i-1 is smaller or equal to the one of R i-1 because deg(U i-1 ) = D + 1 -deg(R i-2 ) < D + 1 -N 2 = N 1 + 1, and so deg(U i-1 ) ≤ N 1 ≤ N 2 = deg(R i-1 ).
Hence, by Remark 16,

(U i-1 , R i-1 ) is a rational reconstruction of A modulo x D+1 such that deg(U i-1 ) ≤ N 1 and deg(R i-1 ) = N 2 . So, max(deg(U i-1 ), deg(R i-1 ) + 1) = N 2 + 1 and, by Re- mark 12, there is a kernel polynomial P ω = U i-1 ( y x )x N 2 +1 of degree N 2 +1 such that x N 2 +1-deg(U i-1 ) divides P ω . As deg(U i-1 ) ≤ N 1 , x N 2 +1-N 1 divides x N 2 +1-deg(U i-1
) and so, it divides P ω . We assumed that the degree of U i is strictly bigger than the one of R i , and so x does not divide P v (Remark 12). Hence, there is no binary form Q of degree N 2 -N 1 such that x N 2 -N 1 +1 divides Q P v . Therefore, by Corollary 8, we can consider P w = P ω .

(C). Assume that the degree of R i is bigger or equal to the one of U i , that is deg(R

i ) ≥ deg(U i ). Hence, deg(R i ) + 1 = deg(P v ) = N 1 + 1 (Remark 12
), and so deg(R i ) = N 1 . In particular, R i 0, and so the (i + 1)-th row of the Extended Euclidean Algorithm,

(U i+1 ,V i+1 , R i+1 ), is defined. The degree of U i+1 is N 2 + 1 because, by Lemma 15, deg(U i+1 ) = D + 1 -deg(R i ) = N 2 + 1 (Proposition 4
). The degree of R i+1 is strictly smaller than the one of R i (Remark 13), which is

N 1 . Hence, the degree of R i+1 is smaller than the degree of U i+1 because deg(R i+1 ) < N 1 ≤ N 2 < deg(U i+1 ). Therefore, P ω = U i+1 ( y x )x N 2 +1 is a kernel polynomial in Ker(H N 2 +1 ) (Lemma 11). By Remark 12, as deg(R i+1 ) < deg(U i+1 ),
x does not divide P ω . Also, the maximal power of x that divides P v is x deg(R i )+1-deg(U i ) , and, as we assumed deg(R

i ) ≥ deg(U i ), x divides P v . Hence, every polynomial in {Q P v : deg(Q) = N 2 -N 1 } is divisible by x,
and so, by Corollary 8, we can consider P w = P ω . tive. Hence, (µ 0 , . . . , µ N 2 -N 1 ) T ∈ {µ ∈ K N 2 -N 1 +1 : Q µ has double roots} if and only if, there is

(µ 0 , . . . , µ N 2 -N 1 , α) ∈ K N 2 -N 1 +1 × K such that it is a solutions to the following system of equa- tions      (P µ • P v + P w )(x) = 0 (P µ • P v + P µ • P v + P w )(x) = 0. (11)
In Equation ( 11), µ 0 only appears in P µ (x, 1) with degree 1. We eliminate it to obtain the polynomial

(P v • P µ + P w )P v -P v • P w .
This polynomial is not identically 0 as P v does not divide P v and P v and P w are relative prime. Hence, for each (µ 1 , . . . , µ N 2 -N 1 ), there is a finite number of solutions for this equation, bounded by the degree of the polynomial. Moreover, as the polynomials of Equation ( 11) are linear in µ 0 , each solution of the deduced equation is extensible to a finite number of solutions of Equation ( 11). Hence, there is a

µ ∈ K N 2 -N 1 +1 , such that Q µ is square-free. Therefore, the discriminant of Q µ (x) is not identically zero.
Corollary 23. For every vector (µ 1 , . . . , µ N 2 -N 1 ) ∈ K N 2 -N 1 such that there is a µ 0 ∈ K such that y 2 does not divides Q µ , where µ = (µ 0 , . . . , µ N 2 -N 1 ), there are at most 2D + 2 different values for µ 0 ∈ K such that the polynomial Q µ (x, y) is not square-free.

Proof. If Q µ (x, y) is not square-free, then it has a double root in P 1 (K). This root could be of the form (α, 1) ∈ P 1 (K) or (1, 0) ∈ P 1 (K). We analyze separately these cases.

First, we consider the polynomial Q µ (x, 1) ∈ K[µ 0 , x]. By Theorem 22, the discriminant of Q µ (x, 1) with respect to x is not zero. As Q µ (x, 1) is a polynomial of degree N 2 + 1 with respect to x, and of degree 1 with respect to µ 0 , the degree with respect to µ 0 of the discriminant of Q µ (x, 1) with respect to x is at most (N 2 + 1) + N 2 ≤ 2D + 1. Hence, there are at most 2D + 1 values for µ 0 such that Q µ (x, y) has a root of the form (α, 1) ∈ P 1 (K) with multiplicity bigger than one.

The polynomial Q µ (x, y) has a root of the form (1, 0) ∈ P 1 (K) with multiplicity bigger than one, if and only if y 2 divides Q µ (x, y). If this happens, then the coefficients of the monomials y • x N 2 -N 1 -1 and x N 2 -N 1 in the polynomial Q µ (x, y) are zero. By assumption, these coefficients are not identically zero as polynomials in K[µ 0 ]. As Q µ (x, y) is a linear polynomial with respect to µ 0 , there is at most one value for µ 0 such that y 2 divides Q µ (x, y).

Therefore, there are at most (2D + 1) + 1 values such that Q µ (x, y) is not square-free.

Remark 24. The previous assumption is not restrictive. If y 2 divides Q µ , where µ = (µ 0 , . . . , µ N 2 -N 1 ), then y 2 does not divide

Q (µ 0 ,...,µ N 2 -N 1 +1) = Q µ + x N 2 +1 nor Q (µ 0 ,...,µ N 2 -N 1 -1 +1,µ N 2 -N 1 ) = Q µ + yx N 2 . Moreover, if N 2 -N 1 ≥ 2, y 2 divides (or not) Q µ (x, y)
regardless the value of µ 0 . Conversely, if N 2 -N 1 < 2, there is always a µ 0 such that y 2 does not divide Q µ .

Correctness of Algorithm 3

For computing a decomposition for a binary form f , we need to compute the kernel of a Hankel matrix (Theorem 2). Algorithm 4 computes correctly the polynomials P v and P w that characterize the kernels of the family of Hankel matrices associated to f . Once we obtain these polynomials, step 2 (see Corollary 23) and step 3 computes the coefficients α j , β j , λ j of the decomposition. Hence, we have a decomposition for f , as

f (x, y) = ∑ r j=1 λ j • (αx + β y) D .
Example. Consider f (x, y) = y 4 + 8xy 3 + 18x 2 y 2 + 16x 3 y + 5x 4 . The family of Hankel matrices associated to f are related to the vector a := (1, 2, 3, 4, 5) T , it is denoted by {H k a } k , and it contains the following matrices:

H 1 a =         1 2 2 3 3 4 4 5         H 2 a =      1 2 3 2 3 4 3 4 5      H 3 a =   1 2 3 4 2 3 4 5   H 4 a = 1 2 3 4 5
The kernel H 1 a is trivial, so we compute the one of H 2 a . This kernel is generated by the vector (1, -2, 1) T , so by Proposition 6 we consider v = (1, -2, 1) T . Also, by Proposition 4,

N 1 + 1 = 2 and N 2 = D-N 1 = 3. The kernel polynomial P v = y 2 -2x y+x 2 = (x-y) 2 is not square-free thus,
by Lemma 17, the rank of f (x, y) is N 2 + 1 = 4 and we have to compute the kernel polynomial P w in the kernel of H 4 a . Following Proposition 6, the kernel of

H 4 a is generated by U-chain of v given vectors U 0 2 v = (1, -2, 1, 0, 0) T , U 1 2 v = (0, 1, -2, 1, 0) T ,
and U 2 2 v = (0, 0, 1, -2, 1) T , plus a vector w linear independent with this U-chain. We consider the vector w = (0, 0, 0, 5, -4), which fulfills that assumption. Hence, P v = y 2 -2x y + x 2 and P w = 5yx 3 -4x 4 .

We proceed by computing a square-free polynomial combination of P v and P w . For that, we choose

Q := (44y 2 + 11yx + 149x 2 ) P v + 36P w = (5x -11y)(x -2y)(x + 2y)(x + y).
Finally, we solve the system given by the transposed of a Vandermonde matrix,

        5 4 1 1 1 11 • 5 3 2 -2 -1 11 2 • 5 2 2 2 (-2) 2 (-1) 2 11 3 • 5 2 3 (-2) 3 (-1) 3                 λ 1 λ 2 λ 3 λ 4         =         1 2 3 4         . ( 12
)
The unique solution of the system is (-1 336 , 3, 1 21 , 3 16 ) T , and so we recover the decomposition

f (x, y) = - 1 336 (11x + 5y) 4 + 3(2x + y) 4 + 1 21 (-2x + y) 4 - 3 16 (-x + y) 4 .
Instead of considering incrementally the matrices in the Hankel family we can compute the polynomials P v and P w faster by applying Algorithm 4. For this, we consider the polynomial A := 5x 4 + 4x 3 + 3x 2 + 2x + 1, and the rows of the Extended Euclidean Algorithm for A and x 5 .

j V j U j R j 0 1 0 x 5 1 0 1 5x 4 + 4x 3 + 3x 2 + 2x + 1 2 1 1 25 (5x -4) 1 25 (x 3 + 2x 2 + 3x + 4) 3 -25(5x -6) 25(x -1) 2 25 4 1 25 (5x 4 + 4x 3 + 3x 2 + 2x + 1) -1 25 x 5 0
We need to consider the first j such that deg(R j ) < 5 2 , which is j = 3. Hence,

N 1 = max(deg(U 3 ) -1, deg(R 3 )) = 1 and P v := U 3 y x x max(deg(U 3 ),deg(R 3 )+1) = 25 y x -1 2 x 2 = 25(y -x) 2 . As deg(R 3 ) ≤ deg(U 3 ), we consider N 2 = deg(R 2 ) = 3 and P w := U 2 y x x deg(R 2 )+1 = 1 25 (5yx 3 -4x 4 ). ♦
The real case. When we consider the decomposition of binary forms over R, Algorithm 3 might fail. This happens either when the decomposition over C is not unique or when the decomposition over C is unique but it is not a decomposition over R. The algorithm fails because

• The real rank of the binary form might be bigger than N 2 + 1, see Reznick (2013a).

Lemma 17 does not hold over R and so we cannot find a square-free kernel polynomial Q µ that factors over R.

• Even when the real rank of the binary form is N 2 + 1, we need to perform some extra computations to compute a square-free kernel polynomial Q µ that factors over R. This computations are not taken into account in our algorithm, so we could never find such a decomposition.

Recently, Sylvester's algorithm was extended to the real case (Ansola et al., 2017, Algorithm 2).

This algorithm performs an incremental search over r as in Algorithm 1, but it decides if there is a real decomposition of length r by checking the emptiness of a semi-algebraic set. In the step r-th of the algorithm, the semi-algebraic set is embedded in R dim(Ker(H r )) . Hence, the bottleneck of their algorithm is not the computation of P v and P w as in our case, but deciding the emptiness of the semi-algebraic set. We emphasize that when the decompositions over C and R are unique and the same, our algorithm computes such a decomposition. Moreover, given P v , we can check if the previous condition holds by checking if P v splits over R.

Complexity

In this section we study the algebraic degree of the parameters that appear in the decomposition of a binary form as well as the arithmetic and bit complexity of Algorithm 3.

Algebraic degree of the problem

If we assume that the coefficients of the input binary form Equation (1) are rational numbers then the parameters of the decompositions, α j , β j , and λ j (see Equation ( 3)), are algebraic numbers, that is, roots of univariate polynomials with integer coefficients. The maximum degree of these polynomials is the algebraic degree of the problem. We refer the interested reader to Ba- Recall that, from Lemma 17, the rank of f could be either N 1 + 1 or N 2 + 1. When the polynomial P v is square-free, then the rank is N 1 + 1 and Q = P v . Following the discussion of Section 3.2, we prove that, when the rank of the binary form is N 2 + 1, we can compute a squarefree kernel polynomial Q of this degree such that the largest degree of its irreducible factors is N 1 . Moreover, we prove that for almost all the choices of (N 2 -N 1 + 1) different points in P 1 (K)

(the projective space of K) there is a square-free kernel polynomial of H N 2 +1 which vanishes on these points. This will be our choice for Q.

Lemma 25. Let f be a binary form of rank N 2 + 1. Given (N 2 -N 1 + 1) different points

(α 0 , β 0 ), . . . , (α N 2 -N 1 , β N 2 -N 1 ) ∈ P 1 (K)
such that none of them is a root of P v , then there is a unique binary form P µ of degree N 2 -N 1 , such that the kernel polynomial Q µ := P µ • P v + P w vanishes on these points.

Proof. Without loss of generality, we assume β i = 1, for i ∈ {1, . . . , N 2 -N 1 }. By Proposition 7, for any polynomial P µ of degree N 2 -N 1 , Q µ is a kernel polynomial. Since Q µ (α i , 1) = 0, we can interpolate P µ by noticing that P µ (α j , 1) = -P w (α j ,1) 1) . The degree of P µ is (N 2 -N 1 ) and we interpolate it using (N 2 -N 1 + 1) different points.

P v (α j ,
Hence there is a unique interpolation polynomial P µ . So, Q µ is the unique kernel polynomial of H N 2 +1 vanishing at all these points.

Example (cont.). For the example of Section 3.3, we obtained the square-free kernel polynomial Q by choosing the points (2, 1), (-2, 1) and (-1, 1) ∈ P 1 (K). If we choose other points such that Q is square-free, we will obtain a different decomposition. Hence, f does not have a unique decomposition. This holds in general.

♦ From Lemma 25, we deduce the following well-known result about the uniqueness of the decomposition, see also [START_REF] Helmke | Waring's problem for binary forms[END_REF]; [START_REF] Comas | On the rank of a binary form[END_REF]; [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF][START_REF] Carlini | Waring loci and the Strassen conjecture[END_REF].

Corollary 26. A decomposition is unique if and only if the rank is N 1 + 1 and N 1 < N 2 . A decomposition is not unique if and only if the rank is N 2 + 1.

Theorem 27. Let the rank of f be N 2 + 1. Then there is a square-free kernel polynomial Q such that the largest degree of its irreducible factors is at most N 1 .

Proof. If the rank of f is N 2 + 1, then for each set of N 2 -N 1 + 1 different points in P 1 (K), following the assumptions of Lemma 25, there is a unique kernel polynomial. There is a rational map that realizes this relation (see the proof of Lemma 25). Let this map be Q [α] , where can write this polynomial as

Q α 0 ,β 0 (x, y) = - P w (α 0 , β 0 ) P v (α 0 , β 0 ) M(x, y) M(α 0 , β 0 ) + N 2 -N 1 ∑ i=1 β 0 x -α 0 y α 0 β i -α i β 0 E i (x, y) P v (x, y) + P w (x, y)
Where E i (x, y) := -P w (α i ,β i ) P v (α i ,β i ) ∏ j {0,i} β j x-α j y α i β j -α j β i and M(x, y) := ∏ N 2 -N 1 j=1 (β j xα j y).2 For each (α j , β j ) ∈ S, we characterize the possible (α 0 , β 0 ) ∈ P 1 (K) such that (α j , β j ) is a root of Q α 0 ,β 0 of multiplicity bigger than one. Then, we study the (α 0 , β 0 ) ∈ P 1 (K) such that (α 0 , β 0 ) is a root of Q α 0 ,β 0 with multiplicity bigger than one. Finally, we reduce every case to the previous ones.

To study the multiplicities of the roots, we use the fact that (α 0 , β 0 ) is a double root of a binary form P if and only if

P(α 0 , β 0 ) = ∂ P ∂ x (α 0 , β 0 ) = ∂ P ∂ y (α 0 , β 0 ) = 0. Hence, for each (α 0 , β 0 ) ∈ P 1 (K), we consider ∂ Q α 0 ,β 0 ∂ x and ∂ Q α 0 ,β 0 ∂ y , where ∂ Q α 0 ,β 0 ∂ x = - P w (α 0 , β 0 ) P v (α 0 , β 0 ) 1 M(α 0 , β 0 ) ∂ M ∂ x P v + M ∂ P v ∂ x (x, y)+ (13) N 2 -N 1 ∑ i=1 1 β 0 α i -α 0 β i ∂ ((β 0 x -α 0 y)E i P v ) ∂ x (x, y) + ∂ P w ∂ x (x, y)
Let O α 0 ,β 0 x (x, y) be the product between the last line of Equation ( 13) and M(α 0 , β 0 ), that is

O α 0 ,β 0 x (x, y) := N 2 -N 1 ∑ i=1 M(α 0 , β 0 ) β 0 α i -α 0 β i ∂ ((β 0 x -α 0 y)E i P v ) ∂ x (x, y) + M(α 0 , β 0 ) ∂ P w ∂ x (x, y)
In what follows, instead of considering the pair (α 0 , β 0 ) as a point in P 1 (K), we consider it as a pair of variables. Hence, for every

(α i , β i ) ∈ S, (β 0 α i -α 0 β i ) divides M(α 0 , β 0 ), as polynomials in K[α 0 , β 0 ], so O α 0 ,β 0 x (x, y) is a polynomial in K[α 0 , β 0 ][x, y].
The derivative of Q α 0 ,β 0 with respect to x is a rational function in K(α 0 , β 0 )[x, y], that we can write as

∂ Q α 0 ,β 0 ∂ x = T α 0 ,β 0 (x,y) P v (α 0 ,β 0 )M(α 0 ,β 0 )
where

T α 0 ,β 0 (x, y) := -P w (α 0 , β 0 ) ∂ M ∂ x P v + M ∂ P v ∂ x (x, y) + O α 0 ,β 0 x (x, y)P v (α 0 , β 0 ) ∈ K[α 0 , β 0 ][x, y]
Lemma 30. For each (α i , β i ) ∈ S, there are at most N 2 + 1 possible ( ᾱ0 , β0 ) ∈ P 1 (K) such that ( ᾱ0 , β0 ) S, P v ( ᾱ0 , β0 ) 0 and that (α i , β i ) is a root of multiplicity bigger than 1 in Q ᾱ0 , β0 .

Proof. If (α i , β i ) is a root of multiplicity bigger than 1 in Q ᾱ0 , β0 , then ∂ Q ᾱ0 , β0 ∂ x (α i , β i ) = 0. Hence, we are looking for the ( ᾱ0 , β0 ) such that T ᾱ0 , β0 (α i , β i ) = 0. The polynomial T α 0 ,β 0 (α i , β i ) belongs to K[α 0 , β 0 ], so if it is not identically zero, then there are a finite number of points ( ᾱ0 , β0 ) ∈ P 1 (K) such that T ᾱ0 , β0 (α i , β i ) = 0. Moreover, the degree of the polynomial T α 0 ,β 0 (α i , β i ) is at most max(deg(P w ), deg(O

α 0 ,β 0 x (α i , β i ))+deg(P v )) = N 2 +1. Hence, if the polynomial is not zero, this finite number is at most N 2 + 1. The polynomial T α 0 ,β 0 (α i , β i ) ∈ K[α 0 , β 0 ] is not zero. Observe that as M is square-free, M(α i , β i ) = 0 and P v (α i , β i ) 0, then ∂ M ∂ x P v + M ∂ P v ∂ x (α i , β i ) 0.
Hence, as P w and P v are coprime, if ( ᾱ0 , β0 ) ∈ P 1 (K) and P v ( ᾱ0 , β0 ) = 0, then T ᾱ0 , β0 (α i , β i ) 0. That is, T α 0 ,β 0 (α i , β i ) does not vanish in the roots of P v .

Lemma 31. There are at most 2N 2 + 1 possible ( ᾱ0 , β0 ) ∈ P 1 (K) such that ( ᾱ0 , β0 ) S, P v ( ᾱ0 , β0 ) 0 and ( ᾱ0 , β0 ) is a root of multiplicity bigger than 1 in Q ᾱ0 , β0 .

Proof. Following the proof of Lemma 30, we study

T α 0 ,β 0 (α 0 , β 0 ) ∈ K[α 0 , β 0 ]. T α 0 ,β 0 (α 0 , β 0 ) = -P w (α 0 , β 0 ) ∂ M ∂ x P v + M ∂ P v ∂ x (α 0 , β 0 ) + O α 0 ,β 0 x (α 0 , β 0 )P v (α 0 , β 0 ) = -P w M ∂ P v ∂ x (α 0 , β 0 ) + O α 0 ,β 0 x -P w ∂ M ∂ x (α 0 , β 0 )P v (α 0 , β 0 )
Note that P w and P v are coprime. Also, M and P v are coprime. Hence, the polynomial T α 0 ,β 0 (α 0 , β 0 )

is not zero because P v does not divide P w M ∂ P v ∂ x . We conclude the proof by noting that the degree of T α 0 ,β 0 (α 0 , β 0 ) is bounded by 2N 2 + 1.

Lemma 32. Let ( ᾱ0 , β0 ), ( α0 , β0 ) ∈ P 1 (K) such that ( ᾱ0 , β0 ), ( α0 , β0 ) S, P v ( ᾱ0 , β0 ) 0. Hence, Q ᾱ0 , β0 ( α0 , β0 ) = 0 if and only if Q ᾱ0 , β0 (x, y) = Q α0 , β0 (x, y).

Proof. Assume that Q ᾱ0 , β0 ( α0 , β0 ) = 0. Following Lemma 25, we write Q ᾱ0 , β0 = P μ P v + P w and Q α0 , β0 = P μ P v + P w . As P v and P w are coprime and Q ᾱ0 , β0 ( α0 , β0 ) = 0, then P v ( α0 , β0 ) 0.

Consider Q ᾱ0 , β0 -Q α0 , β0 = P v (P μ -P μ ). This polynomial belongs to K[x, y] and it vanishes over P 1 (K) at the N 1 + 1 roots of P v , at the N 2 -N 1 points on S, and at ( α0 , β0 ) ∈ P 1 (K). Hence, Q ᾱ0 , β0 -Q α0 , β0 = 0 as it is a binary form in K[x, y] of degree at most N 2 + 1 with N 2 + 2 different roots over P 1 (K). Therefore, Q ᾱ0 , β0 (x, y) = Q α0 , β0 (x, y).

To prove the second case, note that, by definition, Q α0 , β0 ( α0 , β0 ) = 0. Hence, if we assume that Q ᾱ0 , β0 (x, y) = Q α0 , β0 (x, y), then we have Q ᾱ0 , β0 ( α0 , β0 ) = 0.

Proof of Proposition 29. We want to bound the number of different points ( ᾱ0 , β0 ) ∈ P 1 (K) such that Q ᾱ0 , β0 (x, y) is not a square-free binary form over K[x, y]. If the binary form Q ᾱ0 , β0 (x, y) is not square-free, then it has a root over P 1 (K) with multiplicity bigger than one. If such a root is (α i , β i ) ∈ S, we can bound the possible number of different values for ( ᾱ0 , β0 ) ∈ P 1 (K) by (N 2 + 1) (Lemma 31). Hence, if there is a i such that (α i , β i ) ∈ S has multiplicity bigger than one as a root of Q ᾱ0 , β0 (x, y), we can bound the possible number of different values for

( ᾱ0 , β0 ) ∈ P 1 (K) by #S • (N 2 + 1) = (N 2 -N 1 )(N 2 + 1).
If Q ᾱ0 , β0 is not square-free and the multiplicity of every root (α i , β i ) ∈ S is one, then there must be a root ( α0 , β0 ) ∈ P 1 (K) such that ( α0 , β0 ) S and its multiplicity as a root of Q ᾱ0 , β0 is bigger strictly than one. By Lemma 32, Q α0 , β0 (x, y) = Q ᾱ0 , β0 (x, y), and so ( α0 , β0 ) ∈ P 1 (K) has multiplicity bigger than one as a root of Q α0 , β0 (x, y). Hence, P v ( α0 , β0 ) 0 and, by Lemma 31, we can bound the possible number of different values for ( α0 , β0 ) ∈ P 1 (K) by 2N 2 + 1. As Q α0 , β0 (x, y) has N 1 + 1 roots over P 1 (K) \ S then, by Lemma 32, there are

N 1 + 1 different ( ᾱ0 , β0 ) ∈ P 1 (K) such that Q α0 , β0 (x, y) = Q ᾱ0 , β0 (x, y).
Hence, for each ( α0 , β0 ) ∈ P 1 (K) such that ( α0 , β0 ) has multiplicity bigger than one as a root of Q α0 , β0 (x, y), there are N 1 + 1 points ( ᾱ0 , β0 ) ∈ P 1 (K) such that ( α0 , β0 ) has multiplicity bigger than one as a root of Q ᾱ0 , β0 (x, y).

Therefore, the number of different values for ( ᾱ0 , β0 ) ∈ P 1 (K) such that Q ᾱ0 , β0 (x, y) has a root in P 1 (K) \ S with multiplicity bigger than one is bounded by (N 1 + 1)(2N 2 + 1).

Joining these bounds, we deduce that there are at most

(N 2 -N 1 )(N 2 + 1) + (N 1 + 1)(2N 2 + 1) different ( ᾱ0 , β0 ) ∈ P 1 (K) such that Q ᾱ0 , β0 is not square-free. Recalling that N 1 = D - N 2 and N 2 ≤ D (Proposition 4), we can bound (N 2 -N 1 )(N 2 + 1) + (N 1 + 1)(2N 2 + 1), by D 2 + 3D + 1.
We can also relate Proposition 29 to the Waring locus of the binary form f (x, y), see [START_REF] Carlini | Waring loci and the Strassen conjecture[END_REF]. For example, this proposition improves (Carlini et al., 2017, Proposition 3.8) by a factor of two. Moreover, it shows that the uniqueness condition from (Carlini et al., 2017, Proposition 3.8) misses some assumptions to hold3 .

Complexity of computing λ

We compute the coefficients λ j of the decomposition by solving a linear system involving a transposed Vandermonde matrix (Step 3 of Algorithm 3). We follow [START_REF] Kaltofen | Improved sparse multivariate polynomial interpolation algorithms[END_REF] to write the solution of Equation ( 5) as the evaluation of a rational function over the roots of a univariate polynomial.

Definition 33. Given a polynomial P(x) := ∑ n i=0 a i x i and 0 < k ≤ n, let

Quo(P(x), x k ) := n ∑ i=k a i x i-k .
Proposition 34 (Kaltofen and Yagati, 1989, Sec. 5). If α j α i , for all i j, then there is a unique solution to the system of Equation ( 14).

    1 1 • • • 1 α 1 α 2 • • • α r . . . . . . . . . α r-1 1 α r-1 2 • • • α r-1 r     λ =     a 0 a 1 . . . a r-1     (14) Moreover, if the solution is λ = (λ 1 , . . . , λ r ) T then, λ j = T Q (α j ) where Q (x) is the derivative of Q(x) := ∏ r i=1 (x -α i ), R(x) := ∑ r i=1 a r-i x i-1 and T (x) := Quo Q(x) • R(x), x r .
Lemma 35. Given a binary form f (x, y) := ∑ D i=0 D i a i x i y D-i , let Q be a square-free kernel polynomial of degree r, obtained after step 3 of Algorithm 3. Assume that y does not divide Q. Let α j be the j-th roots of Q(x), Q (x) be the derivative of Q(x) and the polynomial T (x) := Quo Q(x) • R(x), x r , with R(x) := ∑ r i=1 a r-i x i-1 . Then, each λ j from step 3 in Algorithm 4 can be written as λ j = T Q (α j ).

Proof. As y does not divide Q, we can write it as Q(x, y) = ∏ i (xα i y), where all the α i are different. Hence, as the r × r leading principal submatrix of Equation ( 5) is invertible, we can restrict the problem to solve that r × r leading principal subsystem. This system is Equation ( 14).

Therefore, the proof follows from Proposition 34.

Proposition-Definition 36 (Symbolic decomposition). Let Q be a square-free kernel polynomial related to a minimal decomposition of a binary form f of degree D, such that y does not divide Q. In this case, we can write f as

f (x, y) = ∑ {α∈K|Q(α)=0} T Q (α) • (αx + y) D .
which computes these rows in O(M(D) • log(D)). For a detailed reference of how this algorithm works see Bostan et al. (2017, Chapter 6.3) or Gathen and Gerhard (2013, Chapter 11).

Lemma 42 (Complexity of computing Q). Given the kernel polynomials P v and P w from Proposition 7, we compute a square-free polynomial Q µ := P µ • P v + P w such that the maximal degree of its irreducible factors is bounded by Theorem 28 in O(M(D) • log(D)).

Proof. To compute the vector µ, we choose randomly N 2 -N 1 + 1 linear forms in K[x, y] and we proceed as in Lemma 25. The complexity bound is due to multi-point evaluation and interpolation of a univariate polynomial (Gathen and Gerhard, 2013, Chapter 10).

Theorem 43. When the decomposition is unique, that is when the rank is N 1 + 1, then Algorithm 3 computes deterministically a symbolic decomposition (Proposition-Definition 36) of a binary form in O(M(D) log(D)).

When the decomposition is not unique, that is when the rank is N 2 +1 and N 1 < N 2 , then Algorithm 3 is a Monte Carlo algorithm that computes a symbolic decomposition of a binary form in O(M(D) log(D)).

Proof. The first step of the algorithm, in both cases, is to compute the kernel polynomials P v and P w of Proposition 7 using Algorithm 4. By Lemma 41, we compute them deterministically in O(M(D) • log(D)).

If P v is square-free, which means that the decomposition is unique, then Q = P v . Otherwise, in step 2, we need to choose some random values to construct the square-free polynomial Q (from the kernel polynomials P v and P w ) in O(M(D) • log(D)) as in Lemma 42. This is the step that makes the algorithm a Monte Carlo one, as we might fail to produce a square-free polynomial Q.

In both cases, at step 3 we compute the rational function that describes the solution of the system in Equation ( 5 all the roots up to accuracy 2 -1 in O B (D 2 σ + D 1 ) (Pan and Tsigaridas, 2017a). That is, we compute absolute approximations of α j , say α j , such that |α jα j | ≤ 2 -1 .

The next step consists in solving the (transposed) Vandermonde system, V ( α) T λ = a, where V ( α) is the Vandermonde matrix we construct with the roots of Q, λ is a vector contains the coefficients of decomposition, and a is a vector containing the coefficients of F, see also Equation (5). We know the entries of V ( α) up to 1 bits. Therefore, we can compute the elements To estimate the accuracy of f we need to expand the approximate decomposition and consider it as a polynomial in x. We do not actually perform this operation; we only estimate the accuracy as if we were. First, we expand each ( α j x + y) D . This results polynomials with coeffi- and Tsigaridas, 2017b, Lemma 19). Next, we multiply each such polynomial with λ i , and we collect the coefficients for the various powers of x. Each coefficient is the sum of r ≤ D terms.

The last two operations do not affect, asymptotically, the precision. Therefore, the polynomial f = ∑ r j=1 λ j ( α j x + (1α j t)y) D that corresponds to the approximate decomposition has an absolute approximation such that ff ≤ 2 -1 +O(D lg(D)σ ) . To achieve an accuracy of 2 -in the decomposition, such that ff ≤ 2 -, we should choose

  tion 4.1.3). For polynomials with integer coefficients, we bound the bit complexity, up to polylogarithmic factors, by O B (D + D 4 + D 3 τ), where τ is the maximum bitsize of the coefficients of the input binary form and 2 -is the error of the terms in the decomposition (Theorem 45).

Figure 1 :

 1 Figure 1: Relation between N 1 , N 2 , and D

  as we want to prove. By assumption, deg(R)< N 2 ≤ D and deg(U v Q) = deg(U) ≤ N 2 .By definition, the degree of R v is smaller of equal to N 1 , and so

Lemma 20 .

 20 If N 1 = N 2 , then there is a unique rational decomposition (U, R) such that deg(U)≤ D 2 , deg(R) ≤ D2 and R is monic. In particular, deg(R) = D 2 and we can consider the kernel polynomial related to v ∈ Ker(H N 1 +1 ) (Proposition 6) as P v = U

  jaj (1988);[START_REF] Nie | The algebraic degree of semidefinite programming[END_REF];[START_REF] Draisma | The euclidean distance degree of an algebraic variety[END_REF] for a detailed exposition about the algebraic degree and how it address the complexity of the problem at hand at a fundamental level.4.1.1. The complexity of computing Q

  ) in O(M(D) • log(D)) (Lemma 35). At step 4 of the algorithm, we return the decomposition.We can bound the probability of error of Algorithm 3 using Proposition 29, which bounds the number of bad values that lead us to a non square-free polynomial Q. Moreover, we can introduce a Las Vegas version of Algorithm 3 by checking if the values that we choose to construct a polynomial Q result indeed a square-free polynomial. We can do this check in O(M(D) • log(D)), by computing the GCD between the Q and its derivatives.

  of the solution vector λ with an absolute approximation correct up to2 = 1 -O(D lg(D)σlg ∏ j ∆(α j )) = 1 -O(D lg(D)σ ) bits(Pan and Tsigaridas, 2017b, Theorem 29). That is, we compute λ j 's such that |λ jλ j | ≤ 2 -2 . At this point we have obtained the approximate decomposition f (x, y) := r ∑ j=1 λ j ( α j x + y) D .

  cients correct up 3 = 2 -O(Dσ ) = 1 -O(D lg(D)σ ) -O(Dσ ) = 1 -O(D lg(D)σ ) bits (Pan

  1 = + O(D lg(D)σ ). Thus, all the computations should be performed with precision of + O(D lg(D)σ ) bits. The bit complexity of computing the decomposition of f up to bits is dominated by the solving and refining process and it is O B (D + D 2 σ ). If we substitute the value for σ , then we arrive at the complexity bound of O B (D + D 4 + D 3 τ). Theorem 45. Let f ∈ Z[x, y] be a homogeneous polynomial of degree D and maximum coefficient bitsize τ. We compute an approximate decomposition of accuracy 2 -in O B (D + D 4 + D 3 τ) bit operations.

+1)×( D 2 +1) is square and, as Ker(H N 1 ) = 0 (Proposition 6), it is invertible. If r N 1 = 0, that is deg(R v ) < N 1 ,then the polynomial U is zero and so (U, R) is not a

For each 0 ≤ i ≤ N 2 -N 1 , Q α 0 ,β 0 (x, y) is a rational function of degree 0 with respect to (α i , β i ). Hence, it is well defined the evaluation of the variables (α i , β i ) in Q α 0 ,β 0 (x, y) at points of P 1 (K).

The authors are not taking into account that the lambdas that they use in their proof are not unique, and so they give us more degrees of freedom that we can use to fix more terms in the decomposition.

Acknowledgments

The authors thank George Labahn and Vincent Neiger for pointing out important references to derandomize the computation of P v and P w . The authors thank the anonymous reviewers for their helpful comments. Matías Bender thanks Joos Heintz for supervising his Master's thesis.

The authors are partially supported by ANR JCJC GALOP (ANR-17-CE40-0009), the PGMO grant ALMA and the PHC GRAPE.

Algorithm 4 COMPUTE PV AND PW

Input: A sequence (a 0 , . . . , a D ).

Output: Polynomials P v and P w as 7.

1. i ← first row of EGCD(x D+1 , ∑ i=0 a i x i ) such that R i < D+1 2 .

2.

.

N 2 ← deg(R i-1 ).

4. Return P v and P w 3.2. Computing a square-free polynomial Q

We can compute Q at step 2 of Algorithm 3 in different ways. If P v is square-free, then we set Q equal to P v . If P v is not square-free, by Lemma 17, we need to find a vector µ ∈ K (N 2 -N 1 +1)

such that Q µ := P µ • P v + P w is square-free. By Proposition 7, P v and P w are relative prime. Thus, if we take a random vector µ, generically, Q µ would be square-free. For this to hold, we have to prove that the discriminant of Q µ is not identically zero. To simplify notation, in the following theorem we dehomogenize the polynomials.

Theorem 22. Given two relative prime univariate polynomials P v (x) and P w (x) of degrees N 1 +1

and N 2 + 1 respectively, let Q µ (x) := P µ (x, 1)

Proof. The zeros of the discriminant of Q µ (x) with respect to x over K correspond to the set {µ ∈ K N 2 -N 1 +1 : Q µ has double roots}. We want to prove that the discriminant is not zero.

A univariate polynomial is square-free if and only if it does share any root with its deriva-

. The image of the map is contained in

Given a kernel polynomial Q(x, y), there is a finite number of distinct points (α, β ) ∈ P 1 (K)

such that Q(α, β ) = 0. Hence, the pre-image of an element in the image of

Therefore, the dimension of the image and the dimension of the domain are the same.

By Theorem 22, the non-square-free kernel polynomials form a hypersurface in the space of kernel polynomials of the shape P µ • P v + P w . If we consider the pre-image of the intersection between this hypersurface and the image of the rational map, then its dimension is smaller than

Therefore, generically, for

a square-free kernel polynomial. As K is the algebraic closure of K ⊂ C, the same holds over K.

Theorem 28. Given a binary form f of rank r and degree D, there is a square-free kernel polynomial of degree r such that the biggest degree of its irreducible factors is min(r, Dr + 1).

Proof. If the rank is r = N 2 + 1, then min(r, Dr + 1) = N 1 . By Theorem 27, such a squarefree kernel polynomial exists. If the rank is r = N 1 + 1 and N 1 < N 2 , by Lemma 17, there is a square-free kernel polynomial of degree min(r,

The previous result is related to the decomposition of tensors of the same border rank (Comas and Seiguer, 2011, Theorem 2); (Bernardi et al., 2011, Theorem 23); [START_REF] Blekherman | Typical real ranks of binary forms[END_REF].

We can also bound the number of possible bad choices in the proof of Theorem 27.

Proposition 29. Let f be a binary form of rank N 2 + 1. For every set S ⊂ P 1 (K) of cardinal

0 and the unique kernel polynomial Q µ := P µ • P v + P w that vanish over S and (α 0 , β 0 ) (Lemma 25) is not square-free.

To prove this proposition we use Lagrange polynomials to construct the maps and varieties of the proof of Theorem 27.

be the set of Proposition 29. For each (α 0 , β 0 ) ∈ P 1 (K) such that (α 0 , β 0 ) S and P v (α 0 , β 0 ) 0 we consider the unique kernel polynomial Q α 0 ,β 0 which vanishes at S and (α 0 , β 0 ), see Lemma 25. Using Lagrange polynomial, we Remark 37. If the square-free kernel polynomial related to a decomposition of rank r is divisible by y, we can compute {λ j } j<r of Equation (5) as in Lemma 35, by taking Q y as the kernel polynomial. It is without loss of generality to consider Q = P (u 0 ,...,-1,0) T , because Q is square-free, and so y 2 can not divide it. Hence, λ r = a D -∑ r-2 i=0 u i a D-r+i+1 (Reznick, 2013a, Equation 2.12).

To summarize the section, given a binary form f of rank r, there is a square-free kernel polynomial Q of the degree r, such that the largest degree of its irreducible factors is bounded by

where T and Q are polynomials whose coefficients belong to K and whose degrees are bounded by r, defined in Lemma 35. When y divides Q, the form is similar.

Lower bounds on the algebraic degree

In this section we analyze the tightness of the bound of Theorem 28. To do so, we construct families of examples where the bound is tight. We present two families of examples. In the first one, the decomposition is unique. In the second one, it is not.

Proposition 38 (Heinig and Rost, 1984, Theorem 5.2). For every pair of relatively prime binary forms, Pv and Pw , of degrees N1 + 1 and N2 + 1, N1 ≤ N2 , there is a sequence a = (a 0 , . . . , a N1 + N2 )

such that N a 1 = N1 , N a 2 = N2 , and we can consider the polynomials Pv and Pw as the kernel polynomials P v and P w from Proposition 7 with respect to the family of Hankel matrices {H k a } k .

Corollary 39. If there is an irreducible binary form of degree N1 + 1 in K[x, y], then for every D > 2 N1 , there is a binary form f of degree D such that its decomposition is unique, its rank N1 + 1, and the degree of the biggest irreducible factor of the polynomial Q from Algorithm 3 in the decomposition is min( N1 + 1, D -N1 ) = N1 + 1. That is, the algebraic degree of the minimal decomposition over K is N1 + 1 and the bound of Theorem 28 is tight.

Proof. Let Pv be a irreducible binary form of degree N1 + 1. Let Pw be any binary form of degree N2 + 1 := D -N1 + 1 relatively prime with Pv . Consider the sequence a = (a 0 , . . . , a N1 + N2 ) of Proposition 38 with respect to Pv and Pw , and the binary form f (x, y) := ∑ D i=0 D i a i x i y D-i . As K is of characteristic zero, K is a perfect field, and so, as Pv is irreducible, it is square-free. Then, by Lemma 17, the rank of the decomposition is N a 1 + 1 = N1 + 1, and by Corollary 26 the decomposition is unique. Following Algorithm 3, the polynomial Q is equal to Pv , which is an irreducible polynomial of degree N1 + 1. As D > 2 N1 , then min( N1 + 1, D -N1 ) = N1 + 1 and the bound of Theorem 28 is tight.

Lemma 40. Let K = Q and p ∈ N a prime number. Then, there is a binary form f of degree 2(p -1) whose decomposition is not unique and the bound of Theorem 28 is tight.

Proof. Consider the binary form f (x, y) := 2(p-1) p-1 x p-1 y p-1 . Using Algorithm 4, we obtain P v = -y p and P w = x p , N 1 = N 2 = p -1. The polynomial P v is not square-free, so we have to consider a square-free kernel polynomial in Ker(H N 2 +1 ). Moreover, the rank of the decomposition is N 2 + 1 = p. Every kernel polynomial in Ker(H N 2 +1 ) in Q[x, y] can be written as µ w x pµ v y p for some µ w , µ v ∈ Q. We are interested in the zeros of these polynomials (step 3 of Algorithm 3), thus we consider coprime µ w , µ v ∈ Z, as the zeros do not change. As we want to consider squarefree kernel polynomials, neither µ w nor µ v can be zero, and so (1, 0) ∈ P 1 (Q) is not a root of any of these polynomials. Hence, we rewrite our polynomial as 1 µ v y p ( µ w x p µ v y p -1), and so we look for the factorization over Q[z] of µ w µ v z p -1, where z = x y . We can use the Newton's polygon criterion, e.g., Cassels (1986, Chapter 6.3)

y] and so the degree of its biggest irreducible factor is p > min(p, 2 (p -1) -p + 1) = p -1 (Theorem 28). If this is not the case, then p | µ w µ v | ∈ Q, and so we can factor it as

The second factor is irreducible because there is an automorphism in Q[x] (given by z → p | µ v µ w |z) that transforms it into the p-th cyclotomic polynomial, which is irreducible as p is prime. Hence, the biggest irreducible factor of this polynomial has degree p -1 = min(p, 2 (p -1) -p + 1) and the bound of Theorem 28 is tight.

Arithmetic complexity

Lemma 41 (Complexity of Algorithm 4). Given a binary form f = ∑ D i=0 D i a i x i y D-i of degree D, Algorithm 4 computes P v and P w in O(M(D) • log(D)) arithmetic operations.

Proof. The complexity of the algorithm is the complexity of computing the rows (i + 1), i and (i -1) of the Extended Euclidean algorithm between ∑ i=0 a i x i and x D+1 , where the i-th row is the first row i such that deg(R i ) < D+1 2 (Lemma 21). This can be done using the Half-GCD algorithm,

Remark 44. If we want to output an approximation of the terms of the minimal decomposition, with a relative error of 2 -ε , we can use Pan's algorithm (Pan, 2002) (McNamee andPan, 2013, Theorem 15.1.1) to approximate the roots of Q. In this case the complexity becomes O D log 2 (D) log 2 (D) + log(ε) .

Bit complexity

Let f ∈ Z[x, y] be a binary form as in Equation ( 1), of degree D and let τ be the maximum bitsize of the coefficients a i . We study the bit complexity of computing suitable approximations of the α j 's, β j 's, and λ j 's of Equation ( 3), say α j , β j and λ j respectively, that induce an approximate decomposition correct up to bits. That is a decomposition such that f -∑ j λ j ( α j x + β j y) D ∞ ≤ 2 -. We need to estimate an upper bound on the number of bits that are necessary to perform all the operations of the algorithm.

The first step of the algorithm is to compute P v and P w , via the computation of three rows of the Extended GCD of two polynomials of degree D and D + 1 with coefficients of maximal sized τ. This can be achieved in O B (D 2 τ) bit operations (Gathen and Gerhard, 2013, Corollary 11.14.B), and the maximal bit size of P v and P w is O(Dτ). We check if P v is a square-free polynomial in O(D 2 τ), via the computation of the GCD of P v (x, 1) and its derivative (Gathen and Gerhard, 2013, Corollary 11.14.A), and by checking if y 2 divides it.

If P v is square-free polynomial, then Q = P v . If P v is not square-free, then we can compute Q by assigning values to the coefficients of P µ . We assume that y 2 does not divide P w , if this does not hold, we replace P w by the kernel polynomial x N 2 -N 1 P v + P w , which is coprime to P v , and so not divisible by y, as P v and the original P w are coprime (Proposition 7). We set all the coefficients of P µ to zero, except the constant term. Then Q = µ 0 y N 2 -N 1 P v + P w . Now we have to choose µ 0 so that Q is square-free. As y N 2 -N 1 P v and P w are coprime, there are at most 2D + 2 forbidden values for µ 0 such that Q is not square-free (Corollary 23), thus at least one of the first 2D + 3 integer fits our requirements. We test them all. Each test corresponds to a GCD computation, that costs O B (D 2 τ) and so the overall cost is O B (D 3 τ).

Let σ = O(Dτ) be the maximal bit size of Q. By Remark 37, we can assume that y does not divide Q, consider y = 1 and treat Q as an univariate polynomial.

Let {α j } j be the roots of Q. We isolate them in O B (D 2 σ ) [START_REF] Pan | Univariate polynomials: Nearly optimal algorithms for numerical factorization and root-finding[END_REF]. For the (aggregate) separation bound of the roots it holds thatlg ∏ j ∆(α j ) = O(Dσ + D lg(D)). We approximate