---
abstract: |
    Let $\mathbf{f}=(f_1,\ldots,f_m)$ and $\mathbf{g}=(g_1,\ldots,g_m)$ be two sets of $m\geq 1$ nonlinear polynomials in $\mathbb{K}[x_1,\ldots,x_n]$ ($\mathbb{K}$ being a field). We consider the computational problem of finding – if any – an invertible transformation on the variables mapping $\mathbf{f}$ to $\mathbf{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 (<span style="font-variant:small-caps;">ISSAC</span> 1997) to be equivalent to finding an invertible matrix in the linear space $\mathbb{K}^{n \times n}$ of $n \times n$ matrices over $\mathbb{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 $\mathbb{K}^{n \times n}$, for various fields (including finite fields), as a product of two matrices. Each coefficient of these matrices lie in an extension field of $\mathbb{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 solving `IP` when $\mathbf{f}=(x_1^d,\ldots,x_n^d)$ is presented. From an algorithmic point of view, the problem boils down to factoring the determinant of a linear matrix (*i.e.* a matrix whose components are linear polynomials). This extends to `IP` a result of Kayal obtained for `PolyProj`.
address:
- 'Sorbonne Universités, <span style="font-variant:small-caps;">UPMC</span> Univ Paris 06, Équipe <span style="font-variant:small-caps;">PolSys</span>, <span style="font-variant:small-caps;">LIP6</span>, F-75005, Paris, France'
- '<span style="font-variant:small-caps;">CNRS</span>, <span style="font-variant:small-caps;">UMR 7606</span>, <span style="font-variant:small-caps;">LIP6</span>, F-75005, Paris, France'
- '<span style="font-variant:small-caps;">INRIA</span>, Équipe <span style="font-variant:small-caps;">PolSys</span>, Centre Paris – Rocquencourt, F-75005, Paris, France'
author:
- Jérémy Berthomieu
- 'Jean-Charles Faugère'
- Ludovic Perret
bibliography:
- 'main.bib'
title: 'Polynomial-Time Algorithms for Quadratic Isomorphism of Polynomials: The Regular Case'
---

Quadratic forms ,computer algebra ,polynomial isomorphism ,multivariate cryptography ,module isomorphism

12Y05 ,94A60 ,68W20 ,68W30 ,68Q25

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.* [@DBLP:conf/eurocrypt/Patarin96; @DBLP:conf/nss/TangX12; @Tang2014; @DBLP:conf/ispec/YangTY11]), arithmetic complexity (*via* Geometric Complexity Theory (GCT) for instance, see [@DBLP:conf/coco/Burgisser12; @DBLP:conf/stoc/Kayal12; @DBLP:journals/cacm/Mulmuley12; @DBLP:journals/siamcomp/MulmuleyS01]), testing low degree affine-invariant properties  ([@DBLP:conf/soda/BhattacharyyaFL13; @DBLP:journals/cdm/GreenT09; @DBLP:journals/eccc/BhattacharyyaGRS11], $\ldots$). As we will see, the notion of equivalence can come with different flavours that impact the intrinsic hardness of the problem considered.

In [@DBLP:conf/stacs/AgrawalS06; @Saxena2006], 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 ${{\mathbb}{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 <span style="font-variant:small-caps;">ID</span> scheme of [@DBLP:conf/eurocrypt/Patarin96], and more recently group/proxy signature schemes ([@DBLP:conf/nss/TangX12; @Tang2014; @DBLP:conf/ispec/YangTY11]). Note that there is a subtle difference between the equivalence problem considered in [@DBLP:conf/stacs/AgrawalS06; @Kayal2011; @Saxena2006] and the one considered in cryptographic applications.

Whilst [@DBLP:conf/stacs/AgrawalS06; @Kayal2011; @Saxena2006] restrict their attention to $m=1$, arbitrary $m \geq 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.

Equivalence of multivariate polynomials is also a fundamental problem in Multivariate Public-Key Cryptography ([`MPKC`]{}). This is a family of asymmetric (encryption and signature) schemes whose public-key is given by a set of $m$ multivariate equations ([@C*; @DBLP:conf/eurocrypt/Patarin96]). To minimize the public-key storage, the multivariate polynomials considered are usually quadratic. The basic idea of [`MPKC`]{}is to construct a public-key which is equivalent to a set of quadratic multivariate polynomials with a specific structure (see for instance [@WP11]). Note that the notion of equivalence considered in this context is more general than the one considered for [`PolyEquiv`]{}or [`IP1S`]{}. Indeed, the equivalence is induced by an invertible linear change of variables and an invertible linear combination on the polynomials. The corresponding equivalence problem is known ([@DBLP:conf/eurocrypt/Patarin96]) as *Isomorphism of Polynomials* ([`IP`]{}or `IP2S`).

[`PolyEquiv`]{}, [`IP`]{}, and [`IP1S`]{}are not NP-Hard unless the polynomial-hierarchy collapses, [@DBLP:journals/eccc/ECCC-TR04-116; @DBLP:conf/eurocrypt/PatarinGC98]. However, the situation changes drastically when considering the equivalence for more general linear transformations (in particular, not necessarily invertible). In this context, the problem is called [`PolyProj`]{}. At [ SODA’11]{}, [@Kayal2011] showed that [`PolyProj`]{}is NP-Hard. This may be due to the fact that various fundamental questions in arithmetic complexity can be re-interpreted as particular instances of [`PolyProj`]{}(see [@DBLP:conf/coco/Burgisser12; @DBLP:conf/stoc/Kayal12; @DBLP:journals/cacm/Mulmuley12; @DBLP:journals/siamcomp/MulmuleyS01]).

Typically, the famous <span style="font-variant:small-caps;">VP</span> vs <span style="font-variant:small-caps;">VNP</span> question ([@DBLP:journals/tcs/Valiant79]) can be formulated as an equivalence problem between the determinant and permanent polynomials. Such a link is in fact the core motivation of Geometric Complexity Theory. The problem of computing the symmetric rank ([@DBLP:journals/jsc/BernardiGI11; @DBLP:journals/siammax/ComonGLM08]) of a symmetric tensor also reduces to an equivalence problem involving a particular multivariate polynomial ([@DBLP:conf/stoc/Kayal12]). To mention another fundamental problem, the task of minimizing the cost of computing matrix multiplication reduces to a particular equivalence problem ([@DBLP:conf/stoc/BurgisserI11; @DBLP:conf/stoc/BurgisserI13; @DBLP:conf/soda/CohnU13; @DBLP:conf/stoc/Kayal12]).

Organization of the Paper and Main Results {#organization-of-the-paper-and-main-results .unnumbered}
------------------------------------------

Let ${{\mathbb}{K}}$ be a field, ${{\mathbf}{f}}$ and ${{\mathbf}{g}}$ be two sets of $m$ polynomials over ${{\mathbb}{K}}[x_1,\ldots,x_n]$. The Isomorphism of Polynomials ([`IP`]{}) problem, introduced by Patarin ([@DBLP:conf/eurocrypt/Patarin96]), is as follows:\
**Isomorphism of Polynomials ([`IP`]{})**\
**Input:** $\big( ({{\mathbf}{f}}=(f_1,\ldots,f_m)$ and ${{\mathbf}{g}}=(g_1,\ldots,g_m)\big) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m$.\
**Question:** Find – if any – $(A,B) \in\operatorname{GL}_n({{\mathbb}{K}}) \times
\operatorname{GL}_m({{\mathbb}{K}})$ such that: $${{\mathbf}{g}}({{\mathbf}{x}})=B\cdot {{\mathbf}{f}}(A \cdot {{\mathbf}{x}}), \mbox{ with
  ${{\mathbf}{x}}=(x_1,\dots,x_n)^{\operatorname{T}}$}.$$ While [`IP`]{}is a fundamental problem in multivariate cryptography, there are quite few algorithms, such as [@DBLP:conf/eurocrypt/PatarinGC98; @DBLP:conf/eurocrypt/BouillaguetFV13; @FaugerePerret06], solving [`IP`]{}. In particular, [@FaugerePerret06] proposed to solve [`IP`]{}by reducing it to a system of nonlinear equations whose variables are the unknown coefficients of the matrices. It was conjectured in [@FaugerePerret06], but never proved, that the corresponding system of nonlinear equations can be solved in polynomial time as soon as the [`IP`]{}instances considered are not homogeneous. 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, [@DBLP:conf/eurocrypt/BouillaguetFV13] 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 \[prop:homo\]). As a side remark, we mention that there exist more efficient methods to handle the affine case; typically by considering the homogeneous components, see [@FaugerePerret06]. 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 [@DBLP:conf/pkc/BouillaguetFFP11], 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, [@MacariotRatPG2013] 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 ${{\mathbf}{f}}\in{{\mathbb}{K}}[{{\mathbf}{x}}]^m$. That is, compute $A\in\operatorname{GL}_n({{\mathbb}{K}})$ such that ${{\mathbf}{g}}({{\mathbf}{x}})={{\mathbf}{f}}(A\cdot{{\mathbf}{x}})$ is easier to solve. In this setting, <span style="font-variant:small-caps;">Ridge</span> algorithm (see [@BerthomieuHM2010; @Hironaka1970; @Giraud1972]) and <span style="font-variant:small-caps;">MinVar</span> algorithm (see [@Carlini05reducingthe; @Kayal2011]) reduce to the best the number of variables of the system. More generally, for a given homogeneous polynomial system ${{\mathbf}{f}}$, the *Functional Decomposition Problem* is the problem of computing $\mathbf{h}=(h_1,\ldots,h_s)$ homogeneous and ${{\mathbf}{g}}$ such that ${{\mathbf}{f}}({{\mathbf}{x}})={{\mathbf}{g}}(\mathbf{h}({{\mathbf}{x}}))$.

To simplify the presentation in this introduction, we mainly deal with fields of characteristic $\not =2$. Results for fields of characteristic $2$ are also given later in this paper. Now, we define formally [`IP1S`]{}:

Let $\big( {{\mathbf}{f}}=(f_1,\ldots,f_m),{{\mathbf}{g}}=(g_1,\ldots,g_m)\big) \in
  {{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m$. We shall say that ${{\mathbf}{f}}$ and ${{\mathbf}{g}}$ are equivalent, denoted ${{\mathbf}{f}}\sim {{\mathbf}{g}}$, if $\exists$ $A
  \in \operatorname{GL}_n({{\mathbb}{K}})$ such that: $${{\mathbf}{g}}({{\mathbf}{x}})={{\mathbf}{f}}(A \cdot {{\mathbf}{x}}).$$ [`IP1S`]{}is then the problem of finding – if any – $A \in \operatorname{GL}_n({{\mathbb}{K}})$ that makes ${{\mathbf}{g}}$ equivalent to ${{\mathbf}{f}}$ (*i.e.* $A \in \operatorname{GL}_n({{\mathbb}{K}})$ such that ${{\mathbf}{g}}({{\mathbf}{x}})={{\mathbf}{f}}(A \cdot {{\mathbf}{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`]{} {#canonical-form-of-ip1s .unnumbered}

Let $\big( {{\mathbf}{f}}=(f_1,\ldots,f_m),{{\mathbf}{g}}=(g_1,\ldots,g_m)\big) \in
{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m$ be homogeneous quadratic polynomials. Let $H_1,\ldots,H_m$ be the Hessian matrices of $f_1,\ldots,f_m$ (resp. $H^{\prime}_1,\ldots,H^{\prime}_m$ be the Hessian matrices of $g_1,\ldots,g_m$). Recall that the Hessian matrix associated to a $f_i$ is defined as $H_i=\left(\frac{\partial^2
    f_i}{\partial x_k \partial x_{\ell}}\right)_{k,\ell}\in
\mathbb{K}^{n\times n}$. Consequently, [`IP1S`]{}for quadratic forms is equivalent to finding $A \in \operatorname{GL}_n({{\mathbb}{K}})$ such that: $$\label{eq:orig}
  H_i^{\prime}=A^{\operatorname{T}}\cdot H_i \cdot A,   \text{ for all } i, 1 \leq i \leq m.$$ Assuming $H_j$ is invertible, and thus so is $H_j'$, one has $H_j'^{-1}=A^{-1}H_j^{-1}A^{-\operatorname{T}}$. Combining this with equation  yields $H_j'^{-1}H_i'=A^{-1}\cdot
H_j^{-1}H_i\cdot 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:

 \[as:regular\] Let ${{\mathbf}{f}}=(f_1,\ldots,f_m)\in{{\mathbb}{K}}[x_1,\ldots,x_n]$. We assume that a linear combination over ${{\mathbb}{K}}$ of the quadratic forms $f_1,\ldots,f_m$ is not degenerate[^1]. In particular, we assume that $|{{\mathbb}{K}}| > n$.

Taking as variables the entries of $A$, we can see that  naturally yields a nonlinear system of equations. However, we show that one can essentially linearize equations . To this end, we prove in Section \[s:normalization\] that under Assumption \[as:regular\] 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:

\[th:normal\] Let ${{\mathbb}{K}}$ be a field of $\operatorname{char}{{\mathbb}{K}}\neq 2$. There exists a randomized polynomial-time algorithm which given a regular quadratic homogeneous instance of [`IP1S`]{}returns “[<span style="font-variant:small-caps;">NoSolution</span>]{}” only if the two systems are not equivalent or a canonical form $${\left(\big(\sum_{i=1}^{n}d_ix_{i}^{2},f_2,\ldots,f_m\big),
    \big(\sum_{i=1}^{n}d_ix_{i}^{2},g_2,\ldots,g_m\big)\right)},$$ where the $d_i$ are equal to $1$ or a nonsquare in ${{\mathbb}{K}}$, $f_i$ and $g_i$ are *nondegenerate* homogeneous quadratic polynomials in ${{\mathbb}{K}}[x_1,\ldots,x_n]$. Any solution on ${{\mathbb}{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 [@Saxena2006 Chapter 3], [@Wallenborn2013 Chapter 1]. As such, one could only assume that the first quadratic form of the second set is $\sum_{i=1}^nd_i'x_i^2$. This does not fundamentally change the algorithms presented in this paper, beside some matrices denoted by $D$ which could be changed into $D'=\operatorname{Diag}(d_1',\ldots,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\in{{\mathbb}{L}}^{n\times n}$ with ${{\mathbb}{L}}$ an extension of ${{\mathbb}{K}}$, one can amplify the success probability over a small field by using the fact that matrices are conjugate over ${{\mathbb}{K}}$ if and only if they are conjugate over an algebraic extension ${{\mathbb}{L}}$ (see [@Pazzis2010]). Thus, one can search linear change of variables with coefficients in some algebraic extension ${{\mathbb}{L}}\supseteq{{\mathbb}{K}}$ (but of limited degree).

### Conjugacy Problem {#conjugacy-problem .unnumbered}

When [`IP1S`]{}is given in canonical form, equations  can be rewritten as $A^{\operatorname{T}}\,D\, A = D$ with $D=\operatorname{Diag}(d_1,\ldots,d_n)$ and $H_i^{\prime}=A^{\operatorname{T}}\cdot H_i \cdot A=D\,A^{-1}\,D^{-1}\cdot H_i \cdot
A$ for all $i, 2 \leq i \leq m$. Our task is now to solve the following problem:

Let ${{\mathbb}{K}}^{n\times n}$ be the set of $n\times n$ matrices with entries in ${{\mathbb}{K}}$. Let $\{H_1,\ldots,H_m\}$ and $\{H^{\prime}_1,\ldots,H^{\prime}_m\}$ be two families of matrices in ${{\mathbb}{K}}^{n\times n}$. The [$D$-`OSMC`]{}problem is the task to recover – if any – a $D$-orthogonal matrix $X \in {{\mathbb}{L}}^{n\times n}$, *i.e.* $X^{\operatorname{T}}DX=D$, with ${{\mathbb}{L}}$ being an algebraic extension of ${{\mathbb}{K}}$, such that: $$\begin{aligned}
    X^{-1}\, H_i \, X=H^{\prime}_i, & \quad \forall\,i, 1 \leq i
    \leq m,\end{aligned}$$

[@DBLP:conf/issac/ChistovIK97] show that [$D$-`OSMC`]{}with $D={{\mathrm}{Id}}$ is equivalent to:

1.  Solving the Simultaneous Matrix Conjugacy problem ([`SMC`]{}) between $\{H_{i}\}_{1 \leq i \leq m}$ and $\{H^{\prime}_i\}_{1 \leq
        i \leq m}$, that is to say finding an invertible matrix $Y \in
      \operatorname{GL}_n({{\mathbb}{K}})$ such that: $$\begin{aligned}
    \label{eq:smc0} Y^{-1}\cdot H_i \cdot Y=H^{\prime}_i
        & \text{ and } & Y^{-1}\cdot H_i^{\operatorname{T}} \cdot Y={H^{\prime}_i}^{\operatorname{T}}
        \quad \forall\,i, 1 \leq i \leq m.
      \end{aligned}$$

2.  Computing the square-root $W$ of the matrix $Z=Y \cdot
      Y^{\operatorname{T}}$. Then, the solution of the [$D$-`OSMC`]{}problem is given by $X=Y\,
      W^{-1}$.

In our context, $D=\operatorname{Diag}(d_1,\ldots,d_n)$ is any diagonal invertible matrix. So, we extend [@DBLP:conf/issac/ChistovIK97] and show that [$D$-`OSMC`]{}is equivalent to

1.  Finding an invertible matrix $Y \in \operatorname{GL}_n({{\mathbb}{K}})$ such that: $$\begin{aligned}
     \label{eq:smc} Y^{-1}\cdot H_i \cdot Y=H^{\prime}_i
        & \text{ and } &D\,Y^{-1}D^{-1}\cdot H_i^{\operatorname{T}} \cdot
        D\,YD^{-1}={H^{\prime}_i}^{\operatorname{T}} \quad \forall\,i, 1 \leq i \leq m.
      \end{aligned}$$

2.  Computing the square-root $W$ of the matrix $Z=D\,Y \cdot
      Y^{\operatorname{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^{\prime}_i$’s) are symmetric (Hessian matrices). Thus, condition  yields a system of *linear* equations and one polynomial inequation: $$\label{eq:lin}
  H_1\cdot Y=Y\cdot H^{\prime}_1,\ldots,H_m \cdot Y = Y\cdot H^{\prime}_m
  \text{ and }\det(Y)\neq 0.$$ From now on, we shall denote by ${\mathcal}{O}_n({{\mathbb}{L}}, D)$ the set of $D$-orthogonal matrices with coefficients in ${{\mathbb}{L}}$.

Let $V \subset {{\mathbb}{K}}^{n \times 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 ([@edm67]). Note that, if the representation of the algebra spanned by $\{H_{1}^{-1}H_{i}\}_{1\leq i\leq m} $ is *irreducible*, we know that $V$ has dimension at most $1$ (Schur’s Lemma, see [@Lang2002 Chap. XVII, Proposition 1.1] and [@Newman1967 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 ([@DL78; @Z79]), a solution to [$D$-`OSMC`]{}with overwhelming probability as soon as ${{\mathbb}{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 [@Pazzis2010]). 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 [@DBLP:conf/issac/ChistovIK97] and provide a deterministic polynomial-time algorithm for solving .

### Characteristic $2$ {#characteristic-2 .unnumbered}

Let us recall that in characteristic $2$, the associated matrices $H_1,\ldots,H_m,H_1',\ldots,H_m'$ to quadratic forms can be chosen as upper triangular. In this context, we show in Section \[ss:char\_2\] that [`IP1S`]{}can still be reduced to a $(H_1+H_1^{\operatorname{T}})$-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 \[ss:bench\]. 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 $\leq 30$, for efficiency reasons).

### Matrix Square Root Computation {#matrix-square-root-computation .unnumbered}

It is well known that computing square roots of matrices can be done efficiently using numerical methods (for instance, see [@grant]). On the other hand, it seems difficult to control the bit complexity of numerical methods. In [@DBLP:conf/issac/ChistovIK97 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 [@Cai1994]. However, by applying a small trick to the proof of [@DBLP:conf/issac/ChistovIK97], one can compute a solution in polynomial-time for various field of characteristic $\neq 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 [@DBLP:conf/issac/ChistovIK97 Section 3] is presented in Section \[ss:computation\]. To do so, we adapt the technique of [@Cai1994] 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,\ldots,H_m$ and $H_1',\ldots,H_m'$ are indeed equivalent. We then present a second simpler method based on the *generalized Jordan normal form* (see Section \[ss:sqrt\_char\_p\]) 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:

\[th:probmain\] Let ${{\mathbb}{K}}$ be a field with $\operatorname{char}{{\mathbb}{K}}\neq
  2$. Under Assumption \[as:regular\], there exists a randomized polynomial-time algorithm for solving quadratic-[`IP1S`]{}over an extension field of ${{\mathbb}{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 [@BettaleFP2013], the authors recover an equivalent key over an extension for the multi-HFE scheme.

In addition, under some nondegeneracy assumption, Theorem \[th:probmain\] can be turned into a deterministic algorithm solving [`IP1S`]{}over the base field ${{\mathbb}{K}}$ or an extension thereof. That is:

\[th:main\] Under Assumption \[as:regular\] 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 ${{\mathbb}{K}}$ of polynomial degree in $n$. Furthermore, if the space of matrices satisfying equations  has dimension $1$, then the algorithm can solve quadratic-[`IP1S`]{}over ${{\mathbb}{K}}$.

Let us note that assuming that one of the Hessian matrix is invertible is not a strong assumption when the size of ${{\mathbb}{K}}$ is not too small. Indeed, the probability of picking a random invertible symmetric matrix over ${{\mathbb}{F}}_q$ is $$\frac{\prod_{i=1}^{n}(1-q^{-i})}
{\prod_{i=1}^{\lfloor n/2\rfloor}(1-q^{-2i})}
=\prod_{i=1}^{\lceil n/2\rceil}(1-q^{-2i+1}),$$ see [@Carlitz1954 Equations 4.7 and 4.8].

If $m\geq 3$, for random matrices $H_1,\ldots,H_m$, the set of solutions of equations  is a $1$-dimensional matrix space. This allows us to solve quadratic-[`IP1S`]{}in polynomial-time over ${{\mathbb}{K}}$. In Section \[ss:bench\], we present our timings for solving [`IP1S`]{}. These experiments confirm that for randomly chosen matrices, our method solves [`IP1S`]{}over ${{\mathbb}{K}}$. We remark also that our method succeeds to solve [`IP1S`]{}over ${{\mathbb}{F}}_2$ for public-keys whose sizes are much bigger than practical ones.

In Section \[s:enumeration\], we consider the counting problem \#[`IP1S`]{}associated to [`IP1S`]{}for quadratic (homogeneous) polynomials in its canonical form (as defined in Theorem \[th:normal\]). 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 [@DFPW11; @WP11]. Given homogeneous quadratic polynomials $\big( {{\mathbf}{f}}=(f_1,\ldots,f_m),{{\mathbf}{g}}=(g_1,\ldots,g_m)\big) \in
{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m$, we want to count the number of invertible matrices $A \in \operatorname{GL}_n({{\mathbb}{K}})$ such that $
{{\mathbf}{g}}({{\mathbf}{x}})={{\mathbf}{f}}(A \cdot {{\mathbf}{x}}).$ To do so, we define:

Let ${{\mathbf}{f}}=(f_1,\ldots,f_m) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^m$, we shall call *automorphism group of ${{\mathbf}{f}}$* the set: $$\mathcal{G}_{{{\mathbf}{f}}}=\{ A \in \operatorname{GL}_n({{\mathbb}{K}}) \mid {{\mathbf}{f}}(A\cdot {{\mathbf}{x}}) ={{\mathbf}{f}}({{\mathbf}{x}})
  \}.$$

If ${{\mathbf}{f}}\sim {{\mathbf}{g}}$, the automorphism groups of ${{\mathbf}{f}}$ and ${{\mathbf}{g}}$ are similar. Thus, the size of the automorphism group of ${{\mathbf}{f}}$ allows us to count the number of invertible matrices mapping ${{\mathbf}{f}}$ to ${{\mathbf}{g}}$. For quadratic homogeneous polynomials, the automorphism group coincides with the subset of regular matrices in the centralizer $\operatorname{{\mathcal}{C}}({{\mathcal}{H}})$ of the Hessian matrices ${\mathcal}{H}$ associated to ${{\mathbf}{f}}$. Taking $\alpha$ an algebraic element of degree $m$ over ${{\mathbb}{K}}={{\mathbb}{F}}_q$, let us assume the Jordan normal form of $H=\sum_{i=i}^m\alpha^{i-1}\,H_i$ has Jordan blocks of sizes $s_{i,1}\leq\cdots\leq s_{i,d_i}$ associated to eigenvalue $\zeta_i$, for $i,1\leq i\leq r$. Then, as a consequence of [@Singla2010 Lemma 4.11], we prove that, if $q$ is an odd prime power, then the number of solutions of quadratic-[`IP1S`]{}in ${{\mathbb}{F}}_q^{n\times n}$ is bounded from above by: $$q^{{\left(\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i}(2\,d_i-2\,j+1)s_{i,j}\right)}}-1.$$

Open Question: The Irregular Case {#open-question-the-irregular-case .unnumbered}
---------------------------------

Given a quadratic instance of [`IP1S`]{}, a nondegenerate instance is an instance wherein the matrix whose rows are all the rows of $H_1,\ldots,H_m$ has rank $n$. In paragraph \[par:variables\], 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,\ldots,H_m$ does not contain a nondegenerate matrix. This situation is well illustrated by the following example $f_1=x_1x_3,f_2=x_2x_3$. Any linear combination of $f_1,f_2$ is degenerate, while ${{\mathbf}{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`]{} ([@DBLP:conf/pkc/BouillaguetFFP11; @MacariotRatPG2013]) will fail on these instances; making the hardness of the irregular case intriguing and then an interesting open question.

Special case of [`IP`]{} {#special-case-of-ip .unnumbered}
------------------------

In our quest of finding instances of [`IP`]{}solvable in polynomial-time, we take a first step in Section \[s:ip\]. We consider [`IP`]{}for a specific set of polynomials with $m=n$. In the aforementioned Section \[s:ip\], we prove the following:

\[th:ip\] Let ${{\mathbf}{g}}=(g_1,\ldots,g_n) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$ be given in dense representation, and ${{\mathbf}{f}}={{\mathrm}{{\mathbf}{POW}}}_{n,d}=(x_1^d,\ldots,x_n^d) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$ for some $d>0$. Whenever $\operatorname{char}{{\mathbb}{K}}=0$, let $e=d$ and $\tilde{{{\mathbf}{g}}}={{\mathbf}{g}}$. Otherwise, let $p=\operatorname{char}{{\mathbb}{K}}$, let $e$ and $r$ be integers such that $d=p^r\,e$, $p$ and $e$ coprime, and let $\tilde{{{\mathbf}{g}}}\in{{\mathbb}{K}}[{{\mathbf}{x}}]^n$ be such that ${{\mathbf}{g}}({{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}}^{p^r})$. Let $L$ be a polynomial size of an arithmetic circuit to evaluate the determinant of the Jacobian matrix of $\tilde{{{\mathbf}{g}}}$. If the size of ${{\mathbb}{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)\in\operatorname{GL}_n({{\mathbb}{K}})\times\operatorname{GL}_n({{\mathbb}{K}})$ such that: $${{\mathbf}{g}}=B\cdot{{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\cdot{{\mathbf}{x}}).$$

This extends a similar result of [@Kayal2011 Section 5] who considered [`PolyEquiv`]{}for a sum of $d$-power polynomials. We show that solving [`IP`]{}for ${{\mathrm}{{\mathbf}{POW}}}_{n,d}$ reduces to factoring the determinant of a Jacobian matrix (in [@Kayal2011], the Hessian matrix is considered). This illustrates, how powerful partial derivatives can be in equivalence problems ([@DBLP:journals/fttcs/ChenKW11; @DBLP:conf/eurocrypt/Perret05]). To go along with the proof of Theorem \[th:ip\], we design Algorithm \[al:ip\] at the end of Section \[s:ip\].

Normalization - Canonical form of [`IP1S`]{} {#s:normalization}
============================================

In this section, we prove Theorem \[th:normal\]. In other words, we explain how to reduce, under Assumption \[as:regular\], any quadratic homogeneous instance $({{\mathbf}{f}},{{\mathbf}{g}}) \in
{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,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.

#### Homogenization {#par:homo}

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\in{{\mathbb}{K}}[{{\mathbf}{x}}]$ of degree $2$, we denote by $p^{\star}(x_0,x_1,\ldots,x_n)=x_0^2\,p(x_1/x_0,\ldots, x_n/x_0)$ its *homogenization*.

\[prop:homo\] ${\texttt{IP1S}\xspace}$ 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.

Let $({{\mathbf}{f}},{{\mathbf}{g}})\in{{\mathbb}{K}}[{{\mathbf}{x}}]^m \times {{\mathbb}{K}}[{{\mathbf}{x}}]^m$ be inhomogeneous polynomials of degree $2$. We consider the transformation which maps $({{\mathbf}{f}},{{\mathbf}{g}})$ to $\big({{\mathbf}{f}}^{\star}=(f_0^{\star}=x_0^2,f_1^{\star},\ldots,f_m^{\star}),
  {{\mathbf}{g}}^{\star}=(g_0^{\star}=x_0^2,g_1^{\star},\ldots,g_m^{\star})\big)$. This clearly transforms polynomials of degree $2$ to homogeneous quadratic polynomials. We can write $f_i({{\mathbf}{x}})={{\mathbf}{x}}^{\operatorname{T}}\,H_i\,{{\mathbf}{x}}+L_i\,{{\mathbf}{x}}+c_i$ with $H_i\in{{\mathbb}{K}}^{n\times n}$, $L_i\in{{\mathbb}{K}}^n$ and $c_i\in{{\mathbb}{K}}$, then $f_i(A{{\mathbf}{x}}+b)=(A{{\mathbf}{x}}+b)^{\operatorname{T}}\,H_i\,(A{{\mathbf}{x}}+b)+L_i\,(A{{\mathbf}{x}}+b)+c_i$ and its homogenization is $(A{{\mathbf}{x}}+bx_0)^{\operatorname{T}}\,H_i\,(A{{\mathbf}{x}}+bx_0)+L_i\,(A{{\mathbf}{x}}+bx_0)\,x_0+c_i\,x_0^2
  =f^{\star}_i(A'{{\mathbf}{x}}^{\star})$, with ${{\mathbf}{x}}^{\star}=(x_0,x_1,\ldots,x_n)^{\operatorname{T}}$. If $(A,b) \in
  \operatorname{GL}_{n}({{\mathbb}{K}}) \times {{\mathbb}{K}}^n$ is an affine transformation solution on the inhomogeneous instance then $A^{\prime}={\left(\begin{smallmatrix}1 &{\mathbf}{0}\\
      b &A\end{smallmatrix}\right)}$ is a solution for the homogenized instance. Conversely, a solution $A' \in \operatorname{GL}_{n+1}({{\mathbb}{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_0^2$ and $g_0=x_0^2$ and setting $C'=A'/a_{0,0}'$, with $a_{0,0}'=\pm 1$. One can see that $C'$ is of the form ${\left(\begin{smallmatrix}1 &{\mathbf}{0}\\d &C\end{smallmatrix}\right)}$, and $(C,d)\in \operatorname{GL}_n({{\mathbb}{K}}) \times {{\mathbb}{K}}^n$ is a solution for $({{\mathbf}{f}},{{\mathbf}{g}})$.

#### Redundant Variables {#par:variables}

As a first preliminary natural manipulation, we first want to eliminate – if any – *redundant variables* from the instances considered. Thanks to [@Carlini05reducingthe] (and reformulated in [@Kayal2011]), this task can be done in randomized polynomial time:

([@Carlini05reducingthe; @Kayal2011]) \[prop:carl\] Let $f \in {{\mathbb}{K}}[x_1\ldots,x_n]$ be a polynomial. We shall say that $f$ has $s$ *essential variables* if $\exists \,M \in \operatorname{GL}_n({{\mathbb}{K}})$ such that $f(M
  {{\mathbf}{x}})$ depends only on the first $s$ variables $x_1,\ldots,x_s$. The remaining $n-s$ variables $x_{s+1}\ldots,x_n$ will be called *redundant* variables. If $\operatorname{char}{{\mathbb}{K}}=0$ or $\operatorname{char}{{\mathbb}{K}}>\deg f$, and $f$ has $s$ essential variables, then we can compute in randomized polynomial time $M \in \operatorname{GL}_n({{\mathbb}{K}})$ such that $f(M\,{{\mathbf}{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.

\[def:ess\] The number of *essential variables* of ${{\mathbf}{f}}=(f_1,\ldots,f_m)\in{{\mathbb}{K}}[x_1,\ldots,x_n]^m$ is the smallest $s$ such that ${{\mathbf}{f}}$ can be decomposed as: $${{\mathbf}{f}}=\tilde{{{\mathbf}{f}}}(\ell_1,\ldots,\ell_{s})$$ with $\ell_1,\ldots,\ell_{s}$ being linear forms in $x_1,\ldots,x_n$ of rank $s$ and $\tilde{{{\mathbf}{f}}}\in{{\mathbb}{K}}[y_1,\ldots,y_{s}]^m$.

The linear forms $\ell_1,\ldots,\ell_s$ can be easily computed thanks to Proposition \[prop:carl\] when the characteristic of ${{\mathbb}{K}}$ is zero or greater than the degrees of $f_1,\ldots,f_m$. In characteristic $2$, when ${{\mathbb}{K}}$ is perfect (which is always true if ${{\mathbb}{K}}$ is finite for instance) the linear forms can also be recovered in polynomial time (see [@BerthomieuHM2010; @Giraud1972; @Hironaka1970] for instance). Below, we show that we can restrict our attention to only essential variables. Namely, solving [`IP1S`]{}on $({{\mathbf}{f}},{{\mathbf}{g}})$ reduces to solving [`IP1S`]{}on instances having only essential variables.

\[lm:rank\] Let $({{\mathbf}{f}},{{\mathbf}{g}})\in{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m $ be two sets of quadratic polynomials. If ${{\mathbf}{f}}\sim {{\mathbf}{g}}$, then their numbers of essential variables must be the same. Let $s$ be the number of essential variables of ${{\mathbf}{f}}$. Finally, let $(\tilde{{{\mathbf}{f}}},\tilde{{{\mathbf}{g}}})\in{{\mathbb}{K}}[y_1,\ldots,y_s]^m \times
  {{\mathbb}{K}}[y_1,\ldots,y_s]^m$ be such that: $${{\mathbf}{f}}=\tilde{{{\mathbf}{f}}}(\ell_1,\ldots,\ell_{s}) \mbox{ and }
  {{\mathbf}{g}}=\tilde{{{\mathbf}{g}}}(\ell'_1,\ldots,\ell'_{s}),$$ with $\ell_1,\ldots,\ell_{s}$ (resp. $\ell'_1,\ldots,\ell'_{s}$) linear forms in ${{\mathbf}{x}}$ of rank $s$ and $\tilde{{{\mathbf}{f}}},
  \tilde{{{\mathbf}{g}}}\in{{\mathbb}{K}}[y_1,\ldots,y_{s}]^m$. It holds that: $${{\mathbf}{f}}\sim {{\mathbf}{g}}\iff \tilde{{{\mathbf}{f}}} \sim \tilde{{{\mathbf}{g}}}.$$

Let $H_1,\ldots,H_m$ be the Hessian matrices of $f_1,\ldots,f_m$ (resp. $H^{\prime}_1,\ldots,H^{\prime}_m$ be the Hessian matrices of $g_1,\ldots,g_m$). Similarly, we define the Hessian matrices $\tilde{H}_1,\ldots,\tilde{H}_m$ (resp. $\tilde{H}^{\prime}_1,\ldots,\tilde{H}^{\prime}_m$) of $\tilde{f}_1,\ldots,\tilde{f}_m$ (resp. $\tilde{g}_1,\ldots,\tilde{g}_m$). Let also $M$ and $N$ be matrices in $\operatorname{GL}_n({{\mathbb}{K}})$ such that $H_i=M^{\operatorname{T}}{\left(\begin{smallmatrix} \tilde{H}_i & {\mathbf}{0}\\{\mathbf}{0}
      &{\mathbf}{0}\end{smallmatrix}\right)}M$ and $H_i'=N^{\operatorname{T}}{\left(\begin{smallmatrix}
      \tilde{H}_i' &{\mathbf}{0}\\{\mathbf}{0} &{\mathbf}{0}\end{smallmatrix}\right)}N$ for all $i,1\leq i\leq
  m$. There exist such $M$ and $N$, as ${{\mathbf}{f}}$ and ${{\mathbf}{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={\left(\begin{smallmatrix}M_1&M_2\\{\mathbf}{0}&{{\mathrm}{Id}}\end{smallmatrix}\right)}$ and $N={\left(\begin{smallmatrix}N_1&N_2\\{\mathbf}{0}&{{\mathrm}{Id}}\end{smallmatrix}\right)}$, with $M_1,N_1\in\operatorname{GL}_s({{\mathbb}{K}})$.

If $\tilde{{{\mathbf}{f}}} \sim \tilde{{{\mathbf}{g}}}$, $\exists \, \tilde{A}\in\operatorname{GL}_s({{\mathbb}{K}})$ such that $A^{\operatorname{T}}\tilde{H}_i\tilde{A}=\tilde{H}_i'$, for all $i,1\leq i\leq m$. Then, for all $B\in{{\mathbb}{K}}^{(n-s)\times s}$ and $C\in\operatorname{GL}_{n-s}({{\mathbb}{K}})$: $$\begin{aligned}
    {\left(\begin{smallmatrix}\tilde{A}^{\operatorname{T}} &B^{\operatorname{T}}\\ {\mathbf}{0}
        &C^{\operatorname{T}}\end{smallmatrix}\right)}
  {\left(\begin{smallmatrix}\tilde{H}_i\vphantom{^{\operatorname{T}}}&{\mathbf}{0}\\{\mathbf}{0} &{\mathbf}{0}\vphantom{^{\operatorname{T}}}
    \end{smallmatrix}\right)}
  {\left(\begin{smallmatrix}\tilde{A}\vphantom{^{\operatorname{T}}} &{\mathbf}{0}\\ B &C\vphantom{^{\operatorname{T}}}
    \end{smallmatrix}\right)}
  &=
  {\left(\begin{smallmatrix}\tilde{H}_i^{\prime}\vphantom{^{\operatorname{T}}}&{\mathbf}{0}\\{\mathbf}{0}
      &{\mathbf}{0}\vphantom{^{\operatorname{T}}}
    \end{smallmatrix}\right)},\\
  N^{\operatorname{T}}{\left(\begin{smallmatrix}\tilde{A}^{\operatorname{T}} &B^{\operatorname{T}}\\ {\mathbf}{0}
      &C^{\operatorname{T}}\end{smallmatrix}\right)} M^{-\operatorname{T}}H_i\,M^{-1}
  {\left(\begin{smallmatrix}\tilde{A}\vphantom{^{\operatorname{T}}} &{\mathbf}{0}\\ B
      &C\vphantom{^{\operatorname{T}}}
    \end{smallmatrix}\right)}N &=H_i^{\prime}.\end{aligned}$$ Therefore, ${{\mathbf}{f}}$ and ${{\mathbf}{g}}$ are equivalent.

Conversely, we assume now that ${{\mathbf}{f}}\sim {{\mathbf}{g}}$, *i.e.* there exists $A\in\operatorname{GL}_n({{\mathbb}{K}})$ such that $A^{\operatorname{T}}\cdot H_i \cdot A=H^{\prime}_i$, for all $i,1 \leq i \leq m$. This implies that: $$N^{-\operatorname{T}}A^{\operatorname{T}}M^{\operatorname{T}}{\left(\begin{smallmatrix}\tilde{H}_i &{\mathbf}{0}\\{\mathbf}{0} &{\mathbf}{0}
  \end{smallmatrix}\right)}M\,A\,N^{-1}={\left(\begin{smallmatrix}\tilde{H}^{\prime}_i
    &{\mathbf}{0}\\{\mathbf}{0} &{\mathbf}{0}
  \end{smallmatrix}\right)}, \forall\, i,1 \leq i \leq m.$$ We then define $\tilde{A}=((MAN^{-1})_{i,j})_{1\leq i,j\leq s}$, so that $\tilde{{{\mathbf}{f}}}(\tilde{A}{{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}})$. As ${{\mathbf}{g}}$ has $s$ essential variables, then $\operatorname{rank}\tilde{A}$ cannot be smaller than $s$, hence $\tilde{A} \in \operatorname{GL}_s({{\mathbb}{K}})$. We then get $\tilde{A}^{\operatorname{T}}\tilde{H}_i\tilde{A}=\tilde{H}_i'$ for all $i,1 \leq i
\leq m$, *i.e.* $\tilde{{{\mathbf}{f}}} \sim \tilde{{{\mathbf}{g}}}$.

According to Proposition \[lm:rank\], there is an efficient reduction mapping an instance $({{\mathbf}{f}},{{\mathbf}{g}})$ of [`IP1S`]{}to an instance $(\tilde{{{\mathbf}{f}}},\tilde{{{\mathbf}{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 ${{\mathbf}{f}}$ and ${{\mathbf}{g}}$.

#### Canonical Form {#par:identity}

We now assume that $\operatorname{char}{{\mathbb}{K}}\neq 2$.

Let ${{\mathbf}{f}}=(f_1,\ldots,f_m)\in{{\mathbb}{K}}[x_1,\ldots,x_n]^m$ be quadratic homogeneous forms with Hessian matrices $H_1,\ldots,H_m$. We shall say that ${{\mathbf}{f}}$ is *regular* if its number of essential variables is $n$ and if $\exists\,\lambda_1,\dots,\lambda_m\in{{\mathbb}{K}}$ such that $\det{\left(\sum_{i=1}^m\lambda_i\,H_i\right)}\neq 0$.

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 ${{\mathbf}{f}}=(x_1x_3,x_2x_3)$ has $3$ essential variables, but any nonzero linear combination $\lambda_1 f_1+\lambda_2 f_2$ has only $2$ essential variables $\lambda_1 x_1+\lambda_2 x_2$ and $x_3$. Similarly, ${{\mathbf}{f}}=(x_1^2+x_2^2+x_3^2,x_2^2+2x_3^2+x_4^2)$ has $4$ essential variables but any nonzero linear combination $\lambda_1f_1+\lambda_2f_2$ over ${{\mathbb}{F}}_3$ has only $3$ essential variables. This explains the additional condition on the previous definition, and our Assumption \[as:regular\].

We are now in a position to reduce quadratic homogeneous instances of [`IP1S`]{}to a first simplified form.

\[prop:first\_reduction\] Let $({{\mathbf}{f}},{{\mathbf}{g}})\in{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m $ be regular quadratic homogeneous polynomials. There is a randomized polynomial-time algorithm which returns “[<span style="font-variant:small-caps;">NoSolution</span>]{}” only if ${{\mathbf}{f}}\not \sim {{\mathbf}{g}}$, or a new instance $$(\tilde{{{\mathbf}{f}}},\tilde{{{\mathbf}{g}}})={\left(\big(\sum_{i=1}^{n}d_ix_{i}^{2},\tilde{f}_2,
    \ldots,\tilde{f}_m\big),\big(\sum_{i=1}^{n}d_ix_{i}^{2},\tilde{g}_2,\ldots,
    \tilde{g}_m\big)\right)} \in {{\mathbb}{K}}[{{\mathbf}{x}}]^m \times {{\mathbb}{K}}[{{\mathbf}{x}}]^m,$$ with $d_1,\ldots,d_n$ being $1$ or nonsquares in ${{\mathbb}{K}}$, such that ${{\mathbf}{f}}\sim{{\mathbf}{g}}\iff \tilde{{{\mathbf}{f}}} \sim \tilde{{{\mathbf}{g}}}$. If ${{\mathbb}{K}}$ is finite, the output of this algorithm is correct with probability at least $1-n/|{{\mathbb}{K}}|$. If $\tilde{{{\mathbf}{f}}} \sim \tilde{{{\mathbf}{g}}}$, invertible matrices $P, Q$ and $A' \in \operatorname{GL}_n({{\mathbb}{K}})$ are returned such that ${{\mathbf}{f}}(P{{\mathbf}{x}})=\tilde{{{\mathbf}{f}}}({{\mathbf}{x}})$, ${{\mathbf}{g}}(Q{{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}})$ and $\tilde{{{\mathbf}{f}}}(A^{\prime}{{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}})$. It then holds that ${{\mathbf}{f}}(PA'Q^{-1}{{\mathbf}{x}})={{\mathbf}{g}}({{\mathbf}{x}})$.

Let $H_1,\ldots,H_m$ be the Hessian matrices associated to $f_1,\ldots,f_m$. According to Schwartz-Zippel-DeMillo-Lipton Lemma ([@DL78; @Z79]), we can compute in randomized polynomial time $\lambda_1,\ldots,\lambda_m \in {{\mathbb}{K}}$ such that ${\varphi}=\sum_{i=1}^{m}\lambda_i \cdot f_i$ is regular, *i.e.* $\det{\left(\sum_{i=1}^m\lambda_i\,H_i\right)}\neq 0$. The probability to pick $(\lambda_1,\ldots,\lambda_m)\in{{\mathbb}{K}}^m$ on which ${\varphi}$ is regular is bounded from above by $n/|{{\mathbb}{K}}|$. We define $\gamma=\sum_{i=1}^{m}\lambda_i \cdot g_i$. Should one reorder the equations, we can assume w.lo.g. that $\lambda_1\neq 0$. We have then: $${{\mathbf}{f}}\sim {{\mathbf}{g}}\iff ({\varphi},f_2,\ldots,f_m) \sim (\gamma,g_2,\ldots,g_m).$$ Now, applying [Gauß]{}’s reduction algorithm to ${\varphi}$, there exists $d_1,\ldots,d_n\in{{\mathbb}{K}}$, each being $1$ or a nonsquare, such that ${\varphi}=\sum_{i=1}^n{d_i\,\ell_i^2}$, where $\ell_1,\ldots,\ell_n$ are independent linear forms in $x_1,\ldots,x_n$. This gives a $P \in\operatorname{GL}_n({{\mathbb}{L}})$ such that $\tilde{{{\mathbf}{f}}}=(\tilde{{\varphi}}=\sum_{i=1}^nd_i\,x_i^2,\tilde{f_2},\ldots,
  \tilde{f_m})=({\varphi}(P{{\mathbf}{x}}),f_2(P{{\mathbf}{x}}),\ldots,f_m(P{{\mathbf}{x}}))$. Clearly, ${{\mathbf}{f}}\sim\tilde{{{\mathbf}{f}}}$, hence, $\tilde{{{\mathbf}{f}}}\sim{{\mathbf}{g}}$.

After that, we can apply once again [Gauß]{}’s reduction algorithm to $\gamma$. If the reduced polynomial is different from $\sum_{i=1}^{n}d_i\,x_i^2$, then ${{\mathbf}{f}}\not \sim {{\mathbf}{g}}$ and we return “[<span style="font-variant:small-caps;">NoSolution</span>]{}”. Otherwise, the reduction is given by a matrix $Q \in
  \operatorname{GL}_n({{\mathbb}{L}})$ such that $\tilde{{{\mathbf}{g}}}=(\tilde{\gamma}=\sum_{i=1}^nd_i\,x_i^2,\tilde{g}_2,
  \ldots\tilde{g}_m)=
  (\gamma(Q{{\mathbf}{x}}),g_2(Q{{\mathbf}{x}}),\ldots,g_m(Q{{\mathbf}{x}}))$ and ${{\mathbf}{g}}\sim\tilde{{{\mathbf}{g}}}$. Thus, $\tilde{{{\mathbf}{f}}}\sim\tilde{{{\mathbf}{g}}}$ if and only if ${{\mathbf}{f}}\sim{{\mathbf}{g}}$.

Now, assume that $ \exists\,A^{\prime} \in
  \operatorname{GL}_n({{\mathbb}{K}})$ such that $\tilde{{{\mathbf}{f}}}(A^{\prime}{{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}})$. Then, ${{\mathbf}{f}}(PA^{\prime}{{\mathbf}{x}})={{\mathbf}{g}}(Q{{\mathbf}{x}})$, *i.e.*  ${{\mathbf}{f}}(PA^{\prime}Q^{-1}{{\mathbf}{x}})={{\mathbf}{g}}({{\mathbf}{x}})$.

Let us recall that whenever ${{\mathbb}{K}}={{\mathbb}{Q}}$, computing the exact same sum of squares for $\tilde{f}_1$ and $\tilde{g}_1$ is difficult, see [@Saxena2006 Chapter 3], [@Wallenborn2013 Chapter 1]. As such, we could only assume that our canonical form is $\tilde{g}_1=\sum_{i=1}^nd_i'x_i^2$. This would merely change the formulation of following Theorem \[th:second\_reduction\].

#### Invertible Hessian Matrices {#par:invertible}

We are now in a position to reduce any regular homogeneous quadratic instances $({{\mathbf}{f}},{{\mathbf}{g}})$ of [`IP1S`]{}to a new form of the instances where all the polynomials are themselves regular assuming we could find one. From Proposition \[prop:first\_reduction\], this is already the case – under randomized reduction – for $f_1$ and thus $g_1$. For the other polynomials, we proceed as follows. For $i,2\leq i\leq m$, if the Hessian matrix $H_i$ of $f_i$ is invertible, then we do nothing. Otherwise, we change $H_i$ into $H_i-\nu_i\,H_1$, with $\nu_i$ not an eigenvalue of $H_i\,H_1^{-1}$. As ${{\mathbb}{K}}$ has at least $n+1$ elements, there exists such a $\nu_i$ in ${{\mathbb}{K}}$. This gives the following result:

\[th:second\_reduction\] Let $({{\mathbf}{f}},{{\mathbf}{g}})\in{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m $ be regular quadratic homogeneous polynomials. There is a randomized polynomial-time algorithm which returns “[<span style="font-variant:small-caps;">NoSolution</span>]{}” only if ${{\mathbf}{f}}\not
  \sim {{\mathbf}{g}}$. Otherwise, the algorithm returns two sets of $n\times n$ invertible symmetric matrices $\{D,\tilde{H}_2\ldots,\tilde{H}_m\}$ and $\{D,\tilde{H}^{\prime}_2,\ldots,\tilde{H}^{\prime}_m\}$, with $D$ diagonal, defined over ${{\mathbb}{K}}$ such that: $${{\mathbf}{g}}({{\mathbf}{x}})={{\mathbf}{f}}(A\,{{\mathbf}{x}}), \text{ for } A \in \operatorname{GL}_n({{\mathbb}{K}})
  \iff
  \begin{array}{l}
    A'^{-1}\,D^{-1}\,\tilde{H}_i\,A^{\prime}=D^{-1}\,\tilde{H}_i',
    \forall \,i, 1\leq i\leq m,\\
    \text{ for } A^{\prime} \in {\mathcal}{O}_n({{\mathbb}{K}},D),
  \end{array}$$ with ${\mathcal}{O}_n({{\mathbb}{K}},D)$ denoting the set of $n \times n$ $D$-orthogonal matrices over ${{\mathbb}{K}}$.

Combining Proposition \[prop:first\_reduction\] and paragraph \[par:invertible\] any regular quadratic homogeneous instance of [`IP1S`]{}can be reduced in randomized polynomial time to “[<span style="font-variant:small-caps;">NoSolution</span>]{}”, only if the two systems are not equivalent, or to a $$(\tilde{{{\mathbf}{f}}},\tilde{{{\mathbf}{g}}})={\left(\big(\sum_{i=1}^{n}d_ix_{i}^{2},\tilde{f}_2,
    \ldots,\tilde{f}_m\big),\big(\sum_{i=1}^{n}d_ix_{i}^{2},\tilde{g}_2,\ldots,
    \tilde{g}_m\big)\right)},$$ where all the polynomials are *nondegenerate* homogeneous quadratic polynomials in ${{\mathbb}{K}}[{{\mathbf}{x}}]$. It follows that $\tilde{{{\mathbf}{f}}}\sim\tilde{{{\mathbf}{g}}} \iff \exists\, A^{\prime} \in\operatorname{GL}_n({{\mathbb}{K}})$ such that $\forall\,i, 1\leq i\leq m$, ${A^{\prime}}^{\operatorname{T}}\,\tilde{H}_i\,A^{\prime}=\tilde{H}_i'$. In particular ${A^{\prime}}^{\operatorname{T}}\, D \, A^{\prime}=D$ and $A'$ is $D$-orthogonal. Hence, ${A^{\prime}}^{\operatorname{T}}\,\tilde{H}_i\,A^{\prime}=D\,A'^{-1}\,D^{-1}\,\tilde{H}_i\,A'=
  \tilde{H}_i',\ \forall \,i, 1\leq i\leq m$.

The proof of this result implies Theorem \[th:normal\].

#### Field Extensions and Jordan Normal Form {#par:jordan}

To amplify the success probability of our results, it will be convenient to embed a field ${{\mathbb}{F}}$ in some finite extension ${{\mathbb}{F}}'$ of ${{\mathbb}{F}}$. This is motivated by the fact that matrices in ${{\mathbb}{F}}^{n\times n}$ are similar if and only if they are similar in ${{\mathbb}{F}}'^{n\times n}$, see [@Pazzis2010]. 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 ${{\mathbb}{F}}={{\mathbb}{F}}_q$ in $O(n\,{\mathsf}{M}(n)\log (q\,n))$ operations in ${{\mathbb}{F}}$, where ${\mathsf}{M}(n)$ is a bound on the number of operations in ${{\mathbb}{F}}$ to multiply two polynomials in ${{\mathbb}{F}}[x]$ of degree at most $n-1$, see [@GaGe1999 Theorem 14.14]. Then, we use [@Storjohann1998]’s algorithm to compute the generalized eigenvectors in $O(n^{\omega})$ operations in ${{\mathbb}{F}}$, with $\omega$ being the exponent of time complexity of matrix multiplication, $2\leq\omega\leq 3$.

Quadratic IP1S {#s:IP1S}
==============

In this section, we present efficient algorithms for solving regular quadratic-[`IP1S`]{}. According to Proposition \[prop:homo\], 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 $\mathcal{H}=\{D,H_2,\ldots,H_m\}$ and $\mathcal{H}'=\{D,H^{\prime}_2,\ldots,H^{\prime}_m\}$ be two families of invertible symmetric matrices in ${{\mathbb}{K}}^{n\times n}$. As explained in Theorem \[th:second\_reduction\], our task reduces – under a randomized process – to finding a $D$-orthogonal matrix $A^{\prime} \in {\mathcal}{O}_n({{\mathbb}{K}},D)$ such that: $$\label{eq:main}
  A'^{-1}\,D^{-1}\,H_i\,A^{\prime}=D^{-1}\,H_i',\  \forall \,i, 1\leq i\leq m.$$ Case $D={{\mathrm}{Id}}$ was studied in [@DBLP:conf/issac/ChistovIK97 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_i^{\operatorname{T}}\,Y=Y\,H_i^{\operatorname{T}}$. In our case, whenever $D={{\mathrm}{Id}}$, the matrices are symmetric. So, the added conditions – with the transpose – are automatically fulfilled. In [@DBLP:conf/issac/ChistovIK97], the authors suggest then to use the polar decomposition of $Y=A\,W$, with $W$ symmetric and $A$ orthogonal. Then, $A$ is an orthogonal solution of .

The main idea to compute $A$ is to compute $W$ as the square root of $Z=Y^{\operatorname{T}}\,Y$ as stated in [@DBLP:conf/issac/ChistovIK97 Section 3]. However, in general $W$ and $A$ are not defined over ${{\mathbb}{K}}$ but over ${{\mathbb}{L}}={{\mathbb}{K}}(\sqrt{\zeta_1},\ldots,\sqrt{\zeta_r})$, where $\zeta_1,\ldots,\zeta_r$ are the eigenvalues of $Z$. Assuming $\zeta_1$ is the root of an irreducible polynomial $P$ of degree $d$, then $\zeta_2,\ldots,\zeta_d$ are also roots of the same polynomial. However, there is no reason for them to be in ${{\mathbb}{K}}[x]/(P)={{\mathbb}{K}}(\zeta_1)$. But they will be the roots of a polynomial of degree $d-1$, in general, over the field ${{\mathbb}{K}}(\zeta_1)$. Then, doing another extension might only add one eigenvalue in the field. Repeating this process yields a field of degree $d!$ over ${{\mathbb}{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.

[@DBLP:conf/issac/ChistovIK97] 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=\sqrt{Y^{\operatorname{T}}\,Y}$ and then to set $A=W^{-1}\,Y$. According to Cai’s work ([@Cai1994]), some coefficients of matrix $A$ may lie in an extension of exponential degree. *Blockwise computation* (see the proof of Proposition \[prop:comput\_ortho\]) can allow us to compute such a matrix. Chistov, Ivanyos and Karpinski set $y_i$ as the restriction of $Y$ to the $i$th eigenspace, associated to $\zeta_i$, of $Y^{\operatorname{T}}\,Y$. Then, $x_i=\sqrt{\zeta_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 $i$th eigenspace of $Y^{\operatorname{T}}\,Y$ is stable by $Y$, as $W^{-1}$ would act as a multiplication by $\sqrt{\zeta_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=\sqrt{Y\,Y^{\operatorname{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 [@DBLP:conf/issac/ChistovIK97 Theorem 4] to determine an orthogonal solution. Using the analogous decomposition, the so-called Generalized Polar Decomposition ([<span style="font-variant:small-caps;">GPD</span>]{}), which depends on $D$, yields a $D$-orthogonal solution, see [@MaMaTi2005]. The [<span style="font-variant:small-caps;">GPD</span>]{}of an invertible matrix $Y$ is the factorization $Y=A\,W$, with $A$ $D$-orthogonal and $W$ in the associated Jordan algebra, *i.e.* $W^{\operatorname{T}}=D\,W\,D^{-1}$. Let us notice that $A$ and $W$ might be defined only over ${{\mathbb}{K}}'$ an algebraic extension of ${{\mathbb}{K}}$ of some degree.

\[prop:pseudoortho\_solution\] Let ${{\mathcal}{K}}=\{K_1,\ldots,K_m\}$ and ${{\mathcal}{K}}'=\{K_1',\ldots,K_m'\}$ be two subsets of $m$ matrices in ${{\mathbb}{K}}^{n\times n}$. Let $D$ be an invertible diagonal matrix. There is a $D$-orthogonal solution $A\in{{\mathbb}{K}}'^{n\times n}$ to the conjugacy problem $K_i\,A=A\,K_i'$ for all $1\leq i\leq m$, if and only if there is an invertible solution $Y\in{{\mathbb}{K}}'^{n\times n}$ to the conjugacy problem $K_i\,Y=Y\,K_i'$ and $K_i^{\operatorname{T}}\,D\,Y\,D^{-1}=D\,Y\,D^{-1}\,K_i'^{\operatorname{T}}$ for all $1\leq i\leq m$. Furthermore, if $Y=A\,W$ is the [<span style="font-variant:small-caps;">GPD</span>]{}of $Y$ with respect to $D$, then $A$ suits.

This proof is a generalization of [@DBLP:conf/issac/ChistovIK97 Section 3]. If $A$ is a $D$-orthogonal solution to the first problem, then as $A^{\operatorname{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^{\operatorname{T}}\,D\,Y$ 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^{\operatorname{T}}=D\,W\,D^{-1}$, see [@MaMaTi2005 Theorem 6.2]. As such, $W$ is a polynomial in $Z$ as proven in Section \[ss:sqrt\] 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^{\operatorname{T}}\,D\,A=W^{-\operatorname{T}}\,Y^{\operatorname{T}}\,D\,Y\,W^{-1}=
  D\,W^{-1}\,D^{-1}\,Y^{\operatorname{T}}\,D\,Y\,W^{-1}=
  D\,W^{-1}\,Z\,W^{-1}=D.\qedhere$$

For the sake of completeness, we present several efficient algorithms for performing the square root computation.

Computing the $D$-Orthogonal Solution {#ss:computation}
-------------------------------------

The goal of this part is to “$D$-orthogonalize” an invertible solution $Y \in\operatorname{GL}_n({{\mathbb}{K}})$ of equation . Instead of computing exactly $A \in {\mathcal}{O}_n({{\mathbb}{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\leq i\leq m$. To be more precise, we prove the following proposition.

\[prop:comput\_ortho\] Let ${{\mathcal}{H}}=\{H_1=D,H_2,\ldots,H_m\}$ and ${{\mathcal}{H}}'=\{H_1'=D,H_2',\ldots,H_m'\}$ be two sets of invertible matrices in ${{\mathbb}{K}}^{n\times n}$. We can compute in polynomial time two matrices representations of matrices $S$ and $T$ defined over an algebraic extension ${{\mathbb}{L}}$ such that $S\,T^{-1}$ is $D$-orthogonal and for all $1\leq i\leq 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 ${{\mathbb}{L}}$. However, matrices $S^{\operatorname{T}}\,H_i\,S$ and $T^{\operatorname{T}}\,H_i'\,T$ can be computed and tested for equality in polynomial time.

Let $Y\in\operatorname{GL}_n({{\mathbb}{K}})$ such that $D^{-1}\,H_i\,Y=Y\,D^{-1}\,H'_i,\ \forall \,i,
  1\leq i\leq m$. We set $Z=D^{-1}\,Y^{\operatorname{T}}\,D\,Y$. 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 [@Cai1994], $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 ${{\mathbb}{K}}(\zeta_i)$, where $\zeta_1,\ldots,\zeta_r$ are the eigenvalues of $Z$.

We shall say that a matrix is *block-wise* (resp. *columnblock-wise*, *rowblock-wise*) *defined over ${{\mathbb}{K}}(\zeta)$* if for all $1\leq i\leq r$, its $i$th block (resp. block of columns, block of rows) is defined over ${{\mathbb}{K}}(\zeta_i)$. The size of the $i$th block being the size of the $i$th Jordan block of $J$.

As $J=T^{-1}\,Z\,T$ is a Jordan normal form, it is block-wise defined over ${{\mathbb}{K}}(\zeta)$. Using the closed formula of Section \[ss:sqrt\], one can compute in polynomial time a square root $G$ of $J$. This matrix is a block diagonal matrix, block-wise defined over ${{\mathbb}{K}}(\sqrt{\zeta})$, 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 ${{\mathbb}{K}}(\zeta)$, see [@Cai1994 Section 4]. Since $Y$ is defined over ${{\mathbb}{K}}$, then $Y\,T$ is columnblock-wise defined over ${{\mathbb}{K}}(\zeta)$. Thus $S=A\,T=Y\,W^{-1}\,T=Y\,T\,G^{-1}$ is columnblock-wise defined over ${{\mathbb}{K}}(\sqrt{\zeta})$. We recall that product $A\,T\cdot
  T^{-1}$ mangles the eigenvalues and make each coefficient defined over ${{\mathbb}{K}}(\sqrt{\zeta_1},\ldots,\sqrt{\zeta_r})$ and thus must be avoided.

Now, to verify that $A^{\operatorname{T}}\,H\,A=H'$, for any $H\in{{\mathcal}{H}}$ and the corresponding $H'\in{{\mathcal}{H}}'$, we compute separately $S^{\operatorname{T}}\,H\,S=T^{\operatorname{T}}\,A^{\operatorname{T}}\,H\,A\,T$ and $T^{\operatorname{T}}\,H'\,T$. For the former, $S=A\,T$ (resp. $S^{\operatorname{T}}=(A\,T)^{\operatorname{T}}$) is columnblock-wise (resp. rowblock-wise) defined over ${{\mathbb}{K}}(\sqrt{\zeta})$ and $H$ is defined over ${{\mathbb}{K}}$. Therefore, the product matrix makes each of the coefficients which are on both the $i$th block of rows and the $j$th block of columns defined over ${{\mathbb}{K}}(\sqrt{\zeta_i},\sqrt{\zeta_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 ${{\mathbb}{K}}(\zeta)$.

Let us assume that the characteristic polynomial of $Z$, of degree $n$, can be factored as $P_1^{e_1}\cdots P_s^{e_s}$ with $P_i$ and $P_j$ coprime whenever $i\neq j$, $\deg P_i=d_i$ and $e_i\geq
  1$. From a computation point of view, one needs to introduce a variable $\alpha_{i,j}$ for each root of $P_i$ and then a variable $\beta_{i,j}$ for the square root of $\alpha_{i,j}$. This yields a total number of $2\sum_{i=1}^sd_i$ variables. In Section \[ss:sqrt\_char\_p\], we present another method which manages to introduce only $2s$ variables in characteristic $p>2$.

Probabilistic and Deterministic Algorithms {#ss:proba_algo}
------------------------------------------

We first describe a simple probabilistic algorithm summarizing the method of Section \[ss:computation\].

\[al:one\_sol\] Probabilistic algorithm.

Input

:   Two sets of invertible symmetric matrices ${{\mathcal}{H}}=\{H_1=D,\ldots,H_m\}\subseteq{{\mathbb}{K}}^{n\times n}$ and ${{\mathcal}{H}}'=\{H_1'=D,\ldots,H_m'\}\subseteq{{\mathbb}{K}}^{n\times
          n}$.

Output

:   A description of the matrix $A \in \operatorname{GL}_n({{\mathbb}{L}})$ such that $H_i'=A^{\operatorname{T}}\,H_i\,A$ for all $1\leq i\leq m$ or “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

1.  Compute the vector subspace ${\mathcal}{Y}=\{Y\ \mid
        \ D^{-1}\,H_i\,Y=Y\,D^{-1}\,H_i',\
        \forall\,1\leq i\leq m\}\subseteq{{\mathbb}{K}}^{n \times n}$.

2.  **If** ${\mathcal}{Y}$ is reduced to the null matrix **then return** “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

3.  Pick at random $Y\in{\mathcal}{Y}$.

4.  Compute $Z=D^{-1}\,Y^{\operatorname{T}}\,D\,Y$ and $J=T^{-1}Z\,T\in{{\mathbb}{L}}^{n\times 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$.

\[proba:algo\] Algorithm \[al:one\_sol\] is correct with probability at least $1-n/|{{\mathbb}{K}}|$ and runs in polynomial time.

The correctness and the polynomial-time complexity of the algorithm come from Section \[ss:computation\]. After computing ${\mathcal}{Y}$ and putting the equations defining its matrices in triangular form, one has to pick at random one matrix $Y
  \in {\mathcal}{Y}$. By sampling the whole field ${{\mathbb}{K}}$ on these free variables, the probability that $\det Y=0$ is upper bounded by $n/|{{\mathbb}{K}}|$ thanks to Schwartz-Zippel-DeMillo-Lipton Lemma ([@DL78; @Z79]).

Let us recall that the conjugacy problem does not depend on the ground field (see [@Pazzis2010]), *i.e.* if there exists $Y\in\operatorname{GL}_n({{\mathbb}{K}}')$, such that $H_i\,Y=Y\,H_i'$, then there exists $Y'\in\operatorname{GL}_n({{\mathbb}{K}})$ such that $H_iY'=Y'H_i'$. This allows us to extend ${{\mathbb}{K}}$ to a finite extension in order to decrease the probability of getting a singular matrix $Y$. Thus the success probability of Algorithm \[al:one\_sol\] can be amplified to $1-n/|{{\mathbb}{K}}'|$ for any extension ${{\mathbb}{K}}'\supseteq{{\mathbb}{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 ${{\mathbb}{K}}'(\sqrt{\zeta_1},\ldots,\sqrt{\zeta_r})$. Notice also that this algorithm can be turned into a deterministic algorithm using [@DBLP:conf/issac/ChistovIK97 Theorem 2]. That is, there is a polynomial-time algorithm allowing to compute an invertible element in ${\mathcal}{Y}$. Furthermore, if one of the original Hessian matrices is already invertible, the computations of the essential variables of paragraph \[par:variables\] and the search of an equation with $n$ essential variables in paragraph \[par:identity\] can be done in a deterministic way. Whence, the whole algorithm is deterministic.

The main Theorem \[th:main\] summarizes this remark together with Theorem \[proba:algo\].

The binary Case {#ss:char_2}
---------------

In this section, we investigate fields of characteristic $2$. Let ${{\mathbb}{K}}={{\mathbb}{F}}_q$ and $({{\mathbf}{f}},{{\mathbf}{g}})\in {{\mathbb}{K}}[{{\mathbf}{x}}]^m \times {{\mathbb}{K}}[{{\mathbf}{x}}]^m$. Instead of Hessian matrices, we consider equivalently upper triangular matrices $H_1,\ldots,H_m$ and $H_1',\ldots,H_m'$ such that: $$f_i({{\mathbf}{x}})={{\mathbf}{x}}^{\operatorname{T}}\,H_i\,{{\mathbf}{x}},\ g_i({{\mathbf}{x}})={{\mathbf}{x}}^{\operatorname{T}}\,H_i'\,{{\mathbf}{x}},\quad\forall\,1\leq
i\leq m.$$

For any matrix $M\in{{\mathbb}{K}}^{n\times n}$, let us denote $\Delta(M)=\operatorname{Diag}(m_{11},\ldots,m_{nn})$ and $\Sigma(M)=M+M^{\operatorname{T}}$. It is classical that if there exists $A \in \operatorname{GL}_n({{\mathbb}{K}})$ such that ${{\mathbf}{g}}({{\mathbf}{x}})={{\mathbf}{f}}(A\cdot{{\mathbf}{x}})$, then we also have $$\begin{aligned}
  \Sigma(H_i')&=A^{\operatorname{T}}\,\Sigma(H_i)\,A,
  \label{eq:char2sigma}\\
  \Delta(H_i')&=\Delta(A^{\operatorname{T}}\,H_i\,A),\quad\forall\,i,\ 1\leq i\leq m.
  \label{eq:char2delta}\end{aligned}$$ It suffices for this to expand ${{\mathbf}{f}}(A\cdot{{\mathbf}{x}})$ and to consider the upper triangular matrices. In a sense, $\Sigma(H_i)$ is the Hessian matrix of $f_i$ and $\Delta(H_i)$ allow us to remember the $x_j^2$ terms in $f_i$. Combining two equations of  yields $\Sigma(H_j')^{-1}\Sigma(H_i')=A^{-1}\Sigma(H_j)^{-1}\Sigma(H_i)A$ as long as $\Sigma(H_j)$ is invertible. Let us notice that $\Sigma(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 $\sum_{\sigma\in{\mathfrak}{S}_n}\prod_{i=1}^ns_{i,\sigma(i)}$ yields, for each nonzero term $\prod_{i=1}^ns_{i,\sigma(i)}$, the term $\prod_{i=1}^ns_{i,\sigma^{-1}(i)}=\prod_{i=1}^ns_{\sigma(i),i}=
\prod_{i=1}^ns_{i,\sigma(i)}$. Dimension $n$ being odd, $\prod_{i=1}^ns_{i,\sigma(i)}$ cannot be the same term as $\prod_{i=1}^ns_{i,\sigma^{-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 ${{\mathbb}{Z}}_q$ the unramified extension of the ring of dyadic integers of degree $\log_2 q$. Let $\tilde{S}\in{{\mathbb}{Z}}_q^{n\times n}$ be antisymmetric such that $\tilde{S}\mapsto S$. Then $\det\tilde{S}=\det \tilde{S}^{\operatorname{T}}=\det(-\tilde{S})=(-1)^n\det\tilde{S}$, hence $\det\tilde{S}=0$ and $\det S=0$.

Therefore, if $n$ is odd, then a linear combination of the $\Sigma(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 \[as:regular\] to assume that a linear combination $\sum_{i=1}^m\lambda_i\,f_i$ is not degenerate, and $\lambda_1,\ldots,\lambda_m$ can be found in randomized polynomial time. Assuming $\lambda_1\neq 0$, we substitute the linear combinations $\sum_{i=1}^m\lambda_i\,H_i$ and $\sum_{i=1}^m\lambda_i\,H_i'$ to $H_1$ and $H_1'$.

As a consequence, we can find linear forms $\ell_1,\ldots,\ell_n$ in ${{\mathbf}{x}}$ such that, see [@finite_fields Theorem 6.30]: $f_1({{\mathbf}{x}})=\ell_1\ell_2+\ell_3\ell_4+\cdots+\ell_{n-1}\ell_n$ or $f_1({{\mathbf}{x}})=\ell_1\ell_2+\ell_3\ell_4+\cdots+\ell_{n-1}\ell_n+\ell_{n-1}^2
+d\ell_n^2$, where $\operatorname{Tr}_{{{\mathbb}{K}}}(d)=d+d^2+\cdots+d^{q/2}=1$. After applying this change of variables, $\Sigma(H_1)$ is always the following invertible block diagonal matrix: $$\Sigma(H_1)=\operatorname{Diag}\left(\begin{pmatrix}0 &1\\1 &0\end{pmatrix},
  \ldots,\begin{pmatrix}0 &1\\1 &0\end{pmatrix}\right).$$ Following paragraph \[par:invertible\], we once again choose $\nu_i$ such that $\Sigma(H_i+\nu_i\,H_1)$ is invertible and replace $H_i$ by $H_i+\nu_i\,H_1$. Thus, Proposition \[prop:first\_reduction\] and Theorem \[th:second\_reduction\] become:

\[prop:char\_2\] Let $n$ be an even integer and ${{\mathbb}{K}}$ be a field of characteristic $2$. Let $({{\mathbf}{f}},{{\mathbf}{g}})\in{{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m $ be regular quadratic homogeneous polynomials. There is a randomized polynomial-time algorithm which returns “[<span style="font-variant:small-caps;">NoSolution</span>]{}” only if ${{\mathbf}{f}}\not\sim{{\mathbf}{g}}$ or a new instance $$(\tilde{{{\mathbf}{f}}},\tilde{{{\mathbf}{g}}})=((\delta,\tilde{f}_2,\ldots,\tilde{f}_m),
  (\delta,\tilde{g}_2,\ldots,\tilde{g}_m)){{\mathbb}{K}}[x_1,\ldots,x_n]^m\times
  {{\mathbb}{K}}[x_1,\ldots,x_n]^m$$ such that ${{\mathbf}{f}}\sim{{\mathbf}{g}}\iff\tilde{{{\mathbf}{f}}}\sim\tilde{{{\mathbf}{g}}}$. Furthermore, denoting $D$ the upper triangular matrix of $\tilde{f}_1=\tilde{g}_1=\delta$, [`IP1S`]{}comes down to a $\Sigma(D)$-Orthogonal Simultaneous Matrix Conjugacy problem, *i.e.* conjugacy by an $\Sigma(D)$-orthogonal matrix under some constraints: $$\begin{aligned}
    A^{\operatorname{T}}\,\Sigma(D)\,A&=\Sigma(D)\text{ and }\forall\,i,2\leq i\leq m,
    \Sigma(D)^{-1}\,\Sigma(H_i')=A^{-1}\,\Sigma(D)^{-1}\,\Sigma(H_i)\,A,\\
    \Delta(A^{\operatorname{T}}\,H_i\,A)&=\Delta(H_i').
  \end{aligned}$$

We mimic the proof of Proposition \[prop:first\_reduction\]. We compute in randomized polynomial time $\lambda_1,\ldots,\lambda_m \in {{\mathbb}{K}}$ such that ${\varphi}=\sum_{i=1}^{m}\lambda_i \cdot f_i$ is regular and we define $\gamma=\sum_{i=1}^{m}\lambda_i \cdot g_i$. Assuming w.lo.g. $\lambda_1\neq 0$. We have then: $${{\mathbf}{f}}\sim {{\mathbf}{g}}\iff ({\varphi},f_2,\ldots,f_m) \sim (\gamma,g_2,\ldots,g_m).$$ Computing $\delta$ the canonical quadratic form equivalent to ${\varphi}$ yields a $P \in\operatorname{GL}_n({{\mathbb}{K}})$ such that $\tilde{{{\mathbf}{f}}}=(\tilde{{\varphi}}=\delta,\tilde{f_2},\ldots,
  \tilde{f_m})=({\varphi}(P{{\mathbf}{x}}),f_2(P{{\mathbf}{x}}),\ldots,f_m(P{{\mathbf}{x}}))$.

Then computing the canonical quadratic form equivalent to $\gamma$ allows us to compare it with $\delta$. If they are different, then ${{\mathbf}{f}}\not \sim {{\mathbf}{g}}$ and we return “[<span style="font-variant:small-caps;">NoSolution</span>]{}”. Otherwise, the reduction is given by a matrix $Q \in
  \operatorname{GL}_n({{\mathbb}{K}})$ such that $\tilde{{{\mathbf}{g}}}=(\tilde{\gamma}=\delta,\tilde{g}_2,\ldots\tilde{g}_m)=
  (\gamma(Q{{\mathbf}{x}}),g_2(Q{{\mathbf}{x}}),\ldots,g_m(Q{{\mathbf}{x}}))$. Thus, $\tilde{{{\mathbf}{f}}}\sim\tilde{{{\mathbf}{g}}}$ if and only if ${{\mathbf}{f}}\sim{{\mathbf}{g}}$.

Finally, equations  $A^{\operatorname{T}}\,\Sigma(H_i)\,A=\Sigma(H_i')$ for all $i$, $1\leq i\leq m$ of can be rewritten as $A^{\operatorname{T}}\,\Sigma(D)\,A=\Sigma(D)$ and $\Sigma(D)^{-1}\,\Sigma(H_i')=A^{-1}\,\Sigma(D)^{-1}\,\Sigma(H_i)A$ for all $i$, $2\leq i\leq m$, while equations  $\Delta(A^{\operatorname{T}}\,H_i\,A)=\Delta(H_i')$ for all $i$, $1\leq i\leq 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$.

\[al:char\_2\] Probabilistic algorithm in characteristic $2$.

Input

:   Two sets of triangular matrices ${{\mathcal}{H}}=\{H_1=D,\ldots,H_m\}\subseteq{{\mathbb}{K}}^{n\times n}$ and ${{\mathcal}{H}}'=\{H_1'=D,\ldots,H_m'\}\subseteq{{\mathbb}{K}}^{n\times n}$ such that $H_1+H_1^{\operatorname{T}}$ and $H_1'+H_1'^{\operatorname{T}}$ are invertible.

Output

:   A description of the matrix $A \in \operatorname{GL}_n({{\mathbb}{L}})$ such that $H_i'=A^{\operatorname{T}}\,H_i\,A$ for all $1\leq i\leq m$ or “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

1.  Compute the vector subspace ${\mathcal}{Y}=\{Y\ \mid
        \ \Sigma(D)^{-1}\,\Sigma(H_i)\,Y=Y\,\Sigma(D)^{-1}\,\Sigma(H_i'),\
        \forall\,1\leq i\leq m\}\subseteq{{\mathbb}{K}}^{n \times n}$.

2.  **If** ${\mathcal}{Y}$ is reduced to the null matrix **then return** “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

3.  Pick at random $Y\in{\mathcal}{Y}$.

4.  Compute $Z=\Sigma(D)^{-1}\,\Sigma(Y)^{\operatorname{T}}\,\Sigma(D)\,\Sigma(Y)$ and $J=T^{-1}Z\,T\in{{\mathbb}{L}}^{n\times n}$, the Jordan normal form of $Z$ together with $T$.

5.  **While** $J$ is not diagonal

    1.  Pick at random $Y\in{\mathcal}{Y}$.

    2.  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 \[ss:sqrt\_char\_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 {#open-question-the-irregular-case-1 .unnumbered}

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 {#ss:bench}
----------

We present in this section some timings of our algorithms over instances of [`IP1S`]{}. We created instances ${\mathcal}{H}=\{H_1,\ldots,H_m\}$ and ${\mathcal}{H}'=\{H_1',\ldots,H_m'\}$ which are randomly alternatively equivalent over ${{\mathbb}{F}}_p$, equivalent over ${{\mathbb}{F}}_{p^2}$ but not ${{\mathbb}{F}}_p$ or not equivalent at all over $\bar{{{\mathbb}{F}}}_p$, the algebraic closure of ${{\mathbb}{F}}_p$, for an odd $p$. We report our timings in the following Table \[tab:char\_odd\] obtained using one core of an <span style="font-variant:small-caps;">Intel Core i7</span> at $2.6 \si{\giga\hertz}$ running <span style="font-variant:small-caps;">Magma 2.19</span>, [@magma], on <span style="font-variant:small-caps;">Linux</span> with $16\si{\giga}{\mathrm}{B}$ of RAM. These timings corresponds to solving the linear system which is the dominant part in our algorithm with complexity $O(n^{2\,\omega})$. 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\leq i\leq
m$. In fact, in practice, $i=2,3$ give enough equations to retrieve $A$ up to one free parameter if ${\mathcal}{H}$ and ${\mathcal}{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^{\operatorname{T}}\,H_1\,A=H_1'$. Let us notice that either all these equations can be solved over ${{\mathbb}{F}}_p$ or none of them can. If they can, then ${\mathcal}{H}$ and ${\mathcal}{H}'$ are equivalent over ${{\mathbb}{F}}_p$ and we have determined $A$ up to a sign, otherwise ${\mathcal}{H}$ and ${\mathcal}{H}'$ are only equivalent over ${{\mathbb}{F}}_{p^2}$ but not ${{\mathbb}{F}}_p$ and we also have computed such an $A$. This yields Algorithm \[al:crypto\].

\[al:crypto\] Simplified Algorithm.

Input

:   Two sets of generic invertible symmetric matrices ${{\mathcal}{H}}=\{H_1,\ldots,H_m\}\subseteq{{\mathbb}{K}}^{n\times n}$ and ${{\mathcal}{H}}'=\{H_1',\ldots,H_m'\}\subseteq{{\mathbb}{K}}^{n\times n}$.

Output

:   A matrix $A \in \operatorname{GL}_n({{\mathbb}{K}})$ such that $H_i'=A^{\operatorname{T}}\,H_i\,A$ for all $1\leq i\leq m$ or “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

1.  Compute the vector subspace ${\mathcal}{Y}=\{Y\ \mid
        \ H_1^{-1}\,H_i\,Y=Y\,{H_1'}^{-1}\,H_i',\
        \forall\,2\leq i\leq 3\}\subseteq{{\mathbb}{K}}^{n \times n}$.

2.  **If** ${\mathcal}{Y}$ is reduced to a space of singular matrices **then return** “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

3.  Determine $Y_0$ such that ${\mathcal}{Y}={\left\{\lambda Y_0 \mid \lambda\in{{\mathbb}{K}}\right\}}$.

4.  Solve in $\lambda$ one equation $\lambda^2(Y_0^{\operatorname{T}}\,H_1\,Y_0)_{i,j}=(H_1')_{i,j}$ for a suitable pair $(i,j)$.

5.  Set $A=\lambda\,Y_0$.

6.  Pick at random $r\in{{\mathbb}{K}}^n$.

7.  Check that $A^{\operatorname{T}}\,H_i\,A\,r=H_i'\,r$ for all $i$, $1\leq i\leq 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\,\omega})$ operations in ${{\mathbb}{K}}$. Then, one needs to determine $\lambda$ by extracting one square root in ${{\mathbb}{K}}$ which can be done in $O((\log q)^3)$ operations in ${{\mathbb}{K}}={{\mathbb}{F}}_q$ with Tonelli–Shanks’s algorithm, [@Shanks1973]. Finally, one can check that $A^{\operatorname{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 ${{\mathbb}{K}}$.

Recall that, in practice, the best matrix multiplication algorithm is due to [@Strassen1969] whose complexity is in $O(n^{\log_2 7})\subseteq 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$.

     $n$         $20$     $30$     $40$    $50$    $60$   $70$   $80$   $90$   $100$ 
  --------- --------- -------- -------- ------- ------- ------ ------ ------ ------- --
   Timings    $0.040$   $0.20$   $0.84$   $2.7$   $7.5$   $17$   $40$   $79$   $130$ 

  : Timings for solving [`IP1S`]{}over ${{\mathbb}{F}}_{65521}$ in $\si{s}$.[]{data-label="tab:char_odd"}

In Table \[tab:char\_2\], we report our timings for solving the linear system of our algorithm in characteristic $2$ presented in Section \[ss:char\_2\]. Our method does not differ much from the one in odd characteristic. We pick at random two sets of $m$ upper triangular matrices over ${{\mathbb}{F}}_2$ which are either equivalent over ${{\mathbb}{F}}_2$ or not equivalent at all over $\bar{{{\mathbb}{F}}}_2$, the algebraic closure of ${{\mathbb}{F}}_2$. We first solve the linear system $\Sigma(H_1)^{-1}\,\Sigma(H_i)\,A=A\,\Sigma(H_1')^{-1}\,\Sigma(H_i')$, for all $i$, $2\leq i\leq 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 ${{\mathcal}{H}}$ and ${{\mathcal}{H}}'$ are equivalent, the linear system yields in general $A$ up to one free parameter if $m\geq 3$ (resp. $m\geq 5$). Otherwise, it yields the zero matrix. Then, it suffices to solve one of the quadratic equations amongst the one given by $A^{\operatorname{T}}\,\Sigma(H_1)\,A=\Sigma(H_1')$ and $\Delta(A^{\operatorname{T}}\,H_i\,A)=\Delta(H_i')$, for all $i$, $1\leq i\leq m$ (see Proposition \[prop:char\_2\]).

We compare the timings of both <span style="font-variant:small-caps;">MAGMA</span> and the C library <span style="font-variant:small-caps;">M4RI</span>, due to [@m4ri].

Once again, our complexity in $O(n^{2\,\omega})$ 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 ${{\mathbb}{F}}_{65521}$ although it would seem a lot of quadratic equations must be solved.

                                $n$                                    $20$      $30$      $40$                $50$     $60$               $70$               $80$     $90$               $100$ 
  --------------------------------------------------------------- --------- --------- --------- ------------------- -------- ------------------ ------------------ -------- ------------------- --
   Timings (<span style="font-variant:small-caps;">MAGMA</span>)    $0.010$   $0.030$   $0.080$   $0.25\phantom{0}$   $0.68$   $1.4\phantom{0}$   $3.2\phantom{0}$   $6.25$   $16\phantom{.00}$ 
   Timings (<span style="font-variant:small-caps;">M4RI</span>)                         $0.010$             $0.030$   $0.06$             $0.14$             $0.27$   $0.51$              $0.91$ 

  : Timings for solving [`IP1S`]{}over ${{\mathbb}{F}}_2$ in $\si{s}$.[]{data-label="tab:char_2"}

Counting the Solutions: \#[`IP1S`]{} {#s:enumeration}
====================================

In this part, we present a method for counting the number of solutions to quadratic-[`IP1S`]{}. The main result is a consequence of [@Singla2010 Lemma 4.11]. According to Proposition \[prop:homo\], 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 $({{\mathbf}{f}},{{\mathbf}{g}}) \in {{\mathbb}{K}}[{{\mathbf}{x}}]^m \times {{\mathbb}{K}}[{{\mathbf}{x}}]^m$.

Let ${{\mathcal}{H}}=\{H_1=D,\ldots,H_m\}$ and ${{\mathcal}{H}}^{\prime}=\{H^{\prime}_1=D,\ldots,H^{\prime}_m\}$ be the Hessian matrices in ${{\mathbb}{K}}^{n\times n}$ of ${{\mathbf}{f}}$ and ${{\mathbf}{g}}$ respectively. Our counting problem is equivalent to enumerating the number of $D$-orthogonal matrices $X$ satisfying: $$\begin{aligned}
\label{eq:enum}
  X^{-1}\,D^{-1}\,H_i \, X=D^{-1}\,H^{\prime}_i,  \quad \forall\,i, 1 \leq i
  \leq m.\end{aligned}$$ In [@Singla2010 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 \[ap:enum\_proof\].

Let us notice that if $X$ and $X'$ are both orthogonal solutions of , then $XX'^{-1}$ commutes with $D^{-1}{{\mathcal}{H}}$ (resp. $X^{-1}X'$ commutes with $D^{-1}{{\mathcal}{H}}^{\prime}$). Therefore, the size of the set of solutions is upper bounded by the number of invertible elements in the centralizer $\operatorname{{\mathcal}{C}}(D^{-1}{{\mathcal}{H}})$ of $D^{-1}{{\mathcal}{H}}$.

Let $\alpha$ be an algebraic element of degree $m$ over ${{\mathbb}{K}}$ and let ${{\mathbb}{K}}'={{\mathbb}{K}}(\alpha)$. We consider the matrix $H=D^{-1}{\left(H_1+\cdots+\alpha^{m-1}H_m\right)}\in{{\mathbb}{K}}'^{n\times n}$. It is clear that a matrix $X\in{{\mathbb}{K}}^{n\times n}$ is such that $X^{-1}\,D^{-1}\,H_i\,X=D^{-1}\,H_i$ for all $i,1\leq i\leq m$ if and only if $X^{-1}HX=H$. Hence, the problem again reduces itself to the computation of the centralizer $\operatorname{{\mathcal}{C}}(H)$ of $H$ intersected with $\operatorname{GL}_n({{\mathbb}{K}})$. To ease the analysis, we consider the subspace ${\mathcal}{V}=\operatorname{{\mathcal}{C}}(H)\cap{{\mathbb}{K}}^{n\times n}$ of matrices in ${{\mathbb}{K}}^{n\times n}$ commuting with $H$. This provides an upper bound on the number of solutions. The dimension of ${\mathcal}{V}$ as a ${{\mathbb}{K}}$-vector space is upper bounded by the dimension of $\operatorname{{\mathcal}{C}}(H)$ as a ${{\mathbb}{K}}'$-vector space. Indeed, ${\mathcal}{V}\otimes{{\mathbb}{K}}'\subseteq\operatorname{{\mathcal}{C}}(H)$, hence $\dim_{{{\mathbb}{K}}}{\mathcal}{V}=\dim_{{{\mathbb}{K}}'}({\mathcal}{V}\otimes{{\mathbb}{K}}')\leq
\dim_{{{\mathbb}{K}}'}\operatorname{{\mathcal}{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 ${{\mathbb}{L}}$.

Let us denote $\zeta_1,\ldots,\zeta_r$ the eigenvalues of $J$. According to [@Singla2010 Lemma 4.11] and Lemma \[lm:dim\_cent\] in \[ap:enum\_proof\], if $J$ is made of Jordan blocks of size $s_{i,1}\leq\cdots\leq s_{i,d_i}$ for $i,1\leq i\leq
r$, then centralizer of $H$ has dimension at most $$\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i}(2d_i-2j+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 ${{\mathbb}{F}}_q^{n\times n}$.

\[cor:enum\] Let $H_1,\ldots,H_m\in{{\mathbb}{F}}_q^{n\times n}$ be symmetric matrices. Let $\alpha$ be algebraic over ${{\mathbb}{F}}_q$ of degree $m$. Let $H=\sum_{i=1}^m\alpha^{i-1}\,H_i\in{{\mathbb}{F}}_{q^m}^{n\times n}$ and let $J$ be its normal Jordan form with eigenvalues $\zeta_1,\ldots,\zeta_r$. Assuming the blocks of $J$ associated to $\zeta_i$ are $J_{\zeta_i,s_{i,1}},\ldots,J_{\zeta_i,s_{i,d_1}}$ with $s_{i,1}\leq\cdots\leq s_{i,d_i}$ for $i,1\leq i\leq r$, then the number of solutions of quadratic-[`IP1S`]{}in ${{\mathbb}{F}}_q^{n\times n}$ on the instance $(H_1,\ldots,H_m)$ is at most $$q^{{\left(\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i}(2\,d_i-2\,j+1)\,s_{i,j}\right)}}-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 [@DFPW11; @WP11]). In particular, in [@DFPW11], 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 [@DFPW11] with our approach.

Special Case of the general IP Problem {#s:ip}
======================================

In this part, we present a randomized polynomial-time algorithm for the following task:\
**Input:** ${{\mathbf}{g}}=(g_1,\ldots,g_n) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$, and ${{\mathrm}{{\mathbf}{POW}}}_{n,d}=(x_1^d,\ldots,x_n^d) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$ for some $d>0$.\
**Question:** Find – if any – $(A,B) \in\operatorname{GL}_n({{\mathbb}{K}}) \times
\operatorname{GL}_n({{\mathbb}{K}})$ such that: $${{\mathbf}{g}}=B \cdot {{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\cdot {{\mathbf}{x}}),
\mbox{ with ${{\mathbf}{x}}=(x_1,\ldots,x_n)^{\operatorname{T}}$.}$$ In [@Kayal2011], 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({{\mathbf}{x}})=f(A\cdot{{\mathbf}{x}})$ with $f({{\mathbf}{x}})=\sum_{i=1}^nx_i^d$.

We generalize this result to $m=n$ with an additional transformation on the polynomials. The main tool of our method is the following theorem.

\[th:Kaltofen\] Let ${{\mathbf}{g}}=(g_1,\ldots,g_n)$ be polynomials of degree $d$ over ${{\mathbb}{K}}[x_1,\ldots,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 ${{\mathbf}{g}}$. If the size of ${{\mathbb}{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 ${{\mathbf}{g}}$ in randomized polynomial time.

For this, we will use [@Kaltofen89factorizationof]’s algorithm to factor a polynomial given by evaluation, the needed size of ${{\mathbb}{K}}$ is a consequence of this. As ${{\mathbf}{g}}$ has at most $n\,\binom{n+d-1}{d}\in O(n^{d+1})$ monomials, it can be evaluated in polynomial time using a multivariate Horner’s scheme. Each $\frac{\partial g_i}{\partial x_j}({\mathbf}{a})$ is obtained as the coefficient in front of $x_j$ of the expansion of $g_i(a_1,\ldots,a_{j-1},a_j+x_j,a_{j+1},\ldots,a_n)$ which is a univariate polynomial of degree at most $d$. By [@BiniP1994 Chapter 1, Section 8], this can be computed as the shift of a polynomial in polynomial time. Hence the Jacobian matrix of ${{\mathbf}{g}}$ at ${\mathbf}{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 ${{\mathbf}{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^{\omega})$ operations, with $\omega$ being the exponent of time complexity of matrix multiplication, $2\leq\omega\leq 3$. Using the arithmetic circuit of polynomial size $L$ to evaluate the determinant of the Jacobian matrix, one can use <span style="font-variant:small-caps;">Kaltofen</span>’s algorithm to factor it in polynomial time.

As in [@Kayal2011] or [@DBLP:conf/eurocrypt/Perret05], we use partial derivatives to extract matrices $A$ and $B$. The idea is to factor the Jacobian matrix (whereas [@Kayal2011] uses the Hessian matrix) of ${{\mathbf}{g}}$ at ${{\mathbf}{x}}$ which is defined as follows: $${{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}}) = {\left(\partial_j g_i=\frac{\partial g_i} {\partial
    x_j}\right)}_{\substack{1\leq i\leq m\\1\leq j\leq n}}.$$ According to the following lemma, the Jacobian matrix is especially useful in our context:

\[lm:Jac\] Let $\big( {{\mathbf}{f}}=(f_1,\ldots,f_m),{{\mathbf}{g}}=(g_1,\ldots,g_m)\big) \in
  {{\mathbb}{K}}[x_1,\ldots,x_n]^m \times {{\mathbb}{K}}[x_1,\ldots,x_n]^m$. If $(A,B) \in\operatorname{GL}_n({{\mathbb}{K}}) \times \operatorname{GL}_m({{\mathbb}{K}})$ are such that ${{\mathbf}{g}}=B \cdot
  {{\mathbf}{f}}(A\cdot {{\mathbf}{x}})$, then $${{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}}) = B \cdot {{\mathrm}{J}}_{{{\mathbf}{f}}}( A\cdot {{\mathbf}{x}}) \cdot A.$$ As a consequence, $\det{{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}})=\det A \cdot \det B
  \cdot \det{{\mathrm}{J}}_{{{\mathbf}{f}}}( A\cdot {{\mathbf}{x}})$.

As long as $\operatorname{char}{{\mathbb}{K}}$ does not divide $d$, the Jacobian matrix of ${{\mathbf}{f}}={{\mathrm}{{\mathbf}{POW}}}_{n,d}({{\mathbf}{x}})$ is an invertible diagonal matrix whose diagonal elements are $\big({{\mathrm}{J}}_{{{\mathbf}{f}}}({{\mathbf}{x}}) \big)_{i,i}= d \cdot
x_i^{d-1}$, $\forall\,i, 1 \leq i \leq n$. Thus: $$\det{{\mathrm}{J}}_{{{\mathrm}{{\mathbf}{POW}}}_{n,d}}({{\mathbf}{x}})=d^n \prod_{i=1}^{n}x_i^{d-1}.$$ This gives

\[lm:det\] Let ${{\mathbf}{g}}=(g_1,\ldots,g_n) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$. Let $d>0$ be an integer, and define ${{\mathrm}{{\mathbf}{POW}}}_{n,d}=(x_1^d,\ldots,x_n^d) \in
  {{\mathbb}{K}}[x_1,\ldots,x_n]^m$. If $(A,B) \in\operatorname{GL}_n({{\mathbb}{K}}) \times
  \operatorname{GL}_n({{\mathbb}{K}})$ are such that ${{\mathbf}{g}}({{\mathbf}{x}})=B \cdot {{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\cdot {{\mathbf}{x}})$, then: $$\det{{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}}) = c \cdot \prod_{i=1}^{n} \ell_i({{\mathbf}{x}})^{d-1},$$ with $c \in {{\mathbb}{K}}\setminus \{ 0\}$, and the $\ell_i$’s are linear forms whose coefficients are the $i$th rows of $A$.

According to Lemma \[lm:Jac\], $\det\big( {{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}})
  \big)=\det(A) \cdot \det(B) \cdot d^{n} \cdot \prod_{i=1}^{n}
  \ell_i({{\mathbf}{x}})^{d-1}$.

From Lemmas \[lm:det\], we can derive a randomized polynomial-time algorithm for solving [`IP`]{}on the instance $({{\mathbf}{f}}={{\mathrm}{{\mathbf}{POW}}}_{n,d},{{\mathbf}{g}}) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n
\times {{\mathbb}{K}}[x_1,\ldots,x_n]^n$ in characteristic $0$. It suffices to use [@Kaltofen89factorizationof]’s algorithm for factoring $\det{{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{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 \[th:ip\] whenever $\operatorname{char}{{\mathbb}{K}}\nmid d$.

#### Small characteristic

If $\operatorname{char}{{\mathbb}{K}}$ divides $d$, we must change a little bit our strategy. Let us write $d=p^re$ with $\operatorname{char}{{\mathbb}{K}}=p$ and $e$ coprime. Then, $$\begin{aligned}
  {{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\,{{\mathbf}{x}}) &=
  {\left({\left(\sum_{j=1}^na_{1,j}\,x_j\right)}^{p^re},\ldots,
    {\left(\sum_{j=1}^na_{n,j}\,x_j\right)}^{p^re}\right)}\\
  {{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\,{{\mathbf}{x}}) &={\left({\left(\sum_{j=1}^na_{1,j}^{p^r}\,x_j^{p^r}\right)}^e,\ldots,
    {\left(\sum_{j=1}^na_{n,j}^{p^r}\,x_j^{p^r}\right)}^e\right)}\\
  {{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\,{{\mathbf}{x}}) &={{\mathrm}{{\mathbf}{POW}}}_{n,e}{\left(A^{(p^r)}\,{{\mathbf}{x}}^{p^r}\right)},\end{aligned}$$ with $A^{(p^r)}={\left(a_{i,j}^{p^r}\right)}_{1\leq i,j\leq n}$ and ${{\mathbf}{x}}^{p^r}={\left(x_1^{p^r},\ldots,x_n^{p^r}\right)}$. Thus ${{\mathbf}{g}}$ is a polynomial in ${{\mathbf}{x}}^{p^r}$ and by replacing ${{\mathbf}{x}}^{p^r}$ by ${{\mathbf}{x}}$, the problem comes down to checking if $\tilde{{{\mathbf}{g}}}=B\cdot{{\mathrm}{{\mathbf}{POW}}}_{n,e}(A^{(p^r)}\cdot{{\mathbf}{x}})$ where ${{\mathbf}{g}}({{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}}^{p^r})$. Now, as $\tilde{{{\mathbf}{g}}}({{\mathbf}{x}})=B\cdot{{\mathrm}{{\mathbf}{POW}}}_{n,e}{\left(A^{(p^r)}\,{{\mathbf}{x}}\right)}$, then $${{\mathrm}{J}}_{\tilde{{{\mathbf}{g}}}}=B\cdot{{\mathrm}{J}}_{{{\mathrm}{{\mathbf}{POW}}}_{n,e}}(A^{(p^r)}\,{{\mathbf}{x}})\cdot A^{(p^r)}.$$ Hence, $\det{{\mathrm}{J}}_{\tilde{{{\mathbf}{g}}}}({{\mathbf}{x}})=\det B\,\det{{\mathrm}{J}}_{{{\mathrm}{{\mathbf}{POW}}}_{n,e}}(A^{(p^r)}\,{{\mathbf}{x}})
\det A^{(p^r)}=e^n\prod_{i=1}^n
{\tilde{\ell}_i({{\mathbf}{x}})}^{e-1}(\det A)^{p^r}\det B$, where the $\tilde{\ell}_i$’s are linear forms whose coefficients are the $i$th rows of $A^{(p^r)}$. Then, to use <span style="font-variant:small-caps;">Kaltofen</span>’s algorithm, one must set a low enough probability of failure ${\varepsilon}$ yielding a big enough set of sampling points, see [@Kaltofen89factorizationof 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 $$\frac{6}{{\varepsilon}}\max{\left(2^{L+2},e(n-1)\,2^{e(n-1)}+e^3\,(n-1)^3,
  2\,(e(n-1)+1)^4\right)},$$ 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{\left(2^{L+2},e(n-1)\,2^{e(n-1)}+e^3\,(n-1)^3,2\,(e(n-1)+1)^4\right)}.$$ 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 \[th:ip\] for any characteristic, as in the introduction.

[th:ip]{} Let ${{\mathbf}{g}}=(g_1,\ldots,g_n) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$ be given in dense representation, and ${{\mathbf}{f}}={{\mathrm}{{\mathbf}{POW}}}_{n,d}=(x_1^d,\ldots,x_n^d) \in {{\mathbb}{K}}[x_1,\ldots,x_n]^n$ for some $d>0$. Whenever $\operatorname{char}{{\mathbb}{K}}=0$, let $e=d$ and $\tilde{{{\mathbf}{g}}}={{\mathbf}{g}}$. Otherwise, let $p=\operatorname{char}{{\mathbb}{K}}$, let $e$ and $r$ be integers such that $d=p^r\,e$, $p$ and $e$ coprime, and let $\tilde{{{\mathbf}{g}}}\in{{\mathbb}{K}}[{{\mathbf}{x}}]^n$ be such that ${{\mathbf}{g}}({{\mathbf}{x}})=\tilde{{{\mathbf}{g}}}({{\mathbf}{x}}^{p^r})$. Let $L$ be a polynomial size of an arithmetic circuit to evaluate the determinant of the Jacobian matrix of $\tilde{{{\mathbf}{g}}}$. If the size of ${{\mathbb}{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)\in\operatorname{GL}_n({{\mathbb}{K}})\times\operatorname{GL}_n({{\mathbb}{K}})$ such that: $${{\mathbf}{g}}=B\cdot{{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\cdot{{\mathbf}{x}}).$$

The computation of such a pair $(A,B)$ is summarized in the following Algorithm.

\[al:ip\][`IP`]{}for ${{\mathrm}{{\mathbf}{POW}}}_{n,d}$ and ${{\mathbf}{g}}$.

Input

:   One set of polynomials ${{\mathbf}{g}}=(g_1,\ldots,g_n)\in{{\mathbb}{K}}[x_1,\ldots,x_n]^n$, homogeneous of degree $d$.

Output

:   Two matrices $A,B\in\operatorname{GL}_n({{\mathbb}{K}})$ – if any – such that ${{\mathbf}{g}}({{\mathbf}{x}})=B\cdot{{\mathrm}{{\mathbf}{POW}}}_{n,d}(A\cdot{{\mathbf}{x}})$ with ${{\mathrm}{{\mathbf}{POW}}}_{n,d}({{\mathbf}{x}})=(x_1^d,\ldots,x_n^d)$ or “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

1.  **If** $\operatorname{char}{{\mathbb}{K}}=p>0$ **then**

    1.  Compute $r,e$ such that $d=p^r\,e$ with $p$ and $e$ coprime.

    2.  Compute $\tilde{{{\mathbf}{g}}}$ such that $\tilde{{{\mathbf}{g}}}({{\mathbf}{x}}^{p^r})={{\mathbf}{g}}({{\mathbf}{x}})$.

2.  **Else** $e=d$ and $\tilde{{{\mathbf}{g}}}={{\mathbf}{g}}$.

3.  Create an arithmetic circuit of polynomial size $L$ to evaluate $\det{{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}})$.

4.  Evaluate $\det{{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{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{{\mathrm}{J}}_{{{\mathbf}{g}}}({{\mathbf}{x}})$ with <span style="font-variant:small-caps;">Kaltofen</span>’s algorithm.

6.  **If** the factorization is $c\,\prod_{i=1}^n\ell_i({{\mathbf}{x}})^{e-1}$ **then** $A={\left(\ell_{i,j}^{e/d}\right)}_{1\leq i,j\leq n}$.

7.  **Else return** “[<span style="font-variant:small-caps;">NoSolution</span>]{}”.

8.  Compute $B$ such that $\tilde{{{\mathbf}{g}}}({{\mathbf}{x}})=B\cdot
        {\left(\ell_1({{\mathbf}{x}})^e,\ldots,\ell_n({{\mathbf}{x}})^e\right)}^{\operatorname{T}}$.

9.  **Return** $(A,B)$.

Square Root of a Matrix {#s:matsqrt}
=======================

In this section, we present further algorithms for computing the square root of a matrix. We use the same notation as in Section \[s:IP1S\]. 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 \[s:IP1S\] for computing the square root of a matrix in polynomial time for any field of characteristic $p\geq2$.

The square root as a polynomial in characteristic not $2$ {#ss:sqrt}
---------------------------------------------------------

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.

Let $Z\in{{\mathbb}{K}}^{n\times n}$ be an invertible matrix whose eigenvalues are $\zeta_1,\ldots,\zeta_r$. Let $\omega_1,\ldots,\omega_r$ be such that $\omega_i^2=\zeta_i$ and $\zeta_i=\zeta_j \Rightarrow
  \omega_i=\omega_j$, for all $1\leq i,j\leq r$. Then, there exists $W\in{{\mathbb}{K}}(\omega_1,\ldots,\omega_r)[Z]$ a square root of $Z$ whose eigenvalues are $\omega_1,\ldots,\omega_r$.

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_{\zeta,d}$ be the Jordan block of size $d$ associated with eigenvalue $\zeta$ and $\omega$ be a square root of $\zeta$. We shall first prove that the square root $G_{\omega,d}$ of $J_{\zeta,d}$ is a polynomial in $J_{\zeta,d}$ with coefficients in ${{\mathbb}{K}}(\omega)$. Matrix $J_{\zeta,d}-\zeta\,{{\mathrm}{Id}}_d$ is nilpotent of degree $d$. Hence, by the classical Taylor expansion of the square root near ${{\mathrm}{Id}}_d$, one can write $$\begin{aligned}
    G_{\omega,d}&=\omega\sum_{k=0}^{d-1}
    \binom{1/2}{k}\,\zeta^{-k}{\left(J_{\zeta,d}-\zeta\,{{\mathrm}{Id}}\right)}^k
    =\sum_{k=0}^{d-1}\binom{1/2}{k}\,\omega^{1-2\,k}
    {\left(J_{\zeta,d}-\zeta\,{{\mathrm}{Id}}\right)}^k
    =Q_{\zeta}(J_{\zeta,d})\label{eq:sqrtmat}\\
    &=\begin{pmatrix}
      \omega  &\binom{1/2}{1}\,\omega^{-1} &\cdots
      &\binom{1/2}{d-1}\,\omega^{3-2\,d}\\
      &\ddots &\ddots &\vdots\\
      &&\ddots &\binom{1/2}{1}\,\omega^{-1}\\
      &&&\omega
    \end{pmatrix},\notag
  \end{aligned}$$ with $Q_{\zeta}(x)=\sum_{k=0}^{d-1}\binom{1/2}{k}\,\omega^{1-2\,k}
  {\left(x-\zeta\right)}^i\in{{\mathbb}{K}}(\omega)[x]$.

It remains to prove that for multiple Jordan blocks, one can find a common polynomial. From equation , we deduce that $G$ is a polynomial in $J=\operatorname{Diag}{\left(J_{\zeta_1,d_1},\ldots,J_{\zeta_r,d_r}\right)}$ if and only if there exists a polynomial $Q$ such that $Q=Q_{\zeta_i} \mod
  (X-\zeta_i)^{d_i}$, for all $i$, $1\leq i\leq r$. By the Chinese Remainder Theorem, this can always be solved as soon as $\zeta_i=\zeta_j$ implies $Q_{\zeta_i}=Q_{\zeta_j}\bmod
  (X-\zeta_i)^{\min(d_i,d_j)}$, which is exactly the condition $\omega_i=\omega_j$.

Let us notice that picking the same square root for two equal eigenvalues is necessary. Indeed, although $W=\left(\begin{smallmatrix}1 &\phantom{-}0\\0
    &-1\end{smallmatrix}\right)$ is a square root of $Z=\left(\begin{smallmatrix}1 &0\\0 &1\end{smallmatrix}\right)$, $W\not\in{{\mathbb}{K}}[Z]$.

Matrices with square roots in characteristic $2$ {#ss:sqrt_char_2}
------------------------------------------------

In this part, we consider the trickier case of computing the square root of a matrix over a field ${{\mathbb}{K}}$ with $\operatorname{char}{{\mathbb}{K}}=2$. Unfortunately, unlike other characteristics, an invertible matrix has not necessarily a square root over $\bar{{{\mathbb}{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 $\binom{1/2}{k}$, involved in the Taylor expansion, are meaningless in characteristic $2$.

Let $Z\in{{\mathbb}{K}}^{n\times n}$ be a Jordan normal form with blocks $J_1,\ldots,J_r$ of sizes $d_1,\ldots,d_r\geq 2$, associated to eigenvalues $\zeta_1,\ldots,\zeta_r$ and blocks of sizes $1$ with eigenvalues $\upsilon_1,\ldots,\upsilon_s$. We assume that $J_1,\ldots,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\leq 1$ and $\zeta_1=\zeta_2$, $d_3-d_4\leq 1$ and $\zeta_3=\zeta_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 $\upsilon_j=\zeta_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 $\zeta$. Denoting $\omega$ the square root of $\zeta$, then $K$ is the square root of $J$ and $K_1',K_2'$ are those of $J'$ for $x,y,z$ any. $$\begin{aligned}
  J&=\begin{pmatrix}\zeta&0&0\\0&\zeta&1\\0&0&\zeta\end{pmatrix},\quad
  K=\begin{pmatrix}\omega&0&x\\\frac{1}{x}&\omega&y\\0&0&\omega
  \end{pmatrix},\\
  J'&=\begin{pmatrix}\zeta&1&0&0\\0&\zeta&0&0\\
    0&0&\zeta&1\\0&0&0&\zeta\end{pmatrix},\quad
  K_1'=\begin{pmatrix}\omega&x&0&y\\0&\omega&0&0\\
    \frac{1}{y}&z&\omega&x\\0&\frac{1}{y}&0&\omega\end{pmatrix},\quad
  K_2'=\begin{pmatrix}\omega&x&y&z\\0&\omega&0&y\\
    0&\frac{1}{y}&\omega&x\\0&0&0&\omega\end{pmatrix}.\end{aligned}$$ 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 $\zeta$ of sizes $1$ and $3$, and $J'''$, with three Jordan blocks associated with $\zeta$ of size $2$. Computing a square root of each of them yields an inconsistent system. $$\begin{aligned}
  J''&=\begin{pmatrix}\zeta&0&0&0\\0&\zeta&1&0\\
    0&0&\zeta&1\\0&0&0&\zeta\end{pmatrix},\quad
  J'''=\operatorname{Diag}{\left(\begin{pmatrix}\zeta&1\\0&\zeta\end{pmatrix},
    \begin{pmatrix}\zeta&1\\0&\zeta\end{pmatrix},
    \begin{pmatrix}\zeta&1\\0&\zeta\end{pmatrix}\right)}.\end{aligned}$$

Let $J$ be a Jordan block of size $d$ associated to eigenvalue $\zeta$. Then $J^2-\zeta^2{{\mathrm}{Id}}={\left(\begin{smallmatrix}{\mathbf}{0}&{{\mathrm}{Id}}_{d-2}\\
      {\mathbf}{0}&{\mathbf}{0}\end{smallmatrix}\right)}$ and one can deduce that $\zeta^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-\zeta{{\mathrm}{Id}})^d=0$ and $(J-\zeta{{\mathrm}{Id}})^e\neq 0$ for all $e<d$, then ${\left(J^2-\zeta^2{{\mathrm}{Id}}\right)}^{{\left\lceild/2\right\rceil}}=0$ and ${\left(J^2-\zeta^2{{\mathrm}{Id}}\right)}^e\neq 0$ for $e<{\left\lceild/2\right\rceil}$, *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 ${\left\lceild/2\right\rceil}$. 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$.

By this result, if $Z$ is a square, then one must be able to pair up its Jordan blocks with same eigenvalue $\zeta$ so that the sizes differ by at most $1$. The blocks that need not be paired being the blocks of size $1$.

Conversely, assuming one can pair up the Jordan blocks of $Z$ with same eigenvalue $\zeta$ 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 $\sqrt{\zeta}$. Furthermore, each lonely block of size $1$ associated with $\zeta$ is the square of the block of size $1$ associated with $\sqrt{\zeta}$.

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,\ldots,J_r$. For any polynomial $P$, $P(J)$ is also block diagonal with blocks $P(J_1),\ldots,P(J_r)$. Thus, if $P(J)^2=J$, then $P(J_i)^2=J_i$ for all $1\leq i\leq r$, which is false, unless $J_i$ has size $1$.

Computation in characteristic $p\geq2$ {#ss:sqrt_char_p}
--------------------------------------

In this part, we present an alternative method to the one presented in Section \[ss:computation\]. 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 $\operatorname{char}{{\mathbb}{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 ${{\mathbb}{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 $\bar{{{\mathbb}{K}}}$, then $J$ is block diagonal, otherwise it is a block upper triangular matrix. Its diagonal blocks are companion matrices ${{\mathcal}{C}}(P_1),\ldots, {{\mathcal}{C}}(P_r)$ of irreducible factors $P_1,\ldots,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 ${{\mathbb}{K}}$, see [@Matthews1992; @Storjohann1998]. Thus, after computing a square root $G$ of $J$, one can retrieve $W$ and $A$ of Section \[ss:computation\] in $O(n^{\omega})$ operations in the field of coefficients of $G$, with $\omega$ being the exponent of the time-complexity of matrix multiplication $2\leq\omega\leq 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 \[ss:sqrt\], $G$ is a polynomial in $J$. In the following, we only show how to determine the square root of a companion matrix ${{\mathcal}{C}}(P)$, for an irreducible $P$.

Let $P=x^d+p_{d-1}\,x^{d-1}+\cdots+p_{0}$, let us recall that the companion matrix of $P$ is $${{\mathcal}{C}}(P)=\begin{pmatrix}0 & & &-p_0\\1 &\ddots & &-p_1\\
  &\ddots &0&\vdots\\ & &1 &-p_{d-1}\end{pmatrix}.$$ If polynomial $P$ can be decomposed as $P(z)=(z-\alpha_0)\,\cdots\,(z-\alpha_{d-1})$, then we want to find a polynomial $Q$ such that $Q(z)=(z-\beta_0)\,\cdots\,(z-\beta_{d-1})$, where $\beta_i^2=\alpha_i$ for all $0\leq i\leq d-1$. Let us notice that $$P(z^2)=(z^2-\alpha_0)\,\cdots\,(z^2-\alpha_{d-1})
=Q(z)\,(z+\beta_0)\,\cdots\,(z+\beta_{d-1})
=(-1)^d\,Q(z)\,Q(-z).$$ As a consequence, the characteristic polynomial of ${{\mathcal}{C}}(Q)^2$ is $$\det(\lambda{{\mathrm}{Id}}-{{\mathcal}{C}}(Q)^2)=\det(\sqrt{\lambda}{{\mathrm}{Id}}-{{\mathcal}{C}}(Q))
\det(\sqrt{\lambda}{{\mathrm}{Id}}+{{\mathcal}{C}}(Q))=(-1)^dQ(\sqrt{\lambda})Q(-\sqrt{\lambda})
=P(\lambda).$$ But since $P$ is irreducible over ${{\mathbb}{K}}$, by the invariant factors theory, then ${{\mathcal}{C}}(Q)^2$ must be similar to the companion matrix ${{\mathcal}{C}}(P)$ over ${{\mathbb}{K}}$.

As $P$ is irreducible over ${{\mathbb}{K}}={{\mathbb}{F}}_q$, up to reindexing the roots of $P$, the conjugates $\alpha_1,\ldots,\alpha_{d-1}$ of $\alpha_0$ are just its iterated $q$th powers. Denoting ${{\mathbb}{L}}={{\mathbb}{K}}[x]/(P(x))={{\mathbb}{F}}_{q^d}$, let us assume that $S(y)=y^2-x$ is reducible in ${{\mathbb}{L}}[y]$, then $\beta_0\in{{\mathbb}{L}}$. As such, one can choose $\beta_i=\beta_0^{q^i}$, the iterated $q$th powers. In that case, the previous equations can be rewritten $$\begin{aligned}
  P(z)&=(z-\alpha_0)\,{\left(z-\alpha_0^q\right)}\,\cdots\,{\left(z-\alpha_0^{q^{d-1}}\right)}
  =(z-x)\,{\left(z-x^q\right)}\,\cdots\,{\left(z-x^{q^{d-1}}\right)},\\
  Q(z)&=(z-\beta_0)\,{\left(z-\beta_0^q\right)}\,\cdots\,{\left(z-\beta_0^{q^{d-1}}\right)}
  =(z-y)\,{\left(z-y^q\right)}\,\cdots\,{\left(z-y^{q^{d-1}}\right)}.\end{aligned}$$ As a consequence, $Q(z)\in{{\mathbb}{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 ${{\mathbb}{L}}$: $$u_0=x, u_1=x^q\bmod P(x),\ldots,u_{d-1}=u_{d-2}^q=x^{q^{d-1}}\bmod
P(x).$$ Then, we simply compute in $d$ operations $Q(z)=(z-u_0)\,(z-u_1)\,\cdots\,(z-u_{d-1})$ and we know that the resulting polynomial is in ${{\mathbb}{K}}[z]$.

Whenever $\alpha_0$ is not a square in ${{\mathbb}{L}}$, that is whenever $S(y)$ is irreducible, then $\beta_0^{q^d}$ is a square root of $\alpha_0$ different from $\beta_0$, thus it is $-\beta_0$. As a consequence, setting $Q(z)=(z-\beta_0)\,(z-\beta_0^q)\,\cdots\,(z-\beta_0^{q^{d-1}})$ would yield a polynomial that is not stable by the Frobenius endomorphism.

As such, we introduce a new variable $y$ to represent the field ${{\mathbb}{L}}'={{\mathbb}{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^{\frac{q^{i}-1}{2}}$, we can compute the following values in $O(d\log q)$ field operations in ${{\mathbb}{L}}$: $$u_0 = 1, u_1 = x^{\frac{q-1}{2}}\bmod P(x),\ldots,
u_{d-1}=u_1\,u_{d-2}^q=x^{\frac{q^{d-1}-1}{2}}\bmod P(x).$$ Consequently, $Q(z)=(z-yu_0)\,(z-yu_1)\,\cdots\,(z-yu_{d-1})$.

As a first step, we compute in $d$ operations, the dehomogenized polynomial in $y$, $$\tilde{Q}(z)=(z-u_{0}) (z-u_{1})\cdots (z- u_{d-1})
=z^d +h_{1} z^{d-1} + \cdots + h_{d-1} z + h_{d}.$$ Then, $Q(z)=z^d +y
h_{1} z^{d-1} + \cdots + y^{d-1}h_{d-1} z + y^{d}h_{d}$. But, denoting by $i_0=i\bmod 2$, we have $y^i=y^{i_0}
y^{i-i_0}=y^{i_0}x^{\frac{i-i_0}{2}}$. Hence we deduce: $$\begin{aligned}
  Q(z)&=z^d +y h_{1} z^{d-1} +xh_{2}z^{d-2}+yxh_{3}z^{d-3}+ \cdots
  +y^{d_0}x^{\frac{d-d_0}{2}}h_d\\
  &=z^d+y\,\sum_{i=0}^{{\left\lfloor\frac{d-1}{2}\right\rfloor}} h_{2\,i+1}\,
  x^{i}\,z^{d-2\,i-1} +\sum_{i=1}^{{\left\lfloor\frac{d}{2}\right\rfloor}}h_{2\,i}\,x^{i}\,
  z^{d-2\,i}.\end{aligned}$$

#### 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 ${{\mathbb}{L}}={{\mathbb}{F}}_{q^d}$, this is also bounded by $O(d\,{\mathsf}{M}(d)\,\log q)$ operations in ${{\mathbb}{K}}={{\mathbb}{F}}_q$, where ${\mathsf}{M}(n)$ is a bound on the number of operations in ${{\mathbb}{K}}$ to multiply two polynomials in ${{\mathbb}{K}}[x]$ of degree at most $n-1$. As a consequence, the computation of $W$ can be done in no more than $O(n^{\omega}+n{\mathsf}{M}(n)\log q)$ operations in ${{\mathbb}{K}}$. Let us assume that the characteristic polynomial of $Z$ has degree $n$ and can be factored as $P_1^{e_1}\cdots P_s^{e_s}$ with $P_i$ and $P_j$ coprime whenever $i\neq j$, $\deg P_i=d_i$ and $e_i\geq 1$. From a computation point of view, in the worst case, one needs to introduce a variable $\alpha_i$ for one root of $P_i$ and a variable $\beta_i$ for the square root of $\alpha_i$, assuming $\alpha_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 $P(z)=z^d+p_{d-1}z^{d-1}+\cdots+p_0=(z-\zeta_1)\cdots(z-\zeta_d)$, we want to compute $Q(z)=z^d+q_{d-1}z^{d-1}+\cdots+q_0=(z-\omega_1)\cdots(z-\omega_d)$, with $\omega_i^2=\zeta_i$ for all $1\leq i\leq d$. As $P(z^2)=Q(z)^2$, this yields $q_i=\sqrt{p_i}=p_i^{q/2}$, for all $1\leq i\leq
d-1$. Thus, $Q$ can be computed in $O(d\log q)$ operations in ${{\mathbb}{K}}$ and as a consequence, $W$ in $O(n^{\omega}+n\log q)$ operations in ${{\mathbb}{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 \[ss:sqrt\_char\_2\].

Acknowledgements {#acknowledgements .unnumbered}
----------------

We would like to thank Gabor <span style="font-variant:small-caps;">Ivanyos</span> for his helpful remarks and references on the irregular case. We wish to thank Gilles <span style="font-variant:small-caps;">Macario-Rat</span> for the many discussions about isomorphism of quadratic polynomials and Nitin <span style="font-variant:small-caps;">Saxena</span> 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 <span style="font-variant:small-caps;">ANR-11-BS02-0013 HPAC</span> project.

Proofs of \#[`IP1S`]{} {#ap:enum_proof}
======================

In this appendix, we shall prove the dimension of the centralizer of a matrix $J$, a Jordan normal form. This dimension, a consequence of [@Singla2010 Lemma 4.11], is used in Section \[s:enumeration\] 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\times s$. Indeed, if $X$ commutes with $J$, $XJ-JX$ is as such $$X\,J-J\,X={\left(\begin{smallmatrix}
    -x_{2,1} &x_{1,1}-x_{2,2} &\ldots &x_{1,n-1}-x_{2,n}\\
    \vdots &\vdots &&\vdots\\
    -x_{n,1} &x_{n-1,1}-x_{n,2}&\cdots&x_{n-1,n-1}-x_{n,n}\\
    0 &x_{n,1} &\ldots &x_{n,n-1}
  \end{smallmatrix}\right)}=0.$$

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

\[lm:cent\] Let $J$ be a Jordan normal form. For $1\leq i\leq r$, let us denote $J_i$ the $i$th block of $J$ and let us assume it is associated with eigenvalue $\zeta_i$ and it is of size $s_i$. Let $X=(X_{i,j})_{1\leq i,j\leq r}$ be a block-matrix, with $X_{i,j}\in {{{\mathbb}{L}}{\left(\zeta_1,\ldots,\zeta_r\right)}}^{s_i\times s_j}$, that commutes with $J$. If $\zeta_i=\zeta_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$.

We assume that $r=2$. If $X\,J-J\,X={\left(\begin{smallmatrix}X_{1,1}\,J_1-J_1\,X_{1,1}
      &X_{1,2}\,J_2-J_1\,X_{1,2}\\
      X_{2,1}\,J_1-J_2\,X_{2,1}
      &X_{2,2}\,J_1-J_1\,X_{2,2}\end{smallmatrix}\right)}=0$, then $X_{1,1}$ commutes with $J_1=J_{\zeta_1,s_1}$ and $X_{2,2}$ with $J_2=J_{\zeta_2,s_2}$. Thus they are upper triangular Toeplitz matrices.

From $X_{2,1}\,J_2-J_1\,X_{2,2}$, one deduces that $(\zeta_1-\zeta_2)\,x_{s_1+s_2,1}=0$, hence either $\zeta_1=\zeta_2$ or $x_{s_1+s_2,1}=0$. If $\zeta_1\neq\zeta_2$, then step by step, one has $X_{1,2}=0$. Assuming $\zeta_1=\zeta_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 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.

\[lm:dim\_cent\] Let $H\in{{\mathbb}{K}}^{n\times n}$ be a matrix and let $J$ be its normal Jordan form. Assuming the blocks of $J$ associated to $\zeta_i$ are $J_{\zeta_i,s_{i,1}},\ldots,J_{\zeta_1,s_{i,d_1}}$ with $s_{i,1}\leq\cdots\leq s_{i,d_i}$ for $i,1\leq i\leq r$, then the centralizer of $H$ is a ${{\mathbb}{K}}$-vector subspace of ${{\mathbb}{K}}^{n\times n}$ of dimension no more than $\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i}(2\,d_i-2\,j+1)s_{i,j}$.

Let ${{\mathbb}{L}}$ be the smallest field over which $J$ is defined. It is clear that the centralizer of $H$ over ${{\mathbb}{L}}$, denoted ${\mathcal}{W}$, contains $\operatorname{{\mathcal}{C}}(H)\otimes{{\mathbb}{L}}$. Hence, $\dim_{{{\mathbb}{K}}}\operatorname{{\mathcal}{C}}(H)=\dim_{{{\mathbb}{L}}}(\operatorname{{\mathcal}{C}}(H)\otimes{{\mathbb}{L}})
  \leq\dim_{{{\mathbb}{L}}}{\mathcal}{W}$.

Now, let $X=(X_{i,j})_{1\leq i,j\leq
    d_1+\cdots+d_r}\in{\mathcal}{V}$. From Lemma \[lm:cent\], there are $\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i} s_{i,j}$ free parameters for the diagonal blocks of $X$ and $2\,\sum_{1\leq i\leq
    r}\sum_{1\leq j<k\leq d_i}\min(s_{i,j}, s_{i,k})
  =2\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i}(d_i-j)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 $\operatorname{{\mathcal}{C}}(H)$ is bounded from above by $$q^{{\left(\sum_{1\leq i\leq r}\sum_{1\leq j\leq d_i}(2\,d_i-2\,j+1)\,s_{i,j}\right)}}-1,$$ as stated in Corollary \[cor:enum\].

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