Jérémy Berthomieu
email: jeremy.berthomieu@lip6.fr

Jean-Charles Faugère
email: jean-charles.faugere@inria.fr

Ludovic Perret
email: ludovic.perret@lip6.fr

Polynomial-Time Algorithms for Quadratic Isomorphism of Polynomials: The Regular Case

Keywords: Quadratic forms, computer algebra, polynomial isomorphism, multivariate cryptography, module isomorphism 2010 MSC: 12Y05, 94A60, 68W20, 68W30, 68Q25

Let f = ( f 1 , . . . , f m ) and g = (g 1 , . . . , g m ) be two sets of m ≥ 1 nonlinear polynomials in K[x 1 , . . . , x n ] (K being a field). We consider the computational problem of finding -if any -an invertible transformation on the variables mapping f to g. The corresponding equivalence problem is known as Isomorphism of Polynomials with one Secret (IP1S) and is a fundamental problem in multivariate cryptography. Amongst its applications, we can cite Graph Isomorphism (GI) which reduces to equivalence of cubic polynomials with respect to an invertible linear change of variables, according to Agrawal and Saxena. The main result is a randomized polynomial-time algorithm for solving

IP1S for quadratic instances -a particular case of importance in cryptography.

To this end, we show that IP1S for quadratic polynomials can be reduced to a variant of the classical module isomorphism problem in representation theory. We show that we can essentially linearize the problem by reducing quadratic-IP1S to test the orthogonal simultaneous similarity of symmetric matrices; this latter problem was shown by Chistov, Ivanyos and Karpinski (ISSAC 1997) to be equivalent to finding an invertible matrix in the linear space K n×n of n × n matrices over K and to compute the square root in a certain representation in a matrix algebra. While computing square roots of matrices can be done efficiently using numerical methods, it seems difficult to control the bit complexity of such methods. However, we present exact and polynomial-time algorithms for computing a representation of the square root of a matrix in K n×n , for various fields (including finite fields), as a product of two matrices. Each coefficient of these matrices lie in an extension field of K of polynomial degree. We then consider #IP1S, the counting version of IP1S for quadratic instances. In particular, we provide a (complete) characterization of the automorphism group of homogeneous quadratic polynomials. Finally, we also consider the more general Isomorphism of Polynomials (IP) problem where we allow an invertible linear transformation on the variables and on the set of polynomials. A randomized polynomial-time algorithm for solv-

Introduction

A fundamental question in computer science is to provide algorithms allowing to test if two given objects are equivalent with respect to some transformation. In this paper, we consider equivalence of nonlinear polynomials in several variables. Equivalence of polynomials has profound connections with a rich variety of fundamental problems in computer science, ranging -among others topics -from cryptography (e.g. [START_REF] Patarin | Hidden Fields Equations (HFE) and Isomorphisms of Polynomials (IP): Two New Families of Asymmetric Algorithms[END_REF]; Tang and[START_REF] Tang | Towards provably secure proxy signature scheme based on isomorphisms of polynomials[END_REF]; [START_REF] Yang | A novel group signature scheme based on mpkc[END_REF]), arithmetic complexity (via Geometric Complexity Theory (GCT) for instance, see [START_REF] Bürgisser | Prospects for geometric complexity theory[END_REF]; [START_REF] Kayal | Affine projections of polynomials: extended abstract[END_REF]; [START_REF] Mulmuley | The GCT program toward the P vs. NP problem[END_REF]; [START_REF] Mulmuley | Geometric Complexity Theory I: An Approach to the P vs. NP and Related Problems[END_REF]), testing low degree affine-invariant properties [START_REF] Bhattacharyya | Testing low complexity affine-invariant properties[END_REF]; [START_REF] Green | The distribution of polynomials over finite fields, with applications to the gowers norms[END_REF]; [START_REF] Grigorescu | Tight lower bounds for testing linear isomorphism[END_REF], . . .). As we will see, the notion of equivalence can come with different flavours that impact the intrinsic hardness of the problem considered.

In [START_REF] Agrawal | Equivalence of F-Algebras and Cubic Forms[END_REF]; [START_REF] Saxena | Morphisms of Rings and Applications to Complexity[END_REF], the authors show that Graph Isomorphism reduces to equivalence of cubic polynomials with respect to an invertible linear change of variables (a similar reduction holds between F-algebra Isomorphism and cubic equivalence of polynomials). This strongly suggests that solving equivalence problems efficiently is a very challenging algorithmic task.

In cryptography, the hardness of deciding equivalence between two sets of m polynomials with respect to an invertible linear change of variables is the security core of several cryptographic schemes: the seminal zero-knowledge ID scheme of [START_REF] Patarin | Hidden Fields Equations (HFE) and Isomorphisms of Polynomials (IP): Two New Families of Asymmetric Algorithms[END_REF], and more recently group/proxy signature schemes (Tang and[START_REF] Tang | Towards provably secure proxy signature scheme based on isomorphisms of polynomials[END_REF]; [START_REF] Yang | A novel group signature scheme based on mpkc[END_REF]). Note that there is a subtle difference between the equivalence problem considered in [START_REF] Agrawal | Equivalence of F-Algebras and Cubic Forms[END_REF]; [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF]; [START_REF] Saxena | Morphisms of Rings and Applications to Complexity[END_REF] and the one considered in cryptographic applications.

Whilst [START_REF] Agrawal | Equivalence of F-Algebras and Cubic Forms[END_REF]; [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF]; [START_REF] Saxena | Morphisms of Rings and Applications to Complexity[END_REF] restrict their attention to m = 1, arbitrary m ≥ 1 is usually considered in cryptographic applications. In the former case, the problem is called Polynomial Equivalence (PolyEquiv), whereas it is called Isomorphism of Polynomials with One Secret (IP1S) problem in the latter case. We emphasize that the hardness of equivalence can drastically vary in function of m. An interesting example is the case of quadratic forms. The problem is completely solved when m = 1, but no polynomial-time algorithm exists for deciding simultaneous equivalence of quadratic forms. In this paper, we make a step ahead to close this gap by presenting a randomized polynomial-time algorithm for solving simultaneous equivalence of quadratic forms over various fields. Indeed, by slicing of the polynomials degree by degree, one can find equations in the coefficients of the transformation allowing one to recover the transformation. More recently, [START_REF] Bouillaguet | Graph-theoretic algorithms for the "isomorphism of polynomials" problem[END_REF] presented exponential (in the number of variables n) algorithms for solving quadratic homogeneous instances of IP over finite fields. This situation is clearly unsatisfying, and suggests that an important open problem for IP is to identify large class of instances which can be solved in (randomized) polynomial time.

An important special case of IP is the IP problem with one secret (IP1S for short), where B is the identity matrix. From a cryptographic point of view, the most natural case encountered for equivalence problems is inhomogeneous polynomials with affine transformations. For IP1S, we show that such a case can be handled in the same way as homogeneous instances with linear transformations (see Proposition 5). As a side remark, we mention that there exist more efficient methods to handle the affine case; typically by considering the homogeneous components, see [START_REF] Faugère | Polynomial Equivalence Problems: Algorithmic and Theoretical Aspects[END_REF]. However, homogenizing the instances allows us to make the proofs simpler and cleaner. As such, we focus our attention to solve IP1S for quadratic homogeneous forms.

When m = 1, the IP1S problem can be easily solved by computing a reduced form of the input quadratic forms. In [START_REF] Bouillaguet | Practical cryptanalysis of the identification scheme based on the isomorphism of polynomial with one secret problem[END_REF], the authors present an efficient heuristic algorithm for solving IP1S on quadratic instances. However, the algorithm requires to compute a Gröbner basis. So, its complexity could be exponential in the worst case. More recently, [START_REF] Lidl | New Insight into the Isomorphism of Polynomial Problem IP1S and Its Use in Cryptography[END_REF] proposed a polynomial-time algorithm for solving IP1S on quadratic instances with m = 2 over fields of any characteristic. We consider here arbitrary m > 1.

In computer algebra, a fundamental and related problem is the simplification of a homogeneous polynomial system f ∈ K[x] m . That is, compute A ∈ GL n (K) such that g(x) = f(A • x) is easier to solve. In this setting, RIDGE algorithm (see [START_REF] Berthomieu | Computing Hironaka's invariants: Ridge and Directrix[END_REF]; [START_REF] Hironaka | Additive groups associated with points of a projective space[END_REF]; [START_REF] Giraud | Étude locale des singularités. U.E.R. Mathématique[END_REF]) and MINVAR algorithm (see [START_REF] Carlini | Reducing the number of variables of a polynomial[END_REF]; [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF]) reduce to the best the number of variables of the system. More generally, for a given homogeneous polynomial system f, the Functional Decomposition Problem is the problem of computing h = (h 1 , . . . , h s ) homogeneous and g such that f(x) = g(h(x)).

To simplify the presentation in this introduction, we mainly deal with fields of characteristic 2. Results for fields of characteristic 2 are also given later in this paper. Now, we define formally IP1S:

Definition 1. Let f = ( f 1 , . . . , f m ), g = (g 1 , . . . , g m ) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . , x n ] m . We shall say that f and g are equivalent, denoted f ∼ g, if ∃ A ∈ GL n (K) such that: g(x) = f(A • x). IP1S is then the problem of finding -if any -A ∈ GL n (K) that makes g equivalent to f (i.e. A ∈ GL n (K) such that g(x) = f(A • x) .
In such a case, we present a randomized polynomial-time algorithm for solving IP1S with quadratic polynomials. To do so, we show that such a problem can be reduced to the variant of a classical problem of representation theory over finite dimensional algebras. In our setting we need, as in the case m = 1, to provide a canonical form of the problem.

Canonical Form of IP1S Let f = ( f 1 , . . . , f m ), g = (g 1 , . . . , g m ) ∈ K[x 1 , . . . , x n ] m ×K[x 1 , . . . ,
x n ] m be homogeneous quadratic polynomials. Let H 1 , . . . , H m be the Hessian matrices of f 1 , . . . , f m (resp. H ′ 1 , . . . , H ′ m be the Hessian matrices of g 1 , . . . , g m ). Recall that the Hessian matrix associated to a f i is defined as

H i = ∂ 2 f i ∂ x k ∂ x ℓ k,ℓ ∈ K n×n . Consequently, IP1S for quadratic forms is equivalent to finding A ∈ GL n (K) such that: H ′ i = A T • H i • A, for all i, 1 ≤ i ≤ m. (1)
Assuming H j is invertible, and thus so is

H ′ j , one has H ′-1 j = A -1 H -1 j A -T . Combining this with equation (1) yields H ′-1 j H ′ i = A -1 • H -1 j H i • A.
If none of the H i 's is invertible, then we look for an invertible linear combination thereof. For this reason, we assume all along this paper:

Assumption 1 (Regularity assumption). Let f = ( f 1 , . . . , f m ) ∈ K[x 1 , . . . , x n ].
We assume that a linear combination over K of the quadratic forms f 1 , . . . , f m is not degenerate1 . In particular, we assume that |K| > n.

Taking as variables the entries of A, we can see that (1) naturally yields a nonlinear system of equations. However, we show that one can essentially linearize equations (1). To this end, we prove in Section 2 that under Assumption 1 any quadratic homogeneous instance IP1S can be reduced, under a randomized process, to a canonical form on which -in particular -all the quadratic forms are nondegenerate. We shall call these instances regular. More precisely: Theorem 1. Let K be a field of char K 2. There exists a randomized polynomial-time algorithm which given a regular quadratic homogeneous instance of IP1S returns "NOSOLUTION" only if the two systems are not equivalent or a canonical form

n ∑ i=1 d i x 2 i , f 2 , . . . , f m , n ∑ i=1 d i x 2 i , g 2 , . . . , g m ,
where the d i are equal to 1 or a nonsquare in K, f i and g i are nondegenerate homogeneous quadratic polynomials in K[x 1 , . . . , x n ]. Any solution on K on the canonical form can be efficiently mapped to a solution of the initial instance (and conversely).

Let us note that over the rationals, computing the exact same sum of squares for the first quadratic forms of each set is difficult, see (Saxena, 2006, Chapter 3), (Wallenborn, 2013, Chapter 1). As such, one could only assume that the first quadratic form of the second set is

∑ n i=1 d ′ i x 2 i .
This does not fundamentally change the algorithms presented in this paper, beside some matrices denoted by D which could be changed into

D ′ = Diag(d ′ 1 , . . . , d ′ n ).
Note that the success probability of the algorithms presented here will depend on the size of the field. If one looks for A ∈ L n×n with L an extension of K, one can amplify the success probability over a small field by using the fact that matrices are conjugate over K if and only if they are conjugate over an algebraic extension L (see [START_REF] De Seguins Pazzis | Invariance of simultaneous similarity and equivalence of matrices under extension of the ground field[END_REF]). Thus, one can search linear change of variables with coefficients in some algebraic extension L ⊇ K (but of limited degree).

Conjugacy Problem

When IP1S is given in canonical form, equations (1) can be rewritten as

A T D A = D with D = Diag(d 1 , . . . , d n ) and H ′ i = A T • H i • A = D A -1 D -1 • H i • A for all i, 2 ≤ i ≤ m.
Our task is now to solve the following problem:

Definition 2 (D-Orthogonal Simultaneous Matrix Conjugacy (D-OSMC)). Let K n×n be the set of n × n matrices with entries in K. Let {H 1 , . . . , H m } and {H ′ 1 , . . . , H ′ m } be two families of matrices in K n×n . The D-OSMC problem is the task to recover -if any -a D-orthogonal matrix X ∈ L n×n , i.e. X T DX = D, with L being an algebraic extension of K, such that: [START_REF] Chistov | Polynomial time algorithms for modules over finite dimensional algebras[END_REF] show that D-OSMC with D = Id is equivalent to:

X -1 H i X = H ′ i , ∀ i, 1 ≤ i ≤ m,
1. Solving the Simultaneous Matrix Conjugacy problem (SMC) between {H i } 1≤i≤m and {H ′ i } 1≤i≤m , that is to say finding an invertible matrix Y ∈ GL n (K) such that:

Y -1 • H i •Y = H ′ i and Y -1 • H T i •Y = H ′ i T ∀ i, 1 ≤ i ≤ m.
(2)

2. Computing the square-root W of the matrix Z = Y • Y T . Then, the solution of the D-OSMC problem is given by X = Y W -1 .

In our context, D = Diag(d 1 , . . . , d n ) is any diagonal invertible matrix. So, we extend [START_REF] Chistov | Polynomial time algorithms for modules over finite dimensional algebras[END_REF] and show that D-OSMC is equivalent to 1. Finding an invertible matrix Y ∈ GL n (K) such that:

Y -1 • H i •Y = H ′ i and DY -1 D -1 • H T i • DY D -1 = H ′ i T ∀ i, 1 ≤ i ≤ m. (3) 2. Computing the square-root W of the matrix Z = DY • Y T D -1 . Then, the solution of the D-OSMC problem is given by X = Y W -1 .
In our case, the H i 's (resp. H ′ i 's) are symmetric (Hessian matrices). Thus, condition (3) yields a system of linear equations and one polynomial inequation:

H 1 •Y = Y • H ′ 1 , . . . , H m •Y = Y • H ′ m and det(Y ) 0. ( 4
)
From now on, we shall denote by O n (L, D) the set of D-orthogonal matrices with coefficients in L. Let V ⊂ K n×n be the linear subspace of matrices defined by these linear equations. The SMC problem is then equivalent to recovering an invertible matrix in V ; in other words we have to solve a particular instance of Edmonds' problem (Edmonds (1967)). Note that, if the representation of the algebra spanned by {H -1 1 H i } 1≤i≤m is irreducible, we know that V has dimension at most 1 (Schur's Lemma, see (Lang, 2002, Chap. XVII, Proposition 1.1) and (Newman, 1967, Lemma 2) for a matrix version of this lemma). After putting the equations in triangular form, randomly sampling over the free variables an element in V yields, thanks to Schwartz-Zippel-DeMillo-Lipton Lemma [START_REF] Demillo | A probabilistic remark on algebraic program testing[END_REF]; [START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF]), a solution to D-OSMC with overwhelming probability as soon as K is big enough. If one accepts to have a solution matrix over an extension field, we can amplify the probability of success by considering a bigger algebraic extension (see de Seguins Pazzis ( 2010)). Whilst a rather "easy" randomized polynomial-time algorithm solves SMC, the task of finding a deterministic algorithm is more delicate. In our particular case, we can adapt the result of [START_REF] Chistov | Polynomial time algorithms for modules over finite dimensional algebras[END_REF] and provide a deterministic polynomial-time algorithm for solving (2).

Characteristic 2

Let us recall that in characteristic 2, the associated matrices H 1 , . . . , H m , H ′ 1 , . . . , H ′ m to quadratic forms can be chosen as upper triangular. In this context, we show in Section 3.4 that IP1S can still be reduced to a (H 1 + H T 1 )-conjugacy problem. Under certain conditions in even dimension, we can solve this conjugacy problem in polynomial-time. These results are well confirmed by some experimental results presented in Section 3.5. We can recover a solution in less than one second for n up to one hundred (cryptographic applications of IP1S usually require smaller values of n, typically ≤ 30, for efficiency reasons).

Matrix Square Root Computation

It is well known that computing square roots of matrices can be done efficiently using numerical methods (for instance, see [START_REF] Gantmacher | The Theory of Matrices[END_REF]). On the other hand, it seems difficult to control the bit complexity of numerical methods. In (Chistov et al., 1997, Section 3), the authors consider the problem of computing, in an exact way, the square root of matrices over algebraic number fields. As presented, it is not completely clear that the method proposed is polynomial-time as some coefficients of the result matrix lie in extensions of nonpolynomial size, see [START_REF] Cai | Computing Jordan Normal forms Exactly for Commuting Matrices in Polynomial Time[END_REF]. However, by applying a small trick to the proof of [START_REF] Chistov | Polynomial time algorithms for modules over finite dimensional algebras[END_REF], one can compute a solution in polynomial-time for various field of characteristic 2. In any case, for the sake of completeness, we propose two polynomial-time algorithms for this task. First, a general method fixing the issue encountered in (Chistov et al., 1997, Section 3) is presented in Section 3.2. To do so, we adapt the technique of [START_REF] Cai | Computing Jordan Normal forms Exactly for Commuting Matrices in Polynomial Time[END_REF] and compute the square root as the product of two matrices in an algebraic extension which can both be computed in polynomial time. The delicate task being to control the size of the algebraic extensions occurring during the algorithm. In here, each coefficient of the two matrices are lying in an extension field of polynomial degree in n. Furthermore, these matrices allow us to test in polynomial time if H 1 , . . . , H m and H ′ 1 , . . . , H ′ m are indeed equivalent. We then present a second simpler method based on the generalized Jordan normal form (see Section 6.3) which works (in polynomial time) over finite fields. In general, it deals with algebraic extensions of smaller degree than the first one. Putting things together, we obtain our main result:

Theorem 2. Let K be a field with char K 2. Under Assumption 1, there exists a randomized polynomial-time algorithm for solving quadratic-IP1S over an extension field of K of polynomial degree in n.

Let us note that the authentication scheme using IP1S requires to find a solution over the base field. However, it is not always necessary to find a solution in the base field (typically, in the context of a key-recovery for multivariate schemes). In [START_REF] Bettale | Cryptanalysis of HFE, Multi-HFE and Variants for Odd and Even Characteristic[END_REF], the authors recover an equivalent key over an extension for the multi-HFE scheme.

In addition, under some nondegeneracy assumption, Theorem 2 can be turned into a deterministic algorithm solving IP1S over the base field K or an extension thereof. That is: Theorem 3. Under Assumption 1 and the assumption that one of the quadratic form is nondegenerate, there is a deterministic polynomial-time algorithm for solving quadratic-IP1S over an extension of K of polynomial degree in n. Furthermore, if the space of matrices satisfying equations (4) has dimension 1, then the algorithm can solve quadratic-IP1S over K.

Let us note that assuming that one of the Hessian matrix is invertible is not a strong assumption when the size of K is not too small. Indeed, the probability of picking a random invertible symmetric matrix over F q is

∏ n i=1 (1 -q -i ) ∏ ⌊n/2⌋ i=1 (1 -q -2i ) = ⌈n/2⌉ ∏ i=1 (1 -q -2i+1 ),
see (Carlitz, 1954, Equations 4.7 and 4.8).

If m ≥ 3, for random matrices H 1 , . . . , H m , the set of solutions of equations ( 4) is a 1-dimensional matrix space. This allows us to solve quadratic-IP1S in polynomial-time over K. In Section 3.5, we present our timings for solving IP1S. These experiments confirm that for randomly chosen matrices, our method solves IP1S over K. We remark also that our method succeeds to solve IP1S over F 2 for public-keys whose sizes are much bigger than practical ones.

In Section 4, we consider the counting problem #IP1S associated to IP1S for quadratic (homogeneous) polynomials in its canonical form (as defined in Theorem 1). Note that such a counting problem is also related to cryptographic concerns. It corresponds to evaluating the number of equivalent secret-keys in MPKC, see [START_REF] Faugère | On enumeration of polynomial equivalence classes and their application to MPKC[END_REF]; [START_REF] Wolf | Equivalent keys in multivariate quadratic public key systems[END_REF]. Given homogeneous quadratic polynomials f

= ( f 1 , . . . , f m ), g = (g 1 , . . . , g m ) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . , x n ] m , we want to count the number of invertible matrices A ∈ GL n (K) such that g(x) = f(A • x).
To do so, we define:

Definition 3. Let f = ( f 1 , . . . , f m ) ∈ K[x 1 , . . . ,
x n ] m , we shall call automorphism group of f the set:

G f = {A ∈ GL n (K) | f(A • x) = f(x)}.
If f ∼ g, the automorphism groups of f and g are similar. Thus, the size of the automorphism group of f allows us to count the number of invertible matrices mapping f to g. For quadratic homogeneous polynomials, the automorphism group coincides with the subset of regular matrices in the centralizer C (H ) of the Hessian matrices H associated to f. Taking α an algebraic element of degree m over K = F q , let us assume the Jordan normal form of H = ∑ m i=i α i-1 H i has Jordan blocks of sizes s i,1 ≤ • • • ≤ s i,d i associated to eigenvalue ζ i , for i, 1 ≤ i ≤ r. Then, as a consequence of (Singla, 2010, Lemma 4.11), we prove that, if q is an odd prime power, then the number of solutions of quadratic-IP1S in F n×n q is bounded from above by:

q (∑1≤i≤r ∑ 1≤ j≤d i (2 d i -2 j+1)s i, j ) -1.
Open Question: The Irregular Case Given a quadratic instance of IP1S, a nondegenerate instance is an instance wherein the matrix whose rows are all the rows of H 1 , . . . , H m has rank n. In paragraph 2.ii, we see how to transform some degenerate instances into nondegenerate instances. However, nondegenerate instances are not always regular instances. There are cases, the so-called irregular cases, such that the vector space of matrices spanned by H 1 , . . . , H m does not contain a nondegenerate matrix. This situation is well illustrated by the following example

f 1 = x 1 x 3 , f 2 = x 2 x 3 . Any linear combination of f 1 , f 2 is degenerate, while f = ( f 1 , f 2 ) is not.
Note that we can decide in randomized polynomial time if an instance of quadratic-IP1S is irregular since it is equivalent to checking if a determinant is identically equal to zero; thus it is a particular instance of polynomial identity testing. In the irregular case, it is clear that our algorithm fails. In fact, it seems that most known algorithms dedicated to quadratic-IP1S [START_REF] Bouillaguet | Practical cryptanalysis of the identification scheme based on the isomorphism of polynomial with one secret problem[END_REF]; Macario-Rat et al. ( 2013)) will fail on these instances; making the hardness of the irregular case intriguing and then an interesting open question.

Special case of IP

In our quest of finding instances of IP solvable in polynomial-time, we take a first step in Section 5. We consider IP for a specific set of polynomials with m = n. In the aforementioned Section 5, we prove the following:

Theorem 4. Let g = (g 1 , . . . , g n ) ∈ K[x 1 , . . . , x n ] n be given in dense representation, and f = POW n,d = (x d 1 , . . . , x d n ) ∈ K[x 1 , . . . ,
x n ] n for some d > 0. Whenever char K = 0, let e = d and g = g. Otherwise, let p = char K, let e and r be integers such that d = p r e, p and e coprime, and let g ∈ K[x] n be such that g(x) = g(x p r ). Let L be a polynomial size of an arithmetic circuit to evaluate the determinant of the Jacobian matrix of g. If the size of K is at least 12 max(2 L+2 , e (n -1) 2 e(n-1) + e 3 (n -1) 3 , 2 (e (n -1) + 1) 4 ), then there is a randomized polynomial-time algorithm which recovers -if any -

(A, B) ∈ GL n (K) × GL n (K) such that: g = B • POW n,d (A • x).
This extends a similar result of (Kayal, 2011, Section 5) who considered PolyEquiv for a sum of d-power polynomials. We show that solving IP for POW n,d reduces to factoring the determinant of a Jacobian matrix (in [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF], the Hessian matrix is considered). This illustrates, how powerful partial derivatives can be in equivalence problems [START_REF] Chen | Partial derivatives in arithmetic complexity and beyond[END_REF][START_REF] Perret | A fast cryptanalysis of the isomorphism of polynomials with one secret problem[END_REF]). To go along with the proof of Theorem 4, we design Algorithm 4 at the end of Section 5.

Normalization -Canonical form of IP1S

In this section, we prove Theorem 1. In other words, we explain how to reduce, under Assumption 1, any quadratic homogeneous instance

(f, g) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . ,
x n ] m of IP1S to a suitable canonical form, i.e. an instance of IP1S where all the Hessian matrices are invertible and the first two equal the same diagonal invertible matrix. We emphasize that the reduction presented is randomized.

2.i. Homogenization. We show here that the equivalence problem over inhomogeneous polynomials with affine transformation on the variables reduces to the equivalence problem over homogeneous polynomials with linear transformation on the variables. To do so, we simply homogenize the polynomials. Let x 0 be a new variable. For any polynomial p ∈ K[x] of degree 2, we denote by

p ⋆ (x 0 , x 1 , . . . , x n ) = x 2 0 p(x 1 /x 0 , . . . , x n /x 0 ) its homogenization.
Proposition 5. IP1S with quadratic polynomials and affine transformation on the variables can be reduced in polynomial-time to IP1S with homogeneous quadratic polynomials and linear transformation on the variables.

Proof. Let (f, g) ∈ K[x] m × K[x]
m be inhomogeneous polynomials of degree 2. We consider the transformation which maps (f, g)

to f ⋆ = ( f ⋆ 0 = x 2 0 , f ⋆ 1 , . . . , f ⋆ m ), g ⋆ = (g ⋆ 0 = x 2 0 , g ⋆ 1 , . . . , g ⋆ m )
. This clearly transforms polynomials of degree 2 to homogeneous quadratic polynomials. We can write

f i (x) = x T H i x + L i x + c i with H i ∈ K n×n , L i ∈ K n and c i ∈ K, then f i (Ax + b) = (Ax + b) T H i (Ax + b) + L i (Ax + b) + c i and its homogenization is (Ax + bx 0 ) T H i (Ax + bx 0 ) + L i (Ax + bx 0 ) x 0 + c i x 2 0 = f ⋆ i (A ′ x ⋆ ), with x ⋆ = (x 0 , x 1 , . . . , x n ) T . If (A, b) ∈ GL n (K) × K n is
an affine transformation solution on the inhomogeneous instance then A ′ = 1 0 b A is a solution for the homogenized instance. Conversely, a solution A ′ ∈ GL n+1 (K) of the homogeneous problem must stabilize the homogenization variable x 0 in order to be a solution of the inhomogeneous problem. This is forced by adding f 0 = x 2 0 and g 0 = x 2 0 and setting

C ′ = A ′ /a ′ 0,0 , with a ′ 0,0 = ±1. One can see that C ′ is of the form 1 0 d C , and (C, d) ∈ GL n (K) × K n is a solution for (f, g).
2.ii. Redundant Variables. As a first preliminary natural manipulation, we first want to eliminate -if anyredundant variables from the instances considered. Thanks to [START_REF] Carlini | Reducing the number of variables of a polynomial[END_REF] (and reformulated in Kayal ( 2011)), this task can be done in randomized polynomial time:

Proposition 6. (Carlini (2005); Kayal (2011)) Let f ∈ K[x 1 . . . ,
x n ] be a polynomial. We shall say that f has s essential variables if ∃ M ∈ GL n (K) such that f (Mx) depends only on the first s variables x 1 , . . . , x s . The remaining ns variables x s+1 . . . , x n will be called redundant variables. If char K = 0 or char K > deg f , and f has s essential variables, then we can compute in randomized polynomial time M ∈ GL n (K) such that f (M x) depends only on the first s variables.

For a quadratic form, s is simply the rank of the associated Hessian matrix. As such, for m = 1, a quadratic instance is regular if and only if the associated Hessian matrix is invertible. For a set of equations, we extend the notion of essential variables as follows.

Definition 4. The number of essential variables of f = ( f 1 , . . . , f m ) ∈ K[x 1 , . . . , x n ] m is the smallest s such that f can be decomposed as:

f = f(ℓ 1 , . . . , ℓ s )
with ℓ 1 , . . . , ℓ s being linear forms in x 1 , . . . , x n of rank s and f ∈ K[y 1 , . . . , y s ] m .

The linear forms ℓ 1 , . . . , ℓ s can be easily computed thanks to Proposition 6 when the characteristic of K is zero or greater than the degrees of f 1 , . . . , f m . In characteristic 2, when K is perfect (which is always true if K is finite for instance) the linear forms can also be recovered in polynomial time (see [START_REF] Berthomieu | Computing Hironaka's invariants: Ridge and Directrix[END_REF]; [START_REF] Giraud | Étude locale des singularités. U.E.R. Mathématique[END_REF]; [START_REF] Hironaka | Additive groups associated with points of a projective space[END_REF] for instance). Below, we show that we can restrict our attention to only essential variables. Namely, solving IP1S on (f, g) reduces to solving IP1S on instances having only essential variables.

Proposition 7. Let (f, g) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . ,
x n ] m be two sets of quadratic polynomials. If f ∼ g, then their numbers of essential variables must be the same. Let s be the number of essential variables of f. Finally, let ( f, g) ∈ K[y 1 , . . . , y s ] m × K[y 1 , . . . , y s ] m be such that:

f = f(ℓ 1 , . . . , ℓ s ) and g = g(ℓ ′ 1 , . . . , ℓ ′ s ),
with ℓ 1 , . . . , ℓ s (resp. ℓ ′ 1 , . . . , ℓ ′ s ) linear forms in x of rank s and f, g ∈ K[y 1 , . . . , y s ] m . It holds that:

f ∼ g ⇐⇒ f ∼ g.
Proof. Let H 1 , . . . , H m be the Hessian matrices of f 1 , . . . , f m (resp. H ′ 1 , . . . , H ′ m be the Hessian matrices of g 1 , . . . , g m ). Similarly, we define the Hessian matrices H1 , . . . , Hm (resp. H′ 1 , . . . , H′ m ) of f1 , . . . , fm (resp. g1 , . . . , gm ). Let also M and N be matrices in GL n (K) such that H i = M T Hi 0 0 0 M and H ′ i = N T H′ i 0 0 0 N for all i, 1 ≤ i ≤ m. There exist such M and N, as f and g have essentially s variables. Up to re-indexing the rows and columns of H i and H ′ i , so that they remain symmetric, one can always choose M and N such that

M = M 1 M 2 0 Id and N = N 1 N 2 0 Id , with M 1 , N 1 ∈ GL s (K). If f ∼ g, ∃ à ∈ GL s (K) such that A T Hi à = H′ i , for all i, 1 ≤ i ≤ m. Then, for all B ∈ K (n-s)×s and C ∈ GL n-s (K): ÃT B T 0 C T Hi 0 0 0 à 0 B C = H′ i 0 0 0 , N T ÃT B T 0 C T M -T H i M -1 à 0 B C N = H ′ i .
Therefore, f and g are equivalent.

Conversely, we assume now that f ∼ g, i.e. there exists

A ∈ GL n (K) such that A T • H i • A = H ′ i , for all i, 1 ≤ i ≤ m.
This implies that:

N -T A T M T Hi 0 0 0 M A N -1 = H′ i 0 0 0 , ∀ i, 1 ≤ i ≤ m.
We then define à = ((MAN -1 ) i, j ) 1≤i, j≤s , so that f( Ãx) = g(x). As g has s essential variables, then rank à cannot be smaller than s, hence à ∈ GL s (K). We then get ÃT Hi à = H′ i for all i, 1 ≤ i ≤ m, i.e. f ∼ g.

According to Proposition 7, there is an efficient reduction mapping an instance (f, g) of IP1S to an instance ( f, g) of IP1S having only essential variables. From now on, we will then assume that we consider instances of IP1S with n essential variables for both f and g.

2.

iii. Canonical Form. We now assume that char K 2.

Definition 5. Let f = ( f 1 , . . . , f m ) ∈ K[x 1 , . . . ,
x n ] m be quadratic homogeneous forms with Hessian matrices H 1 , . . . , H m . We shall say that f is regular if its number of essential variables is n and if

∃ λ 1 , . . . , λ m ∈ K such that det (∑ m i=1 λ i H i ) 0.
Remark 8. Our algorithm requires that amongst all the Hessian matrices, one at least is invertible, the so-called regular case. It is not sufficient to only assume that the number of essential variables is n. Indeed, Ivanyos's irregular example f = (x 1 x 3 , x 2 x 3 ) has 3 essential variables, but any nonzero linear combination λ 1 f 1 + λ 2 f 2 has only 2 essential variables λ 1 x 1 + λ 2 x 2 and x 3 . Similarly, f = (x 2 1 +x 2 2 +x 2 3 , x 2 2 +2x 2 3 +x 2 4 ) has 4 essential variables but any nonzero linear combination λ 1 f 1 +λ 2 f 2 over F 3 has only 3 essential variables. This explains the additional condition on the previous definition, and our Assumption 1.

We are now in a position to reduce quadratic homogeneous instances of IP1S to a first simplified form.

Proposition 9. Let (f, g) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . ,
x n ] m be regular quadratic homogeneous polynomials. There is a randomized polynomial-time algorithm which returns "NOSOLUTION" only if f ≁ g, or a new instance

( f, g) = n ∑ i=1 d i x 2 i , f2 , . . . , fm , n ∑ i=1 d i x 2 i , g2 , . . . , gm ∈ K[x] m × K[x] m , with d 1 , . . . , d n being 1 or nonsquares in K, such that f ∼ g ⇐⇒ f ∼ g. If K is finite, the output of this algorithm is correct with probability at least 1 -n/|K|. If f ∼ g, invertible matrices P, Q and A ′ ∈ GL n (K) are returned such that f(Px) = f(x), g(Qx) = g(x) and f(A ′ x) = g(x). It then holds that f(PA ′ Q -1 x) = g(x).
Proof. Let H 1 , . . . , H m be the Hessian matrices associated to f 1 , . . . , f m . According to Schwartz-Zippel-DeMillo-Lipton Lemma [START_REF] Demillo | A probabilistic remark on algebraic program testing[END_REF]; [START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF]), we can compute in randomized polynomial time λ 1 , . . . , λ m ∈ K such that ϕ = ∑ m i=1 λ i • f i is regular, i.e. det (∑ m i=1 λ i H i ) 0. The probability to pick (λ 1 , . . . , λ m ) ∈ K m on which ϕ is regular is bounded from above by n/|K|. We define γ = ∑ m i=1 λ i • g i . Should one reorder the equations, we can assume w.lo.g. that λ 1 0. We have then:

f ∼ g ⇐⇒ (ϕ, f 2 , . . . , f m ) ∼ (γ, g 2 , . . . , g m ).
Now, applying Gauß's reduction algorithm to ϕ, there exists d 1 , . . . , d n ∈ K, each being 1 or a nonsquare, such that ϕ = ∑ n i=1 d i ℓ 2 i , where ℓ 1 , . . . , ℓ n are independent linear forms in x 1 , . . . , x n . This gives a

P ∈ GL n (L) such that f = ( φ = ∑ n i=1 d i x 2 i , f2 , . . . , f m ) = (ϕ(Px), f 2 (Px), . . . , f m (Px)). Clearly, f ∼ f, hence, f ∼ g.
After that, we can apply once again Gauß's reduction algorithm to γ. If the reduced polynomial is different from ∑ n i=1 d i x 2 i , then f ≁ g and we return "NOSOLUTION". Otherwise, the reduction is given by a matrix

Q ∈ GL n (L) such that g = ( γ = ∑ n i=1 d i x 2 i , g2 , . . . gm ) = (γ(Qx), g 2 (Qx), . . . , g m (Qx)) and g ∼ g. Thus, f ∼ g if and only if f ∼ g. Now, assume that ∃ A ′ ∈ GL n (K) such that f(A ′ x) = g(x). Then, f(PA ′ x) = g(Qx), i.e. f(PA ′ Q -1 x) = g(x).
Let us recall that whenever K = Q, computing the exact same sum of squares for f1 and g1 is difficult, see (Saxena, 2006, Chapter 3), (Wallenborn, 2013, Chapter 1). As such, we could only assume that our canonical form is g1 = ∑ n i=1 d ′ i x 2 i . This would merely change the formulation of following Theorem 10.

2.iv. Invertible Hessian Matrices.

We are now in a position to reduce any regular homogeneous quadratic instances (f, g) of IP1S to a new form of the instances where all the polynomials are themselves regular assuming we could find one. From Proposition 9, this is already the caseunder randomized reduction -for f 1 and thus g 1 . For the other polynomials, we proceed as follows. For i, 2 ≤ i ≤ m, if the Hessian matrix H i of f i is invertible, then we do nothing. Otherwise, we change H i into H i -ν i H 1 , with ν i not an eigenvalue of H i H -1 1 . As K has at least n + 1 elements, there exists such a ν i in K. This gives the following result:

Theorem 10. Let (f, g) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . ,
x n ] m be regular quadratic homogeneous polynomials. There is a randomized polynomial-time algorithm which returns "NOSOLUTION" only if f ≁ g. Otherwise, the algorithm returns two sets of n×n invertible symmetric matrices {D, H2 . . . , Hm } and {D, H′ 2 , . . . , H′ m }, with D diagonal, defined over K such that:

g(x) = f(A x), for A ∈ GL n (K) ⇐⇒ A ′-1 D -1 Hi A ′ = D -1 H′ i , ∀ i, 1 ≤ i ≤ m, for A ′ ∈ O n (K, D),
with O n (K, D) denoting the set of n × n D-orthogonal matrices over K.

Proof. Combining Proposition 9 and paragraph 2.iv any regular quadratic homogeneous instance of IP1S can be reduced in randomized polynomial time to "NOSOLUTION", only if the two systems are not equivalent, or to a

( f, g) = n ∑ i=1 d i x 2 i , f2 , . . . , fm , n ∑ i=1 d i x 2 i , g2 , . . . , gm ,
where all the polynomials are nondegenerate homogeneous quadratic polynomials in K

[x]. It follows that f ∼ g ⇐⇒ ∃ A ′ ∈ GL n (K) such that ∀ i, 1 ≤ i ≤ m, A ′ T Hi A ′ = H′ i . In particular A ′ T D A ′ = D and A ′ is D-orthogonal. Hence, A ′ T Hi A ′ = D A ′-1 D -1 Hi A ′ = H′ i , ∀ i, 1 ≤ i ≤ m.
The proof of this result implies Theorem 1.

2.v. Field Extensions and Jordan Normal Form. To amplify the success probability of our results, it will be convenient to embed a field F in some finite extension F ′ of F. This is motivated by the fact that matrices in F n×n are similar if and only if they are similar in F ′n×n , see de Seguins Pazzis (2010). In this paper, we will need to compute the Jordan normal form J of some matrix H in several situations. The computation of the Jordan normal form is done in two steps. First, we factor the characteristic polynomial, using for instance Berlekamp's algorithm over F = F q in O(n M(n) log(q n)) operations in F, where M(n) is a bound on the number of operations in F to multiply two polynomials in F[x] of degree at most n -1, see (von zur Gathen and Gerhard, 1999, Theorem 14.14). Then, we use Storjohann (1998)'s algorithm to compute the generalized eigenvectors in O(n ω ) operations in F, with ω being the exponent of time complexity of matrix multiplication, 2 ≤ ω ≤ 3.

Quadratic IP1S

In this section, we present efficient algorithms for solving regular quadratic-IP1S. According to Proposition 5, we can w.l.o.g. restrict our attention on linear changes of variables and homogeneous quadratic instances. Let D be a diagonal invertible matrix with 1 or nonsquare elements on the diagonal. Let H = {D, H 2 , . . . , H m } and H ′ = {D, H ′ 2 , . . . , H ′ m } be two families of invertible symmetric matrices in K n×n . As explained in Theorem 10, our task reduces -under a randomized process -to finding a D-orthogonal matrix A ′ ∈ O n (K, D) such that:

A ′-1 D -1 H i A ′ = D -1 H ′ i , ∀ i, 1 ≤ i ≤ m. (5)
Case D = Id was studied in (Chistov et al., 1997, Theorem 4). The authors prove that there is an orthogonal solution A, such that

H i A = A H ′ i if and only if there is an invertible matrix Y such that H i Y = Y H ′ i and H T i Y = Y H T i .
In our case, whenever D = Id, the matrices are symmetric. So, the added conditions -with the transpose -are automatically fulfilled. In [START_REF] Chistov | Polynomial time algorithms for modules over finite dimensional algebras[END_REF], the authors suggest then to use the polar decomposition of Y = AW , with W symmetric and A orthogonal. Then, A is an orthogonal solution of (5).

The main idea to compute A is to compute W as the square root of Z = Y T Y as stated in (Chistov et al., 1997, Section 3). However, in general W and A are not defined over K but over L = K( ζ 1 , . . . , ζ r ), where ζ 1 , . . . , ζ r are the eigenvalues of Z. Assuming ζ 1 is the root of an irreducible polynomial P of degree d, then ζ 2 , . . . , ζ d are also roots of the same polynomial. However, there is no reason for them to be in K[x]/(P) = K(ζ 1 ). But they will be the roots of a polynomial of degree d -1, in general, over the field K(ζ 1 ). Then, doing another extension might only add one eigenvalue in the field. Repeating this process yields a field of degree d! over K. As a consequence, in the worst case, we can have to work over an extension field of degree n!. Therefore, computing W could be the bottleneck of the method. [START_REF] Chistov | Polynomial time algorithms for modules over finite dimensional algebras[END_REF] emphasize that constructing such a square root W in polynomial time is the only serious algorithmic problem. As presented, it is not completely clear that the method proposed is efficient. They propose to compute W = √ Y T Y and then to set A = W -1 Y . According to Cai's work [START_REF] Cai | Computing Jordan Normal forms Exactly for Commuting Matrices in Polynomial Time[END_REF]), some coefficients of matrix A may lie in an extension of exponential degree. Blockwise computation (see the proof of Proposition 12) can allow us to compute such a matrix. Chistov, Ivanyos and Karpinski set y i as the restriction of Y to the ith eigenspace, associated to ζ i , of Y T Y . Then, x i = ζ i -1 y i and they return the block diagonal matrix constructed from the x i 's. However, this construction gives the impression that the ith eigenspace of Y T Y is stable by Y , as W -1 would act as a multiplication by ζ i -1 . As a consequence, the blockwise computation was not ensured.

However, this issue does not happen if one uses the same proof on W = √ Y Y T and A = Y W -1 . In the following subsection, we extend their proof to any invertible diagonal matrix D.

Existence of a D-Orthogonal Solution

The classical polar decomposition is used in (Chistov et al., 1997, Theorem 4) to determine an orthogonal solution. Using the analogous decomposition, the so-called Generalized Polar Decomposition (GPD), which depends on D, yields a D-orthogonal solution, see [START_REF] Mackey | Structured factorizations in scalar product spaces[END_REF]. The GPD of an invertible matrix Y is the factorization Y = AW , with A D-orthogonal and W in the associated Jordan algebra, i.e. W T = DW D -1 . Let us notice that A and W might be defined only over K ′ an algebraic extension of K of some degree.

Proposition 11. Let K = {K 1 , . . . , K m } and K ′ = {K ′ 1 , . . . , K ′ m } be two subsets of m matrices in K n×n . Let D be an invertible diagonal matrix. There is a D-orthogonal solution A ∈ K ′n×n to the conjugacy problem K

i A = A K ′ i for all 1 ≤ i ≤ m, if and only if there is an invertible solution Y ∈ K ′n×n to the conjugacy problem K i Y = Y K ′ i and K T i DY D -1 = DY D -1 K ′ T i for all 1 ≤ i ≤ m. Furthermore, if Y = AW is the GPD of Y with respect to D, then A suits.
Proof. This proof is a generalization of (Chistov et al., 1997, Section 3). If A is a D-orthogonal solution to the first problem, then as A T = D A -1 D -1 , it is clear that A is a solution to the second problem. Conversely, let Y be a solution to the second problem, then Z = D -1 Y T DY commutes with K ′ i . As Y is invertible, so is Z, therefore, given a determination of the square roots of the eigenvalues of Z, there is a unique matrix W with these eigenvalues such that W 2 = Z and W is in the Jordan algebra associated to D, that is W T = DW D -1 , see (Mackey et al., 2005, Theorem 6.2). As such, W is a polynomial in Z as proven in Section 6.1 and commutes with K ′ i . Finally, A = Y W -1 is an D-orthogonal solution of the first problem. As W commutes with

K ′ i , A -1 K i A = W Y -1 K i Y W -1 = W K ′ i W -1 = K ′ i and A T D A = W -T Y T DY W -1 = DW -1 D -1 Y T DY W -1 = DW -1 Z W -1 = D.
For the sake of completeness, we present several efficient algorithms for performing the square root computation.

Computing the D-Orthogonal Solution

The goal of this part is to "D-orthogonalize" an invertible solution Y ∈ GL n (K) of equation ( 5). Instead of computing exactly A ∈ O n (L, D), we compute in polynomial time two matrices whose product is A. These matrices allow us to verify in polynomial time that H i and H ′ i are equivalent for all i, 1 ≤ i ≤ m. To be more precise, we prove the following proposition.

Proposition 12. Let H = {H 1 = D, H 2 , . . . , H m } and H ′ = {H ′ 1 = D, H ′ 2 , .
. . , H ′ m } be two sets of invertible matrices in K n×n . We can compute in polynomial time two matrices representations of matrices S and T defined over an algebraic extension L such that S T -1 is D-orthogonal and for all 1

≤ i ≤ m, D -1 H i (S T -1 ) = (S T -1 ) D -1 H ′ i .
In the worst case, product S T -1 cannot be computable in polynomial time over L. However, matrices S T H i S and T T H ′ i T can be computed and tested for equality in polynomial time.

Proof. Let Y ∈ GL n (K) such that D -1 H i Y = Y D -1 H ′ i , ∀ i, 1 ≤ i ≤ m. We set Z = D -1 Y T DY .
Let us denote by T , the change of basis matrix such that J = T -1 Z T is the Jordan normal form of Z. According to [START_REF] Cai | Computing Jordan Normal forms Exactly for Commuting Matrices in Polynomial Time[END_REF], T , T -1 and J can be computed in polynomial time. Because of the issue of mixing all the eigenvalues of Z, we cannot compute efficiently A in one piece. We will then compute A T and T -1 separately. Indeed, A T (resp. T -1 ) is such that each of its columns (resp. each of its rows) is defined over an extension field K(ζ i ), where ζ 1 , . . . , ζ r are the eigenvalues of Z.

We shall say that a matrix is block-wise (resp. columnblock-wise, rowblock-wise) defined over K(ζ ) if for all 1 ≤ i ≤ r, its ith block (resp. block of columns, block of rows) is defined over K(ζ i ). The size of the ith block being the size of the ith Jordan block of J.

As J = T -1 Z T is a Jordan normal form, it is block-wise defined over K(ζ ). Using the closed formula of Section 6.1, one can compute in polynomial time a square root G of J. This matrix is a block diagonal matrix, block-wise defined over K( ζ ), hence it can be inverted in polynomial time. Should one want W , one would have to compute W = T G T -1 . Let us recall that matrices T and T -1 are respectively columnblock-wise and rowblock-wise defined over K(ζ ), see (Cai, 1994, Section 4). Since Y is defined over K, then Y T is columnblock-wise defined over K(ζ ). Thus

S = A T = Y W -1 T = Y T G -1 is columnblock-wise defined over K( ζ ).
We recall that product A T • T -1 mangles the eigenvalues and make each coefficient defined over K( ζ 1 , . . . , ζ r ) and thus must be avoided. Now, to verify that A T H A = H ′ , for any H ∈ H and the corresponding H ′ ∈ H ′ , we compute separately S T H S = T T A T H A T and T T H ′ T . For the former, S = A T (resp. S T = (A T ) T ) is columnblock-wise (resp. rowblock-wise) defined over K( ζ ) and H is defined over K. Therefore, the product matrix makes each of the coefficients which are on both the ith block of rows and the jth block of columns defined over K( ζ i , ζ j ) and so can be computed in polynomial time. For the latter, the same behaviour occurs on the resulting matrix as T is columnblock-wise defined over K(ζ ).

Let us assume that the characteristic polynomial of Z, of degree n, can be factored as P e 1 1 • • • P e s s with P i and P j coprime whenever i j, deg P i = d i and e i ≥ 1. From a computation point of view, one needs to introduce a variable α i, j for each root of P i and then a variable β i, j for the square root of α i, j . This yields a total number of 2 ∑ s i=1 d i variables. In Section 6.3, we present another method which manages to introduce only 2s variables in characteristic p > 2.

Probabilistic and Deterministic Algorithms

We first describe a simple probabilistic algorithm summarizing the method of Section 3.2.

Algorithm 1. Probabilistic algorithm.

Input Two sets of invertible symmetric matrices H

= {H 1 = D, . . . , H m } ⊆ K n×n and H ′ = {H ′ 1 = D, . . . , H ′ m } ⊆ K n×n . Output A description of the matrix A ∈ GL n (L) such that H ′ i = A T H i A for all 1 ≤ i ≤ m or "NOSOLUTION". 1. Compute the vector subspace Y = {Y | D -1 H i Y = Y D -1 H ′ i , ∀ 1 ≤ i ≤ m} ⊆ K n×n . 2.
If Y is reduced to the null matrix then return "NOSOLUTION". 3. Pick at random Y ∈ Y . 4. Compute Z = D -1 Y T DY and J = T -1 Z T ∈ L n×n , the Jordan normal form of Z together with T . 5. Compute G -1 the inverse of a square root of J. 6. Return Y T G -1 and T .

Theorem 13. Algorithm 1 is correct with probability at least 1 -n/|K| and runs in polynomial time.

Proof. The correctness and the polynomial-time complexity of the algorithm come from Section 3.2. After computing Y and putting the equations defining its matrices in triangular form, one has to pick at random one matrix Y ∈ Y . By sampling the whole field K on these free variables, the probability that detY = 0 is upper bounded by n/|K| thanks to Schwartz-Zippel-DeMillo-Lipton Lemma [START_REF] Demillo | A probabilistic remark on algebraic program testing[END_REF]; [START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF]).

Remark 14. Let us recall that the conjugacy problem does not depend on the ground field (see de Seguins Pazzis (2010)), i.e. if there exists

Y ∈ GL n (K ′ ), such that H i Y = Y H ′ i , then there exists Y ′ ∈ GL n (K) such that H i Y ′ = Y ′ H ′ i .
This allows us to extend K to a finite extension in order to decrease the probability of getting a singular matrix Y . Thus the success probability of Algorithm 1 can be amplified to 1 -n/|K ′ | for any extension K ′ ⊇ K. The probability can be then made overwhelming large by considering extension of degree O(n). In this case, the algorithm returns the description of a solution on K ′ ( ζ 1 , . . . , ζ r ). Notice also that this algorithm can be turned into a deterministic algorithm using (Chistov et al., 1997, Theorem 2). That is, there is a polynomial-time algorithm allowing to compute an invertible element in Y . Furthermore, if one of the original Hessian matrices is already invertible, the computations of the essential variables of paragraph 2.ii and the search of an equation with n essential variables in paragraph 2.iii can be done in a deterministic way. Whence, the whole algorithm is deterministic.

The main Theorem 3 summarizes this remark together with Theorem 13.

The binary Case

In this section, we investigate fields of characteristic 2. Let K = F q and (f,

g) ∈ K[x] m × K[x] m .
Instead of Hessian matrices, we consider equivalently upper triangular matrices H 1 , . . . , H m and H ′ 1 , . . . , H ′ m such that:

f i (x) = x T H i x, g i (x) = x T H ′ i x, ∀ 1 ≤ i ≤ m.
For any matrix M ∈ K n×n , let us denote ∆(M) = Diag(m 11 , . . . , m nn ) and

Σ(M) = M + M T . It is classical that if there exists A ∈ GL n (K) such that g(x) = f(A • x), then we also have Σ(H ′ i ) = A T Σ(H i ) A, (6)
∆(H ′ i ) = ∆(A T H i A), ∀ i, 1 ≤ i ≤ m. ( 7
)
It suffices for this to expand f(A • x) and to consider the upper triangular matrices. In a sense, Σ(H i ) is the Hessian matrix of f i and ∆(H i ) allow us to remember the x 2 j terms in f i . Combining two equations of (6) yields Σ(H

′ j ) -1 Σ(H ′ i ) = A -1 Σ(H j ) -1 Σ(H i )A as long as Σ(H j ) is invertible.
Let us notice that Σ(H i )'s are symmetric matrices with a zero diagonal, thus antisymmetric matrices with a zero diagonal. We would like to stress out that in odd dimension, the determinant of a symmetric matrix S with a zero diagonal is always zero. Indeed, expanding formula ∑ σ ∈S n ∏ n i=1 s i,σ (i) yields, for each nonzero term ∏ n i=1 s i,σ (i) , the term

∏ n i=1 s i,σ -1 (i) = ∏ n i=1 s σ (i),i = ∏ n i=1 s i,σ (i) .
Dimension n being odd, ∏ n i=1 s i,σ (i) cannot be the same term as ∏ n i=1 s i,σ -1 (i) . Hence they cancel each other. One can also see these matrices as projections of antisymmetric matrices over a ring of characteristic 0, namely Z q the unramified extension of the ring of dyadic integers of degree log 2 q. Let S ∈ Z n×n q be antisymmetric such that S → S. Then det S = det ST = det(-S) = (-1) n det S, hence det S = 0 and det S = 0.

Therefore, if n is odd, then a linear combination of the Σ(H i )'s will always be singular. This can be related to the irregular case of the introduction.

Reduction to canonical representations in even dimension. In this setting, we also rely on Assumption 1 to assume that a linear combination ∑ m i=1 λ i f i is not degenerate, and λ 1 , . . . , λ m can be found in randomized polynomial time. Assuming λ 1 0, we substitute the linear combinations ∑ m i=1 λ i H i and ∑ m i=1 λ i H ′ i to H 1 and H ′ 1 . As a consequence, we can find linear forms ℓ 1 , . . . , ℓ n in x such that, see (Lidl and Niederreiter, 1997, Theorem 6.30):

f 1 (x) = ℓ 1 ℓ 2 + ℓ 3 ℓ 4 + • • • + ℓ n-1 ℓ n or f 1 (x) = ℓ 1 ℓ 2 + ℓ 3 ℓ 4 + • • • + ℓ n-1 ℓ n + ℓ 2 n-1 + dℓ 2 n , where Tr K (d) = d + d 2 + • • • + d q/2 = 1.
After applying this change of variables, Σ(H 1 ) is always the following invertible block diagonal matrix:

Σ(H 1 ) = Diag 0 1 1 0 , . . . , 0 1 1 0 .
Following paragraph 2.iv, we once again choose ν i such that Σ(H i + ν i H 1 ) is invertible and replace H i by H i + ν i H 1 . Thus, Proposition 9 and Theorem 10 become:

Proposition 15. Let n be an even integer and K be a field of characteristic 2.

Let (f, g) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . ,
x n ] m be regular quadratic homogeneous polynomials. There is a randomized polynomial-time algorithm which returns "NOSOLUTION" only if f ≁ g or a new instance

( f, g) = ((δ , f2 , . . . , fm ), (δ , g2 , . . . , gm ))K[x 1 , . . . , x n ] m × K[x 1 , . . . , x n ] m
such that f ∼ g ⇐⇒ f ∼ g. Furthermore, denoting D the upper triangular matrix of f1 = g1 = δ , IP1S comes down to a Σ(D)-Orthogonal Simultaneous Matrix Conjugacy problem, i.e. conjugacy by an Σ(D)-orthogonal matrix under some constraints:

A T Σ(D) A = Σ(D) and ∀ i, 2 ≤ i ≤ m, Σ(D) -1 Σ(H ′ i ) = A -1 Σ(D) -1 Σ(H i ) A, ∆(A T H i A) = ∆(H ′ i ).
Proof. We mimic the proof of Proposition 9. We compute in randomized polynomial time λ 1 , . . . , λ m ∈ K such that ϕ = ∑ m i=1 λ i • f i is regular and we define γ = ∑ m i=1 λ i • g i . Assuming w.lo.g. λ 1 0. We have then:

f ∼ g ⇐⇒ (ϕ, f 2 , . . . , f m ) ∼ (γ, g 2 , . . . , g m ).
Computing δ the canonical quadratic form equivalent to ϕ yields a

P ∈ GL n (K) such that f = ( φ = δ , f2 , . . . , f m ) = (ϕ(Px), f 2 (Px), . . . , f m (Px)).
Then computing the canonical quadratic form equivalent to γ allows us to compare it with δ . If they are different, then f ≁ g and we return "NOSOLUTION". Otherwise, the reduction is given by a matrix Q ∈ GL n (K) such that g = ( γ = δ , g2 , . . . gm ) = (γ(Qx), g 2 (Qx), . . . , g m (Qx)). Thus, f ∼ g if and only if f ∼ g.

Finally, equations (6)

A T Σ(H i ) A = Σ(H ′ i ) for all i, 1 ≤ i ≤ m of can be rewritten as A T Σ(D) A = Σ(D) and Σ(D) -1 Σ(H ′ i ) = A -1 Σ(D) -1 Σ(H i )A for all i, 2 ≤ i ≤ m, while equations (7) ∆(A T H i A) = ∆(H ′ i ) for all i, 1 ≤ i ≤ m remain.
As a consequence, we designed the following algorithm to solve regular instances of quadratic-IP1S in even dimension over a field of characteristic 2.

Algorithm 2. Probabilistic algorithm in characteristic 2.

Input Two sets of triangular matrices

H = {H 1 = D, . . . , H m } ⊆ K n×n and H ′ = {H ′ 1 = D, . . . , H ′ m } ⊆ K n×n such that H 1 + H T 1 and H ′ 1 + H ′ T 1 are invertible.
Output A description of the matrix A ∈ GL n (L) such that

H ′ i = A T H i A for all 1 ≤ i ≤ m or "NOSOLUTION". 1. Compute the vector subspace Y = {Y | Σ(D) -1 Σ(H i )Y = Y Σ(D) -1 Σ(H ′ i ), ∀ 1 ≤ i ≤ m} ⊆ K n×n . 2. If Y is reduced to the null matrix then return "NOSOLUTION". 3. Pick at random Y ∈ Y . 4. Compute Z = Σ(D) -1 Σ(Y ) T Σ(D) Σ(Y ) and J = T -1 Z T ∈ L n×n , the Jordan normal form of Z together with T . 5. While J is not diagonal a. Pick at random Y ∈ Y . b.
Compute Z, J and T as above. 6. Compute G -1 the inverse of a square root of J. 7. Return Y T G -1 and T .

The while loop comes from the fact that unlike other characteristics, even if Z is invertible, it might not have any square roots which are polynomials in Z. In Section 6.2, we prove that there exists a square root of Z, which is a polynomial in Z if and only if Z is diagonalizable.

Open Question: The Irregular Case. As stated above, in characteristic 2, the irregular case seems to cover more cases than in other characteristics. Indeed, what is called usually a regular quadratic form in odd dimension falls in the irregular case. However, from the Hessian matrix point of view, an instance is irregular if all linear combinations of the Hessian matrices are singular over the ground field. This allows us to unify our statement about irregularity to all characteristics.

It seems to be an intriguing challenge to solve the binary case on instances with regular quadratic forms, in particular in odd dimension.

Benchmarks

We present in this section some timings of our algorithms over instances of IP1S. We created instances H = {H 1 , . . . , H m } and H ′ = {H ′ 1 , . . . , H ′ m } which are randomly alternatively equivalent over F p , equivalent over F p 2 but not F p or not equivalent at all over Fp , the algebraic closure of F p , for an odd p. We report our timings in the following Table 1 obtained using one core of an INTEL CORE I7 at 2.6GHz running MAGMA 2.19, [START_REF] Bosma | The Magma algebra system. I. The user language[END_REF], on LINUX with 16GB of RAM. These timings corresponds to solving the linear system which is the dominant part in our algorithm with complexity O(n 2 ω ). The code is accessible on the first author's webpage http://www-polsys.lip6.fr/~berthomieu/IP1S.html. To simplify the presentation, we only considered the case when m = n. That is, we only considered n matrices of size n.

Since our matrices are randomly chosen, we apply the following strategy. We first solve the linear system H -1

1 H i A = A H ′-1 1 H ′ i , for all i, 2 ≤ i ≤ m.
In fact, in practice, i = 2, 3 give enough equations to retrieve A up to one free parameter if H and H ′ are indeed equivalent. If the matrices are not equivalent, this linear system will return the zero matrix only.

Then, to determine A, we solve one quadratic equation amongst the ones given by A T H 1 A = H ′ 1 . Let us notice that either all these equations can be solved over F p or none of them can. If they can, then H and H ′ are equivalent over F p and we have determined A up to a sign, otherwise H and H ′ are only equivalent over F p 2 but not F p and we also have computed such an A. This yields Algorithm 3. Algorithm 3. Simplified Algorithm.

Input Two sets of generic invertible symmetric matrices

H = {H 1 , . . . , H m } ⊆ K n×n and H ′ = {H ′ 1 , . . . , H ′ m } ⊆ K n×n . Output A matrix A ∈ GL n (K) such that H ′ i = A T H i A for all 1 ≤ i ≤ m or "NOSOLUTION". 1. Compute the vector subspace Y = {Y | H -1 1 H i Y = Y H ′ 1 -1 H ′ i , ∀ 2 ≤ i ≤ 3} ⊆ K n×n . 2. If Y is reduced to a space of singular matrices then return "NOSOLUTION". 3. Determine Y 0 such that Y = {λY 0 | λ ∈ K}. 4. Solve in λ one equation λ 2 (Y T 0 H 1 Y 0 ) i, j = (H ′ 1 ) i, j for a suitable pair (i, j). 5. Set A = λ Y 0 . 6. Pick at random r ∈ K n . 7. Check that A T H i A r = H ′ i r for all i, 1 ≤ i ≤ m. 8. Return A.
Complexity estimate. Taking the first 3 matrix equations

H -1 1 H i Y = Y H ′ 1 -1 H ′ i ,
in the n 2 unknowns, one can solve this system in O(n 2 ω ) operations in K. Then, one needs to determine λ by extracting one square root in K which can be done in O((log q) 3 ) operations in K = F q with Tonelli-Shanks's algorithm, [START_REF] Shanks | Five number-theoretic algorithms[END_REF]. Finally, one can check that A T H i A = H ′ i for all i, with high probability, by picking up at random a vector and checking that the products of this vector with both sets of matrices coincides. This can be done in O(m n 2 ) operations in K.

Recall that, in practice, the best matrix multiplication algorithm is due to [START_REF] Strassen | Gaussian elimination is not optimal[END_REF] whose complexity is in O(n log 2 7 ) ⊆ O(n 2.807 ). Thus, our complexity is in O(n 5.615 ). This complexity is well confirmed since multiplying by 2 the sizes and the number of matrices multiplies our timings roughly by at most 50.

In Table 2, we report our timings for solving the linear system of our algorithm in characteristic 2 presented in Section 3.4. Our method does not differ much from the one in odd characteristic. We pick at random two sets of m upper triangular matrices over F 2 which are either equivalent over F 2 or not equivalent at all over F2 , the algebraic closure of F 2 . We first solve the linear system

Σ(H 1 ) -1 Σ(H i ) A = A Σ(H ′ 1 ) -1 Σ(H ′ i ), for all i, 2 ≤ i ≤ m.
Let us notice that in dimension 2, the linear system does not yield any information on A. In dimensions 8 or more (resp. 4 and 6), if H and H ′ are equivalent, the linear system yields in general A up to one free parameter if m ≥ 3 (resp. m ≥ 5). Otherwise, it yields the zero matrix. Then, it suffices to solve one of the quadratic equations amongst the one given by

A T Σ(H 1 ) A = Σ(H ′ 1 ) and ∆(A T H i A) = ∆(H ′ i ), for all i, 1 ≤ i ≤ m (see Proposition 15).
We compare the timings of both MAGMA and the C library M4RI, due to [START_REF] Albrecht | The M4RI Library -Version 20121224[END_REF].

Once again, our complexity in O(n 2 ω ) is well confirmed by our timings. Thanks to the linear system which totally determines A up to one free parameter, we just need to set this parameter to 1 to obtain A. This also explains why our timings are better than over F 65521 although it would seem a lot of quadratic equations must be solved.

Counting the Solutions: #IP1S

In this part, we present a method for counting the number of solutions to quadratic-IP1S. The main result is a consequence of (Singla, 2010, Lemma 4.11). According to Proposition 5, this is equivalent to enumerating all the invertible linear transformations on the variables between two sets of quadratic homogeneous polynomials. We provide here an upper bound on the number of solutions. We consider in this part regular quadratic homogeneous instances (f, g)

∈ K[x] m × K[x] m . Let H = {H 1 = D, . . . , H m } and H ′ = {H ′ 1 = D, .
. . , H ′ m } be the Hessian matrices in K n×n of f and g respectively. Our counting problem is equivalent to enumerating the number of D-orthogonal matrices X satisfying:

X -1 D -1 H i X = D -1 H ′ i , ∀ i, 1 ≤ i ≤ m. (8)
In (Singla, 2010, Section 4), the author computes the set of all matrices commuting with a given matrix. In particular, from Lemma 4.11 of the aforementioned paper, we can determine the size of this set and thus our upper bound on the number of solutions to quadratic-IP1S. In order to be self-contained, the proofs of the following lemmas shall be found in Appendix A.

Let us notice that if X and X ′ are both orthogonal solutions of (8), then X X ′-1 commutes with D -1 H (resp. X -1 X ′ commutes with D -1 H ′ ). Therefore, the size of the set of solutions is upper bounded by the number of invertible elements in the centralizer C (D -1 H ) of D -1 H .

Let α be an algebraic element of degree m over K and let K ′ = K(α). We consider the matrix

H = D -1 H 1 + • • • + α m-1 H m ∈ K ′n×n . It is clear that a matrix X ∈ K n×n is such that X -1 D -1 H i X = D -1 H i for all i, 1 ≤ i ≤ m if
and only if X -1 HX = H. Hence, the problem again reduces itself to the computation of the centralizer C (H) of H intersected with GL n (K). To ease the analysis, we consider the subspace V = C (H) ∩ K n×n of matrices in K n×n commuting with H. This provides an upper bound on the number of solutions. The dimension of V as a K-vector space is upper bounded by the dimension of C (H) as a K ′ -vector space. Indeed,

V ⊗ K ′ ⊆ C (H), hence dim K V = dim K ′ (V ⊗ K ′ ) ≤ dim K ′ C (H).
Since we only want the size of the centralizer of H, we can restrict our attention to the centralizer of the Jordan normal form J of H defined over a field L.

Let us denote ζ 1 , . . . , ζ r the eigenvalues of J. According to (Singla, 2010, Lemma 4.11) and Lemma 23 in Appendix A, if J is made of Jordan blocks of size

s i,1 ≤ • • • ≤ s i,d i for i, 1 ≤ i ≤ r, then centralizer of H has dimension at most ∑ 1≤i≤r ∑ 1≤ j≤d i (2d i -2 j + 1)s i, j .
As a consequence, if q is an odd prime power, then the following corollary gives an upper bound on the number of solutions of quadratic-IP1S in F n×n q .

Corollary 16. Let H 1 , . . . , H m ∈ F n×n q be symmetric matrices. Let α be algebraic over F q of degree m. Let H = ∑ m i=1 α i-1 H i ∈ F n×n q m and let J be its normal Jordan form with eigenvalues ζ 1 , . . . , ζ r . Assuming the blocks of J associated to ζ i are J ζ i ,s i,1 , . . . , J ζ i ,s i,d 1 with s i,1 ≤ • • • ≤ s i,d i for i, 1 ≤ i ≤ r, then the number of solutions of quadratic-IP1S in F n×n q on the instance (H 1 , . . . , H m ) is at most q (∑1≤i≤r ∑ 1≤ j≤d i (2 d i -2 j+1) s i, j ) -1.

As mentioned in the introduction, the counting problem considered here is related to cryptographic concerns. It corresponds to evaluating the number of equivalent secret keys in MPKC (see [START_REF] Faugère | On enumeration of polynomial equivalence classes and their application to MPKC[END_REF]; [START_REF] Wolf | Equivalent keys in multivariate quadratic public key systems[END_REF]). In particular, in [START_REF] Faugère | On enumeration of polynomial equivalence classes and their application to MPKC[END_REF], the authors propose an "ad-hoc" method for solving a particular instance of #IP1S. An interesting open question would be to revisit the results from [START_REF] Faugère | On enumeration of polynomial equivalence classes and their application to MPKC[END_REF] with our approach.

Special Case of the general IP Problem

In this part, we present a randomized polynomial-time algorithm for the following task:

Input: g = (g 1 , . . . , g n ) ∈ K[x 1 , . . . , x n ] n , and POW n,d = (x d 1 , . . . , x d n ) ∈ K[x 1 , . . . , x n ] n for some d > 0. Question: Find -if any -(A, B) ∈ GL n (K) × GL n (K) such that: g = B • POW n,d (A • x), with x = (x 1 , . . . , x n ) T .
In [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF], the author proposes a randomized polynomial-time algorithm for solving the problem when B is the identity matrix and m = 1 of finding A such that g(x) = f (A • x) with f (x) = ∑ n i=1 x d i . We generalize this result to m = n with an additional transformation on the polynomials. The main tool of our method is the following theorem.

Theorem 17. Let g = (g 1 , . . . , g n ) be polynomials of degree d over K[x 1 , . . . , x n ] given in dense representation. Let L be a polynomial size of an arithmetic circuit to evaluate the determinant of the Jacobian matrix of g. If the size of K is at least 12 max(2 L+2 , d (n -1) 2 d (n-1) + d 3 (n -1) 3 , 2 (d (n -1) + 1) 4 ), then one can factor the determinant of the Jacobian matrix of g in randomized polynomial time.

Proof. For this, we will use [START_REF] Kaltofen | Factorization of polynomials given by straight-line programs[END_REF]'s algorithm to factor a polynomial given by evaluation, the needed size of K is a consequence of this. As g has at most n n+d-1 d ∈ O(n d+1 ) monomials, it can be evaluated in polynomial time using a multivariate Horner's scheme. Each ∂ g i ∂ x j (a) is obtained as the coefficient in front of x j of the expansion of g i (a 1 , . . . , a j-1 , a j + x j , a j+1 , . . . , a n ) which is a univariate polynomial of degree at most d. By [START_REF] Bini | Polynomial and Matrix Computations[END_REF], Chapter 1, Section 8), this can be computed as the shift of a polynomial in polynomial time. Hence the Jacobian matrix of g at a can be evaluated in polynomial time with an arithmetic circuit of polynomial size L. This circuit can be for instance the evaluation of the Jacobian matrix of g followed by a Gaussian elimination on the matrix to compute the determinant. The determinant of the matrix can be recovered by linear algebra in O(n ω ) operations, with ω being the exponent of time complexity of matrix multiplication, 2 ≤ ω ≤ 3. Using the arithmetic circuit of polynomial size L to evaluate the determinant of the Jacobian matrix, one can use KALTOFEN's algorithm to factor it in polynomial time.

As in [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF] or [START_REF] Perret | A fast cryptanalysis of the isomorphism of polynomials with one secret problem[END_REF], we use partial derivatives to extract matrices A and B. The idea is to factor the Jacobian matrix (whereas [START_REF] Kayal | Efficient algorithms for some special cases of the polynomial equivalence problem[END_REF] uses the Hessian matrix) of g at x which is defined as follows:

J g (x) = ∂ j g i = ∂ g i ∂ x j 1≤i≤m 1≤ j≤n
.

According to the following lemma, the Jacobian matrix is especially useful in our context:

Lemma 18. Let f = ( f 1 , . . . , f m ), g = (g 1 , . . . , g m ) ∈ K[x 1 , . . . , x n ] m × K[x 1 , . . . , x n ] m . If (A, B) ∈ GL n (K) × GL m (K) are such that g = B • f(A • x), then J g (x) = B • J f (A • x) • A. As a consequence, det J g (x) = det A • det B • det J f (A • x).
As long as char K does not divide d, the Jacobian matrix of f = POW n,d (x) is an invertible diagonal matrix whose diagonal elements are J

f (x) i,i = d • x d-1 i , ∀ i, 1 ≤ i ≤ n. Thus: det J POW n,d (x) = d n n ∏ i=1 x d-1 i .

This gives

Lemma 19. Let g = (g 1 , . . . , g n ) ∈ K[x 1 , . . . , x n ] n . Let d > 0 be an integer, and define POW n,d = (x d 1 , . . . ,

x d n ) ∈ K[x 1 , . . . , x n ] m . If (A, B) ∈ GL n (K) × GL n (K) are such that g(x) = B • POW n,d (A • x), then: det J g (x) = c • n ∏ i=1 ℓ i (x) d-1 ,
with c ∈ K \ {0}, and the ℓ i 's are linear forms whose coefficients are the ith rows of A.

Proof. According to Lemma 18, det J g

(x) = det(A) • det(B) • d n • ∏ n i=1 ℓ i (x) d-1
. From Lemmas 19, we can derive a randomized polynomial-time algorithm for solving IP on the instance

(f = POW n,d , g) ∈ K[x 1 , . . . , x n ] n × K[x 1 , . . . ,
x n ] n in characteristic 0. It suffices to use [START_REF] Kaltofen | Factorization of polynomials given by straight-line programs[END_REF]'s algorithm for factoring det J g (x) in randomized polynomial time.

This allows us to recover -if any -the change of variables A. The matrix B can be then recovered by linear algebra, i.e. solving a linear system of equations. This proves the result announced in the introduction for IP, that is Theorem 4 whenever char K ∤ d. with A (p r ) = a p r i, j 1≤i, j≤n and x p r = x p r 1 , . . . , x p r n . Thus g is a polynomial in x p r and by replacing

x p r by x, the problem comes down to checking if g = B • POW n,e (A (p r ) • x) where g(x) = g(x p r ). Now, as g(x) = B • POW n,e A (p r ) x , then

J g = B • J POW n,e (A (p r ) x) • A (p r ) . Hence, det J g(x) = det B det J POW n,e (A (p r ) x) det A (p r ) = e n ∏ n i=1 li (x)
e-1 (det A) p r det B, where the li 's are linear forms whose coefficients are the ith rows of A (p r ) . Then, to use KALTOFEN's algorithm, one must set a low enough probability of failure ε yielding a big enough set of sampling points, see [START_REF] Kaltofen | Factorization of polynomials given by straight-line programs[END_REF], Section 6, Algorithm, Step R). In particular, if the arithmetic circuit for evaluating the determinant of the Jacobian we want to factor has size L, then the size of the sampling set must be greater than 6 ε max 2 L+2 , e(n -1) 2 e(n-1) + e 3 (n -1) 3 , 2 (e(n -1) + 1) 4 , recalling that our polynomial has degree e(n -1). In other words, if the probability of failure is less than 1/2, then one must consider a field of size at least 12 max 2 L+2 , e(n -1) 2 e(n-1) + e 3 (n -1) 3 , 2 (e(n -1) + 1) 4 .

All in all, this allows us to retrieve -if any -the change of variables A (p r ) and thus A. Then B can be recovered by linear algebra. This proves Theorem 4 for any characteristic, as in the introduction.

Theorem 4 (restated). Let g = (g 1 , . . . , g n ) ∈ K[x 1 , . . . , x n ] n be given in dense representation, and

f = POW n,d = (x d 1 , . . . , x d n ) ∈ K[x 1 , . . . ,
x n ] n for some d > 0. Whenever char K = 0, let e = d and g = g. Otherwise, let p = char K, let e and r be integers such that d = p r e, p and e coprime, and let g ∈ K[x] n be such that g(x) = g(x p r ). Let L be a polynomial size of an arithmetic circuit to evaluate the determinant of the Jacobian matrix of g. If the size of K is at least 12 max(2 L+2 , e(n -1) 2 e(n-1) + e 3 (n -1) 3 , 2 (e(n -1) + 1) 4 ), then there is a randomize polynomial-time algorithm which recovers -if any -(A, B) ∈ GL n (K) × GL n (K) such that:

g = B • POW n,d (A • x).
The computation of such a pair (A, B) is summarized in the following Algorithm. b. Compute g such that g(x p r ) = g(x). 2. Else e = d and g = g. 3. Create an arithmetic circuit of polynomial size L to evaluate det J g (x). 4. Evaluate det J g (x) in at least 12 max(2 L+2 , e(n -1) 2 e(n-1) + e 3 (n -1) 3 , 2 (e(n -1) + 1) 4 ) distinct points. 5. Factor det J g (x) with KALTOFEN's algorithm.

If the factorization is

c ∏ n i=1 ℓ i (x) e-1 then A = ℓ e/d i, j 1≤i, j≤n .
7. Else return "NOSOLUTION". 8. Compute B such that g(x) = B • (ℓ 1 (x) e , . . . , ℓ n (x) e ) T . 9. Return (A, B).

Square Root of a Matrix

In this section, we present further algorithms for computing the square root of a matrix. We use the same notation as in Section 3. A square root of a matrix Z is a matrix whose square is Z. In the first subsection, we deal with some properties of the square root of a matrix in characteristic not 2. In particular, we show that an invertible matrix Z always has a square root which is a polynomial in Z. In the second subsection, we consider the case of characteristic 2. We recall that whenever Z is not diagonalizable, then Z might have a square root but it is never a polynomial in Z. We give some examples of such matrices Z. Lastly, we propose an alternative to the method of Section 3 for computing the square root of a matrix in polynomial time for any field of characteristic p 2.

6.1. The square root as a polynomial in characteristic not 2

In this part, we prove that an invertible matrix always has a square root which is a polynomial in considered matrix. More specifically, we shall prove the following result.

Proposition 20. Let Z ∈ K n×n be an invertible matrix whose eigenvalues are ζ 1 , . . . , ζ r . Let ω 1 , . . . , ω r be such that ω 2 i = ζ i and ζ i = ζ j ⇒ ω i = ω j , for all 1 ≤ i, j ≤ r. Then, there exists W ∈ K(ω 1 , . . . , ω r )[Z] a square root of Z whose eigenvalues are ω 1 , . . . , ω r .

Proof. Let T be a matrix of change of basis, such that J

= T -1 Z T is made of Jordan blocks. It is clear that W such that W 2 = Z is a polynomial in Z, i.e. W = Q(Z), if and only if G = T -1 W T satisfies G = Q(J). Let J ζ ,d
be the Jordan block of size d associated with eigenvalue ζ and ω be a square root of ζ . We shall first prove that the square root G

ω,d of J ζ ,d is a polynomial in J ζ ,d with coefficients in K(ω). Matrix J ζ ,d -ζ Id d is nilpotent of degree d.
G ω,d = ω d-1 ∑ k=0 1/2 k ζ -k J ζ ,d -ζ Id k = d-1 ∑ k=0 1/2 k ω 1-2 k J ζ ,d -ζ Id k = Q ζ (J ζ ,d ) (9) =       ω 1/2 1 ω -1 • • • 1/2 d-1 ω 3-2 d . . . . . . . . . . . . 1/2 1 ω -1 ω       , with Q ζ (x) = ∑ d-1 k=0 1/2 k ω 1-2 k (x -ζ ) i ∈ K(ω)[x]
. It remains to prove that for multiple Jordan blocks, one can find a common polynomial. From equation ( 9), we deduce that G is a polynomial in J = Diag J ζ 1 ,d 1 , . . . , J ζ r ,d r if and only if there exists a polynomial Q such that Q = Q ζ i mod (X -ζ i ) d i , for all i, 1 ≤ i ≤ r. By the Chinese Remainder Theorem, this can always be solved as soon as

ζ i = ζ j implies Q ζ i = Q ζ j mod (X - ζ i ) min(d i ,d j )
, which is exactly the condition ω i = ω j .

Let us notice that picking the same square root for two equal eigenvalues is necessary. Indeed, although W = 1 0 0 -1 is a square root of Z = 1 0 0 1 , W K[Z].

Matrices with square roots in characteristic 2

In this part, we consider the trickier case of computing the square root of a matrix over a field K with char K = 2. Unfortunately, unlike other characteristics, an invertible matrix has not necessarily a square root over K. In fact, no Jordan block of size at least 2 has any square root. This is mainly coming from the fact that generalized binomial coefficients 1/2 k , involved in the Taylor expansion, are meaningless in characteristic 2.

Proposition 21. Let Z ∈ K n×n be a Jordan normal form with blocks J 1 , . . . , J r of sizes d 1 , . . . , d r ≥ 2, associated to eigenvalues ζ 1 , . . . , ζ r and blocks of sizes 1 with eigenvalues υ 1 , . . . , υ s . We assume that J 1 , . . . , J r are ordered by decreasing sizes and then eigenvalues. Matrix Z has a square root W if and only if d 1 -d 2 ≤ 1 and ζ 1 = ζ 2 , d 3 -d 4 ≤ 1 and ζ 3 = ζ 4 , etc. and if for each J i of size 2 that is not paired with J i-1 or J i+1 , then there exists a j such that υ j = ζ i .

Furthermore, matrix W is a polynomial in Z if and only if Z is diagonalizable.

Before, proving this result, we give some example of matrices with or without square roots. Following matrices J and J ′ both have two Jordan blocks associated with eigenvalue ζ . Denoting ω the square root of ζ , then K is the square root of J and K ′ 1 , K ′ 2 are those of J ′ for x, y, z any. As one can see, none of K, K ′ 1 and K ′ 2 are polynomials in J or J ′ because of the nonzero subdiagonal elements 1/x and 1/y. Examples of matrices without square roots are J ′′ , with two Jordan blocks associated with ζ of sizes 1 and 3, and J ′′′ , with three Jordan blocks associated with ζ of size 2. Computing a square root of each of them yields an inconsistent system. and one can deduce that ζ 2 is the sole eigenvalue of J 2 but that its geometric multiplicity is 2. Hence the Jordan normal form of J 2 is made of two Jordan blocks. As (J -ζ Id) d = 0 and (J -ζ Id) e 0 for all e < d, then J 2 -ζ 2 Id By this result, if Z is a square, then one must be able to pair up its Jordan blocks with same eigenvalue ζ so that the sizes differ by at most 1. The blocks that need not be paired being the blocks of size 1.

J ′′ =     ζ 0 0 0 0 ζ 1 0 0 0 ζ 1 0 0 0 ζ     , J ′′′ = Diag ζ 1 0 ζ , ζ 1 0 ζ , ζ 1
Conversely, assuming one can pair up the Jordan blocks of Z with same eigenvalue ζ so that the sizes differ by at most 1 and the remaining blocks have sizes 1. Then, each pair of blocks is the Jordan normal form of the square of a Jordan block of size the sum of the sizes and eigenvalue ζ . Furthermore, each lonely block of size 1 associated with ζ is the square of the block of size 1 associated with ζ .

Finally, for the last statement, the if part is easy. It remains the only if part for which we assume W 2 = Z and Z is not diagonalizable. Let J be the Jordan normal form of Z with blocks J 1 , . . . , J r . For any polynomial P, P(J) is also block diagonal with blocks P(J 1 ), . . . , P(J r ). Thus, if P(J) 2 = J, then P(J i ) 2 = J i for all 1 ≤ i ≤ r, which is false, unless J i has size 1.

Computation in characteristic p ≥ 2

In this part, we present an alternative method to the one presented in Section 3.2. We aim at diminishing the number of variables needed in the expression of the square root. However, this method does not work in characteristic 0. For the time being, we consider char K > 2. However, we shall see below how to adapt this method to the characteristic 2.

The idea is still to perform a change of basis T over K so that J = T -1 Z T has an easily computable square root. This matrix J is the generalized Jordan normal form, also known as the primary rational canonical form of Z. As the classical Jordan normal form, if Z is diagonalizable over K, then J is block diagonal, otherwise it is a block upper triangular matrix. Its diagonal blocks are companion matrices C (P 1 ), . . . , C (P r ) of irreducible factors P 1 , . . . , P r of its characteristic polynomial. Superdiagonal blocks are zero matrices with eventually a 1 on the bottom-left corner, if the geometric multiplicity associated to the roots of one the P i is not large enough. In other words, it gathers d conjugated eigenvalues in one block of size d which is the companion matrix of their shared minimal polynomial. Let us note that computing such a normal form can be done in polynomial time and that the change of basis matrix T is defined over K, see [START_REF] Matthews | A rational canonical form algorithm[END_REF]; [START_REF] Storjohann | An O(n 3 ) algorithm for the Frobenius normal form[END_REF]. Thus, after computing a square root G of J, one can retrieve W and A of Section 3.2 in O(n ω ) operations in the field of coefficients of G, with ω being the exponent of the time-complexity of matrix multiplication 2 ≤ ω ≤ 3. Furthermore, computing a square root of J is equivalent to computing the square root of each companion matrix. Finally, using the same argument as for the more classical Jordan normal form in Section 6.1, G is a polynomial in J. In the following, we only show how to determine the square root of a companion matrix C (P), for an irreducible P.

  Small characteristic. If char K divides d, we must change a little bit our strategy. Let us write d = p r e with char K = p and e coprime. Then, POW n,d (A x) d (A x) = POW n,e A (p r ) x p r ,

Algorithm 4 .

 4 IP for POW n,d and g.Input One set of polynomials g = (g 1 , . . . , g n ) ∈ K[x 1 , . . . , x n ] n , homogeneous of degree d. Output Two matrices A, B ∈ GL n (K) -if any -such that g(x) = B•POW n,d (A•x) with POW n,d (x) = (x d 1 , . . . , x d n ) or "NOSOLUTION".1. If char K = p > 0 then a. Compute r, e such that d = p r e with p and e coprime.

  Hence, by the classical Taylor expansion of the square root near Id d , one can write

  Let J be a Jordan block of size d associated to eigenvalue ζ . Then J 2 -ζ 2 Id = 0 Id d-2 0 0

Table 1 :

 1 Timings for solving IP1S over F 65521 in s.

	n	20	30	40	50	60	70	80	90	100
	Timings (MAGMA) 0.010 0.030 0.080 0.25	0.68 1.4	3.2	6.25 16
	Timings (M4RI)			0.010 0.030 0.06 0.14 0.27 0.51	0.91

Table 2 :

 2 Timings for solving IP1S over F 2 in s.

  ⌈d/2⌉ = 0 and J 2 -ζ 2 Id e 0 for e < ⌈d/2⌉, i.e. e < d/2 if d is even and e < (d + 1)/2 if d is odd. Thus the Jordan normal form of J 2 has a block of size exactly ⌈d/2⌉. That is, if d is even, both blocks have size d/2 and if d is odd, one block has size (d + 1)/2 and the other block has size (d -1)/2.

We would like to thank G. Ivanyos for having pointed us this issue in a preliminary version of this paper.

Acknowledgements

We would like to thank Gabor IVANYOS for his helpful remarks and references on the irregular case. We wish to thank Gilles MACARIO-RAT for the many discussions about isomorphism of quadratic polynomials and Nitin SAXENA for those about graph isomorphism.

We thank the anonymous referees for their careful reading and their helpful comments. This work has been partly supported by the French National Research Agency ANR-11-BS02-0013 HPAC project.

Let P = x d + p d-1 x d-1 + • • • + p 0 , let us recall that the companion matrix of P is

.

If polynomial P can be decomposed as P(z) = (z -α 0 ) • • • (z -α d-1 ), then we want to find a polynomial Q such that

, where β 2 i = α i for all 0 ≤ i ≤ d -1. Let us notice that

As a consequence, the characteristic polynomial of

But since P is irreducible over K, by the invariant factors theory, then C (Q) 2 must be similar to the companion matrix C (P) over K.

As P is irreducible over K = F q , up to reindexing the roots of P, the conjugates α 1 , . . . , α d-1 of α 0 are just its iterated qth powers. Denoting L = K[x]/(P(x)) = F q d , let us assume that S(y) = y 2 -x is reducible in L[y], then β 0 ∈ L. As such, one can choose β i = β q i 0 , the iterated qth powers. In that case, the previous equations can be rewritten

As a consequence, Q(z) ∈ K[z] and to compute Q(z), we need to compute y q i effectively. This is done by computing the following values in O(d log q) operations in L:

) and we know that the resulting polynomial is in K

Whenever α 0 is not a square in L, that is whenever S(y) is irreducible, then β q d 0 is a square root of α 0 different from β 0 , thus it is -β 0 . As a consequence, setting

) would yield a polynomial that is not stable by the Frobenius endomorphism. As such, we introduce a new variable y to represent the field L ′ = L[y]/(y 2 -x) and to compute Q(z), we need to compute y q i(d+1) effectively. Since y q i = y y q i -1 = y x q i -1 2 , we can compute the following values in O(d log q) field operations in L:

As a first step, we compute in d operations, the dehomogenized polynomial in y,

But, denoting by i 0 = i mod 2, we have

2 . Hence we deduce:

Complexity analysis. Since the number of operations for computing the square root of a block of size d is bounded by O(d log q) operations in L = F q d , this is also bounded by O(d M(d) log q) operations in K = F q , where M(n) is a bound on the number of operations in K to multiply two polynomials in K[x] of degree at most n -1. As a consequence, the computation of W can be done in no more than O(n ω + nM(n) log q) operations in K. Let us assume that the characteristic polynomial of Z has degree n and can be factored as P e 1 1 • • • P e s s with P i and P j coprime whenever i j, deg P i = d i and e i ≥ 1. From a computation point of view, in the worst case, one needs to introduce a variable α i for one root of P i and a variable β i for the square root of α i , assuming α i is not a square. This yields a total number of 2s variables.

Computation in characteristic 2. The case of characteristic 2 is almost the same. From a polynomial

i , for all 1 ≤ i ≤ d -1. Thus, Q can be computed in O(d log q) operations in K and as a consequence, W in O(n ω + n log q) operations in K.

However, let us recall that D is block diagonal if and only if the Jordan normal form is block diagonal. As such, a square root of D is a polynomial in D if and only if D is block diagonal, see Section 6.2.

Appendix A. Proofs of #IP1S

In this appendix, we shall prove the dimension of the centralizer of a matrix J, a Jordan normal form. This dimension, a consequence of (Singla, 2010, Lemma 4.11), is used in Section 4 to determine an upper bound on the counting problem of quadratic-IP1S. As stated by Singla, the proofs only involve matrix multiplications are given in order for the paper to be self-contained.

First, let us recall that the centralizer of J, a Jordan block of size s is the set of upper triangular Toeplitz matrices of size s × s. Indeed, if X commutes with J, X J -JX is as such

This small result is used in the following Lemma to determine the centralizer of a Jordan normal form.

Lemma 22. Let J be a Jordan normal form. For 1 ≤ i ≤ r, let us denote J i the ith block of J and let us assume it is associated with eigenvalue ζ i and it is of size s i . Let X = (X i, j ) 1≤i, j≤r be a block-matrix, with X i, j ∈ L (ζ 1 , . . . , ζ r ) s i ×s j , that commutes with J. If ζ i = ζ j , then X i, j is an upper triangular Toeplitz matrix whose nonnecessary zero coefficients are the one on the first min(s i , s j ) diagonals. Otherwise, X i, j = 0.

Proof. We assume that r

s 1 and X 2,2 with J 2 = J ζ 2 ,s 2 . Thus they are upper triangular Toeplitz matrices.

From X 2,1 J 2 -J 1 X 2,2 , one deduces that (ζ 1 -ζ 2 ) x s 1 +s 2 ,1 = 0, hence either ζ 1 = ζ 2 or x s 1 +s 2 ,1 = 0. If ζ 1 ζ 2 , then step by step, one has X 1,2 = 0. Assuming ζ 1 = ζ 2 , then step by step, one has x s 1 +i,1 = 0 for i > 1 and since x s 1 +i+1, j+1 -x s 1 +i, j = 0 for all i, j, one has in fact that X 1,2 is a upper 36 triangular Toeplitz matrix with potential nonzero coefficients on the first min(s 1 , s 2 ) diagonals. The same argument applies to X 2,1 .

The case r > 2 is an easy generalization of this result.

From this lemma, we can deduce easily the dimension of the centralizer of a matrix.

Lemma 23. Let H ∈ K n×n be a matrix and let J be its normal Jordan form. Assuming the blocks of J associated to

Proof. Let L be the smallest field over which J is defined. It is clear that the centralizer of H over L, denoted W , contains C (H) ⊗ L. Hence, dim K C (H) = dim L (C (H) ⊗ L) ≤ dim L W . Now, let X = (X i, j ) 1≤i, j≤d 1 +•••+d r ∈ V . From Lemma 22, there are ∑ 1≤i≤r ∑ 1≤ j≤d i s i, j free parameters for the diagonal blocks of X and 2 ∑ 1≤i≤r ∑ 1≤ j<k≤d i min(s i, j , s i,k ) = 2 ∑ 1≤i≤r ∑ 1≤ j≤d i (d ij)s i, j free parameters for the off-diagonal blocks of X . This concludes the proof.

As a consequence, the number of invertible matrices in C (H) is bounded from above by q (∑1≤i≤r ∑ 1≤ j≤d i (2 d i -2 j+1) s i, j ) -1, as stated in Corollary 16.