# Faster Modular Composition

Vincent Neiger, Bruno Salvy, Éric Schost, Gilles Villard

► To cite this version:

Vincent Neiger, Bruno Salvy, Éric Schost, Gilles Villard. Faster Modular Composition. Journal of the ACM (JACM), 2023. hal-03380258

HAL Id: hal-03380258

https://hal.science/hal-03380258

Submitted on 15 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Faster Modular Composition

VINCENT NEIGER, Sorbonne Université, France

BRUNO SALVY, Inria, France

ÉRIC SCHOST, University of Waterloo, Canada

GILLES VILLARD, CNRS, France

A new Las Vegas algorithm is presented for the composition of two polynomials modulo a third one, over an arbitrary field. When the degrees of these polynomials are bounded by $n$, the algorithm uses $O(n^{1.43})$ field operations, breaking through the 3/2 barrier in the exponent for the first time. The previous fastest algebraic algorithms, due to Brent and Kung in 1978, require $O(n^{1.63})$ field operations in general, and $n^{3/2+o(1)}$ field operations in the particular case of power series over a field of large enough characteristic. If using cubic-time matrix multiplication, the new algorithm runs in $n^{5/3+o(1)}$ operations, while previous ones run in $O(n^2)$ operations.

Our approach relies on the computation of a matrix of algebraic relations that is typically of small size. Randomization is used to reduce arbitrary input to this favorable situation.

CCS Concepts: • *Computing methodologies* → *Algebraic algorithms*; • *Theory of computation* → *Algebraic complexity theory*.

Additional Key Words and Phrases: composition of polynomials, complexity

## 1 INTRODUCTION

Many fundamental operations over univariate polynomials of degree at most $n$ with coefficients in a commutative ring $\mathbb{A}$ can be computed in a number of arithmetic operations in $\mathbb{A}$ that is quasi-linear in $n$ [22]. It is the case for multiplication, division with remainder by a monic polynomial, multipoint evaluation, interpolation at points whose differences are units in $\mathbb{A}$, and greatest common divisors when $\mathbb{A}$ is a field.

In contrast with these operations, improving the cost bound for *modular composition* is a long-standing open question. Given three polynomials $a, f \in \mathbb{A}[x]$ and $g \in \mathbb{A}[y]$, with $\deg(a) < n$ and $\deg(g) < n$ where $n = \deg(f)$, and with $f$ monic, this problem is to compute $g(a)$ $\mathrm{rem} f$, where the "rem" operation takes the remainder of the Euclidean division.

This operation arises in a variety of contexts. For instance, with $f = x^n$, it amounts to power series composition. Power series reversion can then be reduced to composition, with a small overhead [15], as can further operations such as solving families of functional equations [28].

The application of certain algebra morphisms also translates to modular composition. Over a field $\mathbb{K}$, for $f$ and $a$ in $\mathbb{K}[x]$, we denote by $a \mod f \in \mathbb{K}[x]/\langle f \rangle$ the class of $a$ modulo $f$. Then, for $e$ and $f$ in respectively $\mathbb{K}[y]$ and $\mathbb{K}[x]$, and for an $\mathbb{K}$-algebra morphism $\phi: \mathbb{K}[y]/\langle e \rangle \to \mathbb{K}[x]/\langle f \rangle$, if $\phi(y \mod e) = a \mod f$ then for $g$ in $\mathbb{K}[y]$, the image $\phi(g \mod e)$ is equal to $g(a) \mod f$.

Over finite fields, with $e$ and $f$ the same polynomial and $\phi$ the Frobenius endomorphism, this results in modular composition playing an important role in algorithms for polynomial factorization [23, 44, 45]. Dedicated algorithms exist for modular composition over finite fields, with quasi-linear complexity (they are discussed later), but there remains a variety of questions that can be considered over arbitrary fields, and which are impacted by modular composition (or closely related operations such as power projection, discussed later as well): computing the minimal

<footnote>Authors' addresses: Vincent Neiger, Sorbonne Université, Laboratoire LIP6 UMR 7606 CNRS, Sorbonne Université, Paris, France, vincent.neiger@lip6.fr; Bruno Salvy, Inria, Laboratoire LIP UMR 5668 Univ. Lyon, CNRS, ENS de Lyon, Inria, UCBL, Lyon, France, bruno.salvy@inria.fr; Éric Schost, University of Waterloo, Cheriton School of Computer Science, Waterloo, ON, N2L 3G1, Canada, eric.schost@uwaterloo.ca; Gilles Villard, CNRS, Laboratoire LIP UMR 5668 Univ. Lyon, CNRS, ENS de Lyon, Inria, UCBL, Lyon, France, gilles.villard@ens-lyon.fr.</footnote>
polynomial of an algebraic number [65–67], normal bases computations [26, 45], arithmetic operations with two algebraic numbers [10], computing with towers of algebraic extensions [32, 61, 62], Riemann-Roch space computations [1, 2], etc.

Modular composition can be performed using Horner's algorithm with modular reduction at each stage, which leads to a complexity in $\tilde{O}(n^2)$ operations if fast polynomial multiplication is used. (The notation $c' = \tilde{O}(c)$ means that $c' = O(c \log^k(c))$ for some $k > 0$. In other words, logarithmic factors are dropped.) In 1978, Brent and Kung gave two algorithms that perform composition modulo $x^n$ (the case of power series) [14, 15]. One relies strongly on Taylor expansion and runs in $\tilde{O}(n^{3/2})$ operations; the other one, using a baby steps/giant steps approach, uses $O(n^{(\omega+1)/2}) + \tilde{O}(n^{3/2})$ operations, where $\omega \le 3$ is a feasible matrix multiplication exponent (two $n \times n$ matrices can be multiplied in $O(n^\omega)$ operations). This second algorithm works verbatim and in the same complexity for composition modulo an arbitrary polynomial $f$ of degree $n$ not restricted to be $x^n$ [23].

Both have remained essentially the best algorithms since then. Huang and Pan used fast rectangular matrix multiplication in the central step of the baby steps/giant steps algorithm to reduce its complexity to $O(n^{\omega_2/2}) + \tilde{O}(n^{3/2})$ where $\omega_2 \le \omega + 1$ is a feasible exponent such that a $n \times n^2$ matrix can be multiplied by a $n^2 \times n$ matrix in $O(n^{\omega_2})$ operations. The currently best known value gives $\omega_2 \approx 3.25$ [52], which makes the previous algebraic complexity bound $O(n^{1.626})$ for modular composition for an arbitrary $f$. Even assuming an optimal matrix multiplication, which means $\omega = 2$, these algorithms do not break the exponent barrier $3/2$.

The open problem 2.4 in the book of Bürgisser, Clausen and Shokrollahi [17] asks whether Brent and Kung's algorithm can be improved substantially. The research problem 12.19 in von zur Gathen and Gerhard's book [22] asks for a complexity in $\tilde{O}(n^{1.5})$ or better. Our main result answers both questions positively when $\mathbb{A}$ is a field, with few extra hypotheses.

**THEOREM 1.1.** Given $a, f \in \mathbb{K}[x]$ and $g \in \mathbb{K}[y]$ with coefficients in a field $\mathbb{K}$, with $\deg(f) = n$, $\deg(a)$ and $\deg(g)$ smaller than $n$, and a tuple $r$ of $O(n^{1+1/3})$ field elements, Algorithm MODULARCOMPOSITION returns either $g(a)$ rem $f$ or FAIL after $\tilde{O}(n^\kappa)$ arithmetic operations in $\mathbb{K}$, with

$$ \kappa = 1 + \frac{1}{\frac{1}{\omega-1} + \frac{2}{\omega_2-2}} < 1.43. \quad (1) $$

It returns FAIL with probability at most $(2n^4 + 18n^2)/\text{card}(S)$ when the entries of $r$ are chosen uniformly and independently from a finite subset $S \subseteq \mathbb{K}$.

Here we use an algebraic model of computation: roughly, basic arithmetic operations $\{+, -, \times, \div\}$ and zero-tests in the base field $\mathbb{K}$ are counted at unit cost; for more details, see Section 2. As usual with probabilistic algorithms of Las Vegas type, the algorithm can be repeated until it succeeds, so that only its running time becomes a random variable.

We assume that the characteristic $p$ of $\mathbb{K}$ is known to the algorithms. For $\mathbb{K}$ finite of small cardinality $q$ (namely, $q \le 2n^4 + 18n^2$), the probability statement becomes vacuous. However, in such cases, one can work in a sufficiently large field, by constructing an extension of $\mathbb{K}$ of degree $O(\log(n))$ efficiently (see [22, Sec. 14.9] and references therein). In this extension, each arithmetic operation can be performed in $\tilde{O}(\log(n))$ arithmetic operations in $\mathbb{K}$, so that the asymptotic complexity estimate is unaffected.

We also give a probabilistic algorithm of the Las Vegas type with the same complexity bound for computing an annihilating polynomial for $a \mod f$, that is, a nonzero polynomial $g \in \mathbb{K}[y]$ such that $g(a)$ rem $f = 0$.

To compute $g(a)$ rem $f$, our approach relies on the following main ingredients. Properties of “generic” inputs are mentioned; throughout the article, genericity is understood in the Zariski sense: a property is generic if it holds outside of a hypersurface of the corresponding parameter space.
(1) The existence, under genericity conditions, of sufficiently many independent elements of “small degree” in the ideal $\langle y - a(x), f(x) \rangle$ generated by the polynomials $y - a(x)$ and $f(x)$ in $K[x, y]$. More precisely, we consider *matrices of relations*: these are nonsingular matrices in $K[y]^{m \times m}$ (the choice of $m$ is optimized below) with entries of degree smaller than a certain integer $d$, and whose columns are vectors of coefficients with respect to $x$ of polynomials $r \in K[x, y]$ such that $r(x, a)$ rem $f = 0$. The “small degree” condition is that $d \le \lceil n/m \rceil$. (See Sections 4 and 7.)

(2) Fast algorithms on polynomial matrices. We use *approximant bases* [27] to compute a matrix of relations from sufficiently many coefficients of the polynomials $x^i a^k$ rem $f$. Once we have such a matrix, linear system solving over $K[x]$ [75] allows us to reduce the univariate $g \in K[y]$ to a bivariate $\tilde{g} \in K[x, y]$, of degrees smaller than $m$ and $d$ in $x$ and $y$, such that $g(a) \equiv \tilde{g}(x, a) \mod f$. These operations have complexity $\tilde{O}(m^\omega d)$. (See Sections 4.2 and 5.3.)

(3) An algorithm that computes the first $m$ coefficients of the $2md$ polynomials $x^i a^k$ rem $f$, for $0 \le i < m$ and $0 \le k < 2d$; those are needed above, to compute matrices of relations (see Section 3.3). It uses $\tilde{O}(m^2 d + c(n, m, d))$ operations in $K$, with

$$c(n, m, d) = (m + n/d)d^{\omega_2/2}. \qquad (2)$$

(4) A generalization due to Nüsken and Ziegler [59] of Brent and Kung's algorithm to the case of a bivariate polynomial $\tilde{g}(x, y) \in K[x, y]$. If $m$ and $d$ are the degrees of $\tilde{g}$ in $x$ and $y$, their algorithm computes the “uni-bivariate” composition $\tilde{g}(x, a)$ rem $f$ using $\tilde{O}(c(n, m, d))$ operations in $K$. (See Section 3.2.)

(5) A *randomized change of basis* to bring $f$ and $a$ to a situation where “small” matrices of relations exist. (See Section 8.)

(6) Correctness of the resulting randomized algorithm is established only for $f$ separable (i.e., with no repeated roots in an algebraic closure $\overline{K}$ of $K$) and $f$ purely inseparable (i.e., with only one root in $\overline{K}$). Modular composition modulo an arbitrary $f$ is reduced to these two extreme cases by *separable decomposition* of $f$ [53], Chinese remainder theorem, and a slight generalization of a technique called *untangling* [30]. The latter allows to transport the composition problem over $K$ modulo a factor of the separable decomposition, to a composition problem over a quotient algebra with purely inseparable modulus. (See Section 9.)

Altogether, the costs of the various parts of the algorithm add up to $\tilde{O}(m^\omega d + c(n, m, d))$ operations in $K$. As mentioned above, for a given $m$, we rely on matrices of relations of degree $d = \lceil n/m \rceil$. Then choosing $m$ and $d$ so as to minimize $m^\omega d + c(n, m, d)$ leads us to $m \sim n^\eta$, where

$$\eta = \frac{1}{1 + \frac{\omega - 1}{(\omega_2 - 2)/2}}, \qquad (3)$$

and to the complexity estimate $\tilde{O}(n^\kappa)$ of Theorem 1.1. Using the known bounds $\omega < 2.37286$ [3, 51] and $\omega_2 < 3.251640$ [52] gives $\eta \approx 0.3131$.

The improvements brought by fast matrix multiplication on one hand and by fast rectangular matrix multiplication on the other hand are made clearer by noting that the exponent $\kappa$ of composition satisfies

$$\frac{4}{3} \le \kappa = 1 + \underbrace{\frac{1}{\frac{1}{\omega-1} + \frac{2}{\omega_2-2}}}_{< 1.429866} \le \underbrace{\frac{\omega+2}{3}}_{< 1.45762} \le \underbrace{\frac{5}{3}}_{< 1.666667},$$

where the first approximation is obtained with the bounds on $\omega$ and $\omega_2$ given above; the second one is obtained when no fast rectangular matrix multiplication is used, so that $\omega_2$ simply becomes $\omega + 1$; the last one is obtained when no fast matrix multiplication is used, thus taking $\omega = 3$. In
the latter case, our algorithm is the first subquadratic one for modular composition. In the other direction, considering the lower bounds $\omega \ge 2$ and $\omega_2 \ge 3$ shows that $\kappa \ge 4/3$, giving a lower bound on the complexity estimate that can be achieved by the algorithm designed in this work.

## 1.1 Algorithmic Tools

Our work builds upon a sequence of earlier algorithmic progress that we now recall. We denote by $\mathbb{K}[x]_{<n}$ the set of univariate polynomials in $x$ with coefficients in $\mathbb{K}$ and degree less than $n$; by $\mathbb{K}[x, y]_{< (r, s)}$ the bivariate polynomials in $x, y$ of bi-degree in $(x, y)$ less than $(r, s)$.

**1.1.1 Baby steps/giant steps.** One of the bottlenecks in algebraic approaches for evaluating $g$ at a modulo $f$ is the computation of successive powers $1, a, a^2, \dots$ modulo $f$, which leads to the question of minimizing the number of powers that are used. The solution used by Brent and Kung relies on a baby steps/giant steps scheme [15, 60], where only

$$1, a, \dots, a^{\lceil \sqrt{n} \rceil} \text{ rem } f \quad \text{and} \quad a^{2^{\lceil \sqrt{n} \rceil}}, a^{3^{\lceil \sqrt{n} \rceil}}, \dots \text{ rem } f$$

are computed. The former group forms the baby steps; the latter forms the giant steps. The problem is then reduced to about $\sqrt{n}$ modular compositions “$g_i(a)$ rem $f$” for $g_i$ of degree about $\sqrt{n}$. These compositions are all obtained simultaneously through the multiplication of two matrices of sizes roughly $\sqrt{n} \times \sqrt{n}$ and $\sqrt{n} \times n$. This is followed by a less expensive Horner evaluation step using the powers of $a^{\lceil \sqrt{n} \rceil}$. See Section 3.1 for a complete description.

**1.1.2 Projection-Reconstruction.** Wiedemann's algorithm [72] finds the minimal polynomial of a matrix $A \in \mathbb{K}^{n \times n}$ by considering the sequence $(v^T A^k w)_{k \ge 0}$, for two vectors $v$ and $w$. This sequence is linearly recurrent and its generating function $h(y) = \sum_{k \ge 0} (v^T A^k w)/y^{k+1}$ is rational; for generic $v$ and $w$, the denominator of $h(y)$ is the minimal polynomial $\mu_A$ of the matrix $A$. Writing $d \le n$ for the degree of $\mu_A$, this polynomial can be reconstructed efficiently from the first $2d$ terms of the sequence by the Berlekamp-Massey algorithm or, equivalently, by the computation of a Padé approximant. Given the expansion in $y^{-1}$ of a rational power series $h(y) = q(y)/\mu_A(y)$ with polynomials $q$ and $\mu_A$ of degree at most $d-1$ and $d$, this reconstructs the fraction $(q, \mu_A)$ as a solution of

$$(h(y) + O(y^{-2d-1}))\mu_A(y) - q(y) = O(y^{-d-1}).$$

If the degree of $\mu_A$ is unknown, one can use this approach with the upper bound $d=n$ instead.

Wiedemann's algorithm can be combined with the baby steps/giant steps paradigm [40, Sec. 3; 65; 45, Algorithm AP]. In particular, when $A$ is the matrix $M_a$ of multiplication by $a$ mod $f$ in the basis $\mathcal{B} = (1, x, \dots, x^{n-1})$ of $\mathbb{K}[x]/\langle f \rangle$, this was used by Shoup to compute the minimal polynomial of the polynomial $a$ modulo $f$ [65–67]. For irreducible $f$, Shoup used the vectors $v = w = 1$ (where 1 is the first column of the identity matrix), in which case the sequence $(v^T A^k w)_k$ becomes the sequence of power projections $(\ell(1), \ell(a), \ell(a^2), \dots)$, where $\ell$ is the linear form which takes the coefficient of 1 of an element of $\mathbb{K}[x]/\langle f \rangle$ written on the basis $\mathcal{B}$. For an arbitrary $f$, Shoup used a random linear form $\ell$, corresponding to a random choice of the vector $v$ and $w=1$.

In either case, the required $2d$ elements of the sequence can be obtained by left multiplication by $v^T$ of a matrix whose columns are the coefficient vectors of $1, a, a^2, \dots$ modulo $f$. Now, the right multiplication of the exact same matrix by a vector of coefficients corresponds to modular composition. Using the transposition principle, Shoup described a baby steps/giant steps algorithm which computes the power projections for an arbitrary linear form $\ell: \mathbb{K}[x]/\langle f \rangle \to \mathbb{K}$ in the same complexity as that of Brent and Kung's algorithm [65–67]. (See Section 3.1.2.) This principle states that the existence of an algebraic algorithm for the multiplication of a matrix by a vector induces the existence of an algorithm for the product of the transpose of that matrix by a vector, both having essentially the same complexity [17, Thm. 13.20; 11].
The same idea is used by Shoup for another operation that we also need. Given a, b, f, the inverse modular composition asks for a polynomial g of least degree such that $g(a) \equiv b \mod f$ or for a proof that no such g exists. This problem reduces to the computation of the power projections

$$ (\ell(1), \ell(a), \dots, \ell(a^{2n-1})) \text{ and } (\ell(b), \ell(ab), \dots, \ell(a^{n-1}b)), $$

again in the same complexity as that of modular composition, followed by the resolution of a linear system of Hankel type [65, Thm. 3.5]. The latter is known to be equivalent to Padé approximation [13], where the equation becomes

$$ \left( \sum_{k \ge 0} \frac{\ell(a^k)}{y^{k+1}} + O(y^{-2n-1}) \right) g(y) - q(y) = \sum_{k \ge 0} \frac{\ell(a^k b)}{y^{k+1}} + O(y^{-n-1}), $$

to be solved for a numerator $q(y) \in \mathbb{K}[y]_{<n}$ and the inverse composition $g(y) \in \mathbb{K}[y]_{\le n}$.

### 1.1.3 Blocks for speed and structure.
Coppersmith introduced a block version of Wiedemann's algorithm [18]. There, the scalar sequence $(v^T A^k w)_{k \ge 0}$ is replaced by the matrix sequence $(V^T A^k W)_{k \ge 0}$ for two matrices $V \in \mathbb{K}^{n \times \ell}$ and $W \in \mathbb{K}^{n \times m}$; the generating function $H(y) = \sum_{k \ge 0} (V^T A^k W)/y^{k+1}$ is a rational $\ell \times m$ matrix.

Such a matrix admits an irreducible matrix fraction description $N(y)D(y)^{-1}$ with $N \in \mathbb{K}[y]^{\ell \times m}$ and $D \in \mathbb{K}[y]^{m \times m}$ two polynomial matrices (see Section 5.1.1), and the columns of the denominator matrix $D$ form a basis of the $\mathbb{K}[y]$-module of polynomial vectors $u \in \mathbb{K}[y]^m$ such that $\sum_{i \le \deg(u)} V^T A^{k+i} W u_i = 0$ for all $k \ge 0$, where $u_i$ denotes the coefficient of $y^i$ in $u$ [46, Lem. 2.8]. For $m=1$, this module is the ideal generated by the minimal polynomial of the sequence in Wiedemann's algorithm. For $1 \le m \le \ell \le n$, the matrix $D$ contains more information: for example, for generic $V$ and $W$, its invariant factors are the $m$ invariant factors of largest degree of the characteristic matrix $yI_n - A$ [46, Thm. 2.12], the highest degree one being the minimal polynomial of $A$. Consequently, the determinant of $D$ has degree the sum $v_m$ of the degrees of these $m$ invariant factors, which implies that $v_m \le n$.

The computation of $D$ can be achieved in two steps, which are matrix versions of the methods used for $m=1$ in Section 1.1.2. Writing $d$ for the degree of $D$, it is sufficient to compute the first $2d$ matrices of the sequence $(V^T A^k W)_{k \ge 0}$, which can be done by a baby steps/giant steps approach [46]. Next, $D$ is obtained by matrix fraction reconstruction, solving

$$ (V^T(yI_n - A)^{-1}W + O(y^{-2d-1}))D(y) - N(y) = O(y^{-d-1}) $$

for the unknown $N \in \mathbb{K}[y]^{\ell \times m}$ and $D \in \mathbb{K}[y]^{m \times m}$ of degrees at most $d-1$ and $d$; this can be done efficiently by a generalization of Padé approximation called minimal approximant bases, whose properties are recalled in Section 5.2. (See [46, 47] for bibliographic pointers to algorithms that compute minimal linear generators of matrix sequences.) The parameter $d$ plays a major role in the efficiency of both steps: it is usually unknown a priori, and might be as large as $\Theta(n)$. Yet, the interest of this block approach lies in the fact that, for generic $V$ and $W$ and $\ell \ge m$, the matrix $D$ has degree $d = \lceil v_m/m \rceil \le \lfloor n/m \rfloor$ [70, Cor. 6.4].

### 1.1.4 Efficient projections and small bivariate polynomials.
Special choices of the matrices $V$ and $W$ above, with identity blocks, lead to efficient projections and have been shown to be effective in the context of black-box matrix inversion [20]. Even simpler matrices, $X = (I_m \ 0)^T$ and $Y = (0 \ I_m)^T$ in $\mathbb{K}^{n \times m}$ with $m \in \{1, \dots, n\}$, have been used by Villard in his fast algorithm for the bivariate resultant of two bivariate polynomials $f$ and $g$ in $\mathbb{K}[x, y]$ [71]. In this context, for generic $f$ and $g$, this choice of $X$ and $Y$ is sufficient to ensure that the denominator matrix $D$ contains $m$ "small" polynomials in the ideal of $\mathbb{K}[x, y]$ generated by $f$ and $g$.
## 1.2 Overview of the core algorithm

When a mod $f$ has a minimal polynomial $\mu_a$ of small degree, $\mu_a$ can be computed efficiently using power projections $(\ell(1), \ell(a), \ell(a^2), \dots)$ by Shoup's algorithm, since few terms in the sequence are needed (see Sections 1.1.2 and 3.1.3). Then, for composition, one uses the identity $g(a) \equiv \hat{g}(a) \mod f$, where $\hat{g} = g \operatorname{rem} \mu_a$. Since $\hat{g}$ has small degree, this reduces the number of powers of a mod $f$ that need be considered.

Our algorithm can be viewed as a block or bivariate version of this approach, replacing the univariate polynomial $\mu_a$ by a collection of $m$ small bivariate polynomials in the ideal generated by $y - a(x)$ and $f(x)$, for a fixed parameter $m$. In a generic situation, while $\mu_a$ has degree $n$, there exists such a collection with degrees $m-1$ and $[n/m]$ in $x$ and $y$. This collection is represented as a matrix in $\mathbb{K}[y]^{m \times m}$ and is found efficiently by exploiting the structure of the matrix of multiplication by a mod $f$.

*Matrices of relations.* Let $M_a \in \mathbb{K}^{n \times n}$ be the matrix of multiplication by a mod $f$ in the basis $(1, x, \dots, x^{n-1})$. Following Section 1.1.3, in the special case where $A = M_a$, if $V$ is a generic matrix in $\mathbb{K}^{n \times l}$, and $W$ is the matrix $X = (I_m \ 0)^T$ with $m \le l$ and $m \le n$, the block Wiedemann approach yields a denominator matrix $D \in \mathbb{K}[y]^{m \times m}$ whose columns represent a basis of the $\mathbb{K}[y]$-module

$$ \mathcal{M}_m^{(a,f)} = \{r(x, y) = r_0(y) + \dots + r_{m-1}(y)x^{m-1} \mid r(x, a(x)) \equiv 0 \mod f(x)\}; \quad (4) $$

this follows for instance from [70, Lem. 4.2]. The elements of this module are algebraic relations of degree less than $m$ in $x$ satisfied by a mod $f$ (Section 5.1). We call *matrix of relations* any nonsingular matrix $R_m^{(a,f)} \in \mathbb{K}[y]^{m \times m}$ whose columns are the coefficients of polynomials in $\mathcal{M}_m^{(a,f)}$ (Section 4.1.1), that is, any nonsingular right multiple of $D$. Given a matrix of relations $R_m^{(a,f)}$, the composition $g(a) \operatorname{rem} f$ is obtained in two steps.

* First, by *polynomial matrix division* [39, Thm. 6.3-15, p. 389], there exist vectors $v, w \in \mathbb{K}[y]^m$ such that

$$ (g(y) \ 0 \ \dots \ 0)^T = R_m^{(a,f)} w + v, \quad (5) $$

where $\deg(v) < d$ and $d$ is an upper bound on $\deg(R_m^{(a,f)})$; finding such vectors takes $\tilde{O}(m^\omega(d+n/m))$ operations (Section 4.2) [75]. Then, by design, the bivariate polynomial

$$ \tilde{g}(x, y) = v_1(y) + \dots + v_m(y)x^{m-1} $$

has degree less than $m$ and $d$ in $x$ and $y$, and is such that $g(a) \equiv \tilde{g}(x, a) \mod f$.

* The polynomial $\tilde{g}$ can then be evaluated at $y = a \mod f$ by the Nüsken-Ziegler algorithm in $\tilde{O}(c(n, m, d))$ operations, with $c(\cdot)$ from Eq. (2) (Proposition 3.4).

*Truncated sequence of projections.* In the block Wiedemann approach, using $X$ as our right projection matrix, we need the first $2d$ elements of the matrix sequence $(VM_a^k X)_k$, which amounts to a type of bivariate power projections (see Section 1.6.2). Unfortunately, we do not know how to obtain them efficiently enough for an arbitrary $V$. Choosing $V = X^T$, we design a baby steps/giant steps algorithm in Section 3.3 that runs in $\tilde{O}(c(n, m, d) + m^2d)$ operations. With this choice, by fraction reconstruction the sequence $(X^T M_a^k X)_k \ge 0$ yields a denominator $D$ which is a basis of the $\mathbb{K}[y]$-module

$$ \mathcal{M}_{m,m}^{(a,f)} = \{r(x, y) \in \mathbb{K}[x, y]_{< (m, \cdot)} \mid [a(x)^k r(x, a(x)) \operatorname{rem} f(x)]_0^{m-1} = 0 \text{ for all } k \ge 0\}, $$

where $[\cdot]_0^{m-1}$ is the projection on $\operatorname{Span}(1, x, \dots, x^{m-1})$. The inclusion $\mathcal{M}_m^{(a,f)} \subseteq \mathcal{M}_{m,m}^{(a,f)}$ holds but may be strict, leading to a denominator $D$ that is not a matrix of relations.
*Matrices of relations of small degree.* For an arbitrary $f$ with $f(0) \neq 0$ (this is not really a restriction, see Remark 3.8) and a generic $a$, two important properties hold (see Section 7.3): the above inclusion of modules is an equality — making the algorithm correct — and a basis $R_m^{(a,f)}$ of degree $d = \lceil n/m \rceil$ of $M_m^{(a,f)}$ can be reconstructed from the first $2d$ elements of the sequence $(X^T M_a^k X)_k \in \mathbb{N}$ — making the algorithm fast.

The reconstruction is done via minimal approximant bases in Sections 5.2 and 5.3. Directly extending Section 1.1.3, we would solve the equation at infinity

$$ (X^T(yI_n - M_a)^{-1}X + O(y^{-2d-1})) R_m^{(a,f)}(y) - N(y) = O(y^{-d-1}), \quad (6) $$

for unknown matrices $N$ and $R_m^{(a,f)}$ of degree at most $d-1$ and $d$. For technical reasons coming from the reconstruction algorithm, we actually use an expansion at $y=0$ rather than at infinity, so that the sequence we use involves powers of $M_a^{-1}$ instead of $M_a$ (see Remark 5.7).

Beyond generic cases, a relevant quantity is

$$ v_m^{(a,f)} = \deg(\sigma_1) + \dots + \deg(\sigma_m), \quad (7) $$

where $\sigma_1, \dots, \sigma_n \in \mathbb{K}[y]$ are the invariant factors of $yI_n - M_a$, ordered by decreasing degree. This quantity is at most $n$, and it is the degree of the determinant of any basis of $M_m^{(a,f)}$ (Proposition 4.1). In favorable situations, working with $d = \lceil v_m^{(a,f)}/m \rceil$, and $a$ fortiori with $\lceil n/m \rceil$, is sufficient to obtain such a basis $M_m^{(a,f)}$.

## 1.3 Probabilistic algorithm for $f$ separable or purely inseparable

Our probabilistic algorithm aims at bringing arbitrary inputs to the favorable situation mentioned above, by means of a random change of basis. For a polynomial $\gamma \in \mathbb{K}[x]$ such that the minimal polynomial $\mu_\gamma$ of $\gamma$ mod $f$ has degree $n$, the powers $(1, \gamma, \dots, \gamma^{n-1})$ mod $f$ form a basis of $\mathbb{A} = \mathbb{K}[x]/\langle f \rangle$. This induces a $\mathbb{K}$-algebra isomorphism:

$$ \phi_\gamma : \mathbb{A} \to \mathbb{K}[y]/\langle \mu_\gamma \rangle $$

that maps $\gamma$ to $y$, and more generally $u \in \mathbb{A}$ to $v$ such that $v(\gamma) \equiv u \mod f$.

Using $\phi_\gamma$ allows us to transport our problem of modular composition to the right-hand side. For $a$ in $\mathbb{K}[x]_{<n}$ and $g$ in $\mathbb{K}[y]$, to find $g(a)$ rem $f$, this boils down to the following (see Algorithm 8.1):

- a forward change of basis: through inverse modular composition, compute $\alpha \in \mathbb{K}[y]_{<n}$ such that $a = \alpha(\gamma)$ rem $f$; this step also determines the minimal polynomial $\mu_\gamma$;

- a modular composition in the new basis: compute $\beta = g(\alpha)$ rem $\mu_\gamma$;

- a backward change of basis: the modular composition $\beta(\gamma)$ rem $f$, which equals $g(a)$ rem $f$.

*Computational aspects.* The second and third steps are modular compositions. They can performed efficiently by the approach of Section 1.2, by finding and using matrices of relations $R_m^{(\gamma,f)}$ and $R_m^{(\alpha,\mu_\gamma)}$, as long as certain genericity assumptions hold; this aspect is discussed below.

The first step, for the forward change of basis, is an instance of inverse modular composition and the calculation of a minimal polynomial. As mentioned in Section 1.1.2, Shoup's solutions recover both $\alpha$ and $\mu_\gamma$ from the power projections $(\ell(1), \ell(\gamma), \dots, \ell(\gamma^{2n-1}))$ and $(\ell(a), \ell(\gamma a), \dots, \ell(\gamma^{n-1}a))$, in the complexity of Brent and Kung's modular composition algorithm. Using matrices of relations we achieve a lower complexity, for a generic $\gamma$, as follows.

(1) *Matrix of relations and minimal polynomial.* Generalizing the power projections of $\gamma$, the algorithm of Section 1.2 computes the first $2d$ terms of $(X^T M_\gamma^k X)_k$, where $d = \lceil n/m \rceil$, and then reconstructs a basis $R_m^{(\gamma,f)}$ of $M_m^{(\gamma,f)}$ by solving Eq. (6) (with $\gamma$ instead of $a$). This basis
gives in particular the minimal polynomial $\mu_Y$, which appears as an entry of the Hermite normal form of this basis (Proposition 4.1).

(2) **Bivariate inverse composition.** The use of projections $(\ell(a), \ell(\gamma a), \dots, \ell(\gamma^{n-1}a))$ is directly generalized by computing the first $2d$ terms of $(X^T M_Y^k M_a \mathbf{1})_k$, where $\mathbf{1}$ is the first column of $X$, and solving

$$ (X^T(yI_n - M_Y)^{-1}X + O(y^{-2d-1})) v_{\tilde{\alpha}}(y) - v_N(y) = X^T(yI_n - M_Y)^{-1}M_a \mathbf{1} + O(y^{-d-1}) \quad (8) $$

for polynomial vectors $v_N$ and $v_{\tilde{\alpha}}$ in $\mathbb{K}[y]^m$ of degree less than $d$; the entries of the vector $v_{\tilde{\alpha}}$ are the coefficients of a bivariate polynomial $\tilde{\alpha}(x, y)$ of small degree such that $\tilde{\alpha}(x, y) \equiv a \mod f$.

As for Eq. (8), we actually work with an expansion at $y = 0$ rather than infinity.

(3) **Bivariate $\tilde{\alpha}$ to univariate $\alpha$.** The situation is now symmetric to that of the composition algorithm of Section 1.2: we consider again Eq. (5), where now $g$ is unknown (it is $\alpha$), $v$ is known (it is $\tilde{\alpha}(x, y)$) and both $R_m^{(\alpha, f)}$ and $v$ have degree at most $d$, so that the polynomial matrix problem can be solved in $\tilde{O}(m^\omega(d + n/m))$ operations.

This approach is detailed in Algorithm CHANGE OF BASIS, with the steps reordered and combined so as to retrieve both $R_m^{(\gamma, f)}$ and $v_{\tilde{\alpha}}$ from a single fraction reconstruction.

**Probabilistic aspects.** For a generic $\gamma$, one has $\deg(\mu_\gamma) = n$, so the isomorphism $\phi_\gamma$ is well defined. Using the Schwartz-Zippel lemma, it is straightforward to control the probability of having $\deg(\mu_\gamma) < n$.

For generic $\gamma$, we can then follow the approach described in Section 1.2 to perform the last step, modular composition by $\gamma$, with the desired complexity. The quantitative aspects can be worked out as well, and similar considerations hold for the first step, inverse modular composition by $\gamma$.

However, the composition in the second step, $g(\alpha)$ rem $\mu_\gamma$, is more delicate to analyze. We need the equality of modules $M_m^{(\alpha, \mu_\gamma)} = M_{m,m}^{(\alpha, \mu_\gamma)}$, and that a matrix of relations in this module can be reconstructed from the first $2\lceil v_m^{(\alpha, \mu_\gamma)}/m \rceil \le 2\lceil n/m \rceil$ elements of the corresponding matrix sequence; the analysis is made difficult by the fact that both $\alpha$ and $\mu_\gamma$ are nonlinear functions of the random element $\gamma$.

We prove that this happens for a generic $\gamma$ in two cases: when $f$ is separable in Section 8.3, and when $f$ is purely inseparable, with extra conditions, in Section 8.4; the latter case covers power series composition with $f = x^n$. In both situations, there is a nonzero polynomial $\Delta$ in $n$ variables such that the constraints above hold if $\Delta$ does not vanish at the coefficients of $\gamma$. We choose a random $\gamma$, and the probability of failure is again bounded by the Schwartz-Zippel lemma.

We do not have a proof that a generic $\gamma$ satisfies our requirements for an arbitrary $f$. Our algorithm for the general case proceeds by reduction to the two extreme cases above, separable and purely inseparable polynomials.

**From Monte Carlo to Las Vegas.** At this stage, we have a probabilistic algorithm of Monte Carlo type, that runs in the announced complexity and returns the correct result with a controlled probability of error. The next question is to modify the algorithm so that it detects and reports the unlucky choices of $\gamma$ for which its result would be incorrect.

In order to certify the result obtained for a random choice of $\gamma \in \mathbb{A}$, it would be sufficient to check the following properties:

(1) the computed matrix $R_m^{(\gamma, f)}$ is a basis of relations of $M_m^{(\gamma, f)}$;

(2) the minimal polynomial of $\gamma$ modulo $f$ has degree $n$;

(3) the computed matrix $R_m^{(\alpha, \mu_\gamma)}$ is a basis of relations of $M_m^{(\alpha, \mu_\gamma)}$.
However, we do not know how to check that all the columns of a matrix belong to the ideal $\langle f(x), y - \gamma(x) \rangle$ or $\langle \mu_{\gamma}(x), y - \alpha \rangle$ in sufficiently low complexity and in a deterministic way. The matrix $R_{m}^{(\gamma,f)}$ is easier to deal with: as it is expected to behave like in the generic case, its expected degree structure is known and the matrix can be certified by degree considerations (Item (ii) of Proposition 5.4, and Proposition 6.1). From there, the minimal polynomial of $\gamma$ can be computed efficiently via the Hermite normal form of $M_{m}^{(\gamma,f)}$, and it remains to check that it has degree $n$.

The other matrix, $R_{m}^{(\alpha,\mu_{\gamma})}$, carries more information about $a$ and cannot be expected to behave as predictably as $R_{m}^{(\gamma,f)}$. Our approach is to extract from its columns two small degree polynomials $r$ and $s$ in $\mathbb{K}[x, y]$. Since only two such polynomials are considered, they can be checked to vanish at $\alpha \mod \mu_{\gamma}$ by the Nüsken-Ziegler algorithm without affecting the asymptotic cost. Then, these two polynomials are used to construct a Sylvester matrix that can be used for composition instead of $R_{m}^{(\alpha,\mu_{\gamma})}$, without increasing the overall complexity (Algorithm 5.2).

*Note.* Equivalently, the randomization of our probabilistic algorithm can be seen as a change of projection. Indeed, let $P \in \mathbb{K}^{n \times n}$ have its jth column formed by the coefficients of $\gamma^{j-1}$ rem $f$. If $\gamma \mod f$ generates $\mathbb{K}[x]/\langle f \rangle$ and $M_{\alpha}$ is the matrix of multiplication by $\alpha \mod \mu_{\gamma}$ with basis $(1, y, \dots, y^{n-1})$, then the multiplications by $\alpha$ and by $a$ are related by

$$M_{\alpha} = P^{-1}M_{a}P. \qquad (9)$$

Hence

$$X^{\mathrm{T}} M_{\alpha}^{k} X = (X^{\mathrm{T}} P^{-1}) M_{a}^{k}(P X),$$

which, for instance on the right side, leads to considering the first $m$ columns of $P$ instead of $X$ for projecting. This amounts to kinds of structured projections $(V^{\mathrm{T}} M_{a}^{k} W)_{k \ge 0}$, i.e. with matrices $V$ and $W$ in a special proper subset of $\mathbb{K}^{n \times m}$.

## 1.4 Algorithm for the general case

The algorithm of Section 1.3 is proved to work when $f$ is either separable, or purely inseparable (for the latter, with extra conditions that are dealt with in Section 8.4). In Section 9, we address the general case, by first computing a separable decomposition of $f$ [53], yielding a factorization into a product into pairwise coprime terms of the form $h_i(x^{p^{e_i}})^{\ell_i}$, with $h_i$ separable and $e_i$, $\ell_i$ integers (here, $p$ is the characteristic of $\mathbb{K}$).

Working modulo each factor separately, we are thus left with the question of composition modulo a polynomial of the form $h(x^{p^e})^{\ell}$, with $h$ separable (all such results are eventually recombined via the Chinese remainder theorem).

For a modulus of the form $h(x)^{\ell}$, van der Hoeven and Lecerf showed how composition can be reduced to $\ell$ compositions modulo $h$, the computation of an annihilating polynomial modulo $h$, and a power series composition at precision $\ell$ with coefficients in $\mathbb{L} = \mathbb{K}[x]/\langle h(x) \rangle$ [30]. We extend this result to the case of moduli of the form $h(x^{p^e})^{\ell}$ in Section 9.4, involving essentially the same steps. The first two operations (compositions and annihilating polynomial modulo $h$) are directly handled by our results so far, but this is not quite the case for the latter, power series composition with coefficients in $\mathbb{L}$.

Our algorithms are written assuming they work over a field, as they perform zero-tests and inversions (compare this with Brent and Kung's algorithms, for instance, which apply over a ring). If $h$ is irreducible, $\mathbb{L}$ is a field, but if $h$ is only assumed to be separable, then $\mathbb{L}$ is only a product of fields. The *dynamic evaluation* paradigm [19] explains how an algorithm written for inputs lying in a field can carry over to inputs in a product of fields, but the original approach induces cost overheads that go beyond our cost target. Using van der Hoeven and Lecerf's efficient dynamic evaluation
strategy [33], we show how our algorithm for power series adapts to this situation (Section 9.2) without affecting the asymptotic runtime.

## 1.5 Previous algorithms in special cases

To compute $g(a)$ rem $f$, previous known improvements upon Brent and Kung's approach all have requirements on the input, either on some of the polynomials $f, g$, and $a$, or on the ring or field of coefficients — possibly with non-algebraic algorithms.

### 1.5.1 Special modulus.

**Power series.** For the special case $f = x^n$ of power series, Brent and Kung's second algorithm relying on Taylor expansion performs composition in only $\tilde{O}(n^{3/2})$ operations, provided $a'(0)$ and $([\sqrt{n \log(n)}])!$ are invertible in $\mathbb{A}$; the assumption on $a'(0)$ can be weakened [28, Sec. 3.4.3]. In more variables, even in the specific case $g(x, a)$ rem $f$ handled by the Nüsken-Ziegler algorithm, we do not know of any algorithm computing composition faster for power series than modulo arbitrary polynomials.

Faster composition in only $\tilde{O}(n)$ operations for $g(a)$ rem $x^n$ is possible for many special cases of $g$: when $g$ is a polynomial of degree $O(1)$, but also when it is a power series solution of a polynomial equation of degree $O(1)$ via Newton's iteration, or when it is a solution of a differential equation (e.g., exp), by first forming a differential equation for $g(a)$ and then solving it by Newton's iteration or other divide-and-conquer algorithms, generally in characteristic 0 or large enough [9, 15, 28, 54; 8, §13.4] Similarly, still in the case when $f = x^n$, if furthermore $a$ has specific properties, then composition of power series can be performed in $\tilde{O}(n)$ operations. This is the case when $a$ is a polynomial [15] of moderate degree (it is a part of Brent and Kung's fast composition algorithm), an algebraic power series [28], but also for a class of truncated power series that can be obtained via shifts, reversals, scalings, multiplications by polynomials, exponentials and logarithms [12].

**Separable polynomials.** Ritzmann observed that for a separable modulus $f(x) = (x - \epsilon_1) \cdots (x - \epsilon_n)$ with distinct $\epsilon_1, \ldots, \epsilon_n$ that are known, modular composition boils down to multipoint evaluation and interpolation [64], which can be computed in $\tilde{O}(n)$ arithmetic operations. When furthermore the ring of coefficients is $\mathbb{Z}$, he uses well-chosen $\epsilon_i$'s to give an efficient algorithm for composing power series, in a non-algebraic model of computation: if $g$ and $a$ over $\mathbb{Z}$ have coefficients bounded in absolute value by $K$, then $g(a)$ rem $x^n$ can be computed using $\tilde{O}(n^2 \log(K))$ bit operations, which is quasi-optimal since the output has bit size $\Omega(n^2 \log(K))$ in general.

**Chinese remainder theorem.** In our work, the cases of power series and of separable polynomials play an important role as well. We use the observation that if a factorization $f = f_1 \cdots f_s$ is known with the $f_i$'s relatively prime, then composing modulo $f$ reduces to composing modulo each $f_i$ and reconstructing the result via the Chinese remainder theorem. Several consequences of this observation have been discussed by van der Hoeven and Lecerf [31].

### 1.5.2 Special rings or fields.
For power series over a ring $\mathbb{A}$ of positive characteristic, Bernstein proposed an algebraic algorithm whose complexity is quasi-linear in $n$, with a constant factor that depends on the characteristic of the ring [6]. In particular, this algorithm is very efficient over rings whose characteristic is a product of small primes; if $\mathbb{A}$ is a ring of prime characteristic $p$ then the algorithm uses $\tilde{O}(np)$ operations in $\mathbb{A}$.

A further step forward was achieved by Umans in 2008 [68], with a new algorithm for modular composition modulo and arbitrary $f$, over finite fields of small characteristic: if $p$ is $n^{o(1)}$, his algorithm uses $n^{1+o(1)}$ base field operations. Later, Kedlaya and Umans introduced new techniques
for composition over finite rings of the form $(\mathbb{Z}/r\mathbb{Z})[z]/\langle h(z)\rangle$, for an integer $r$ and $h$ monic. For a finite field $\mathbb{K} = \mathbb{F}_q$, their algorithm runs in $n^{1+\epsilon} \log^{1+o(1)}(q)$ bit operations [48, Cor. 7.2].

As in Ritzmann's work, a key idea in [48, 68] is to exploit fast multipoint evaluation, but this time in a multivariate setting. The composition $g(a)$ rem $f$ is reduced to the evaluation at suitable points of a multivariate polynomial constructed from $g$ by an inverse Kronecker substitution, decreasing degrees at the expense of increasing the number of variables. Umans' algorithm performs the evaluation using the properties of the Frobenius endomorphism, while Kedlaya and Umans' proceeds by lifting to characteristic zero (which requires working in a bit complexity model). For arbitrary fields, efficient analogues of these multivariate multipoint evaluation algorithms are currently unknown.

## 1.6 Related questions

1.6.1 (Multivariate) multipoint evaluation. For simplicity, we limit the discussion to the case of a field; most of it extends to rings, with minor restrictions. The evaluation of a univariate $g \in \mathbb{K}[x]_{<n}$ at $n$ points in the field $\mathbb{K}$, and conversely the interpolation of a polynomial of degree $< n$ from $n$ values, are computable in quasi-linear complexity [22, Chap. 10]. For polynomials in at least two variables, however, the situation becomes tightly related to modular composition.

The motivation of Nüsken and Ziegler [59] was the evaluation of a polynomial $g \in \mathbb{K}[x, y]_{<(m,d)}$ at $n$ points $(x_k, y_k)_{1 \le k \le n}$ in general position, with $md = O(n)$. Their algorithm first computes a univariate interpolation polynomial such that $a(x_k) = y_k$ for all $k$; then the composition $b = g(x, a(x))$ rem $f$, where $f = \prod_k (x - x_k)$; and concludes by a univariate multipoint evaluation of $b$ at $x_1, \dots, x_n$. Since the univariate evaluation and interpolation are performed in essentially linear time, the complexity is dominated by the “uni-bivariate” modular composition $g(x, a)$ rem $f$.

The case when several points have the same $x$-value can be handled by an affine change of coordinates [59]; another approach, taken by Kedlaya and Umans, is to pick $n$ suitable points $t_1, \dots, t_n$ in $\mathbb{K}$, to compute two interpolation polynomials $a_x$ and $a_y$ in $\mathbb{K}[t]$, and thus reduce the evaluation to the fully bivariate modular composition $g(a_x, a_y)$ rem $f$, where now $f = \prod_k (t - t_k)$. This extends to an arbitrary number of variables and shows that multipoint evaluation in $s$ variables reduces to multivariate modular composition in the same number of variables [48, Thm. 3.3].

As mentioned in Section 1.5.2, Kedlaya and Umans actually make a heavy use of a converse reduction [48, Thm. 3.1]. If $g$ is a polynomial in $\mathbb{K}[x_1, \dots, x_s]$, the composition $g(a_1(x), \dots, a_s(x))$ rem $f$ reduces to a multipoint evaluation of a polynomial of smaller degree in each of its variables, whose number is increased. For the univariate case of composition ($s = 1$) studied here, the smallest possible number of variables for evaluation would be 2, leading to a bivariate evaluation of a polynomial of degree $\sqrt{n}$ at $\Theta(n^{3/2})$ points, which is too large for our target complexity. The next possible choice would be a polynomial of 3 variables in degree $n^{1/3}$ at $\Theta(n^{4/3})$ points. Unfortunately, we are not aware of a sufficiently efficient multipoint evaluation algorithm in 3 or more variables to make this approach succeed in the algebraic model.

1.6.2 Bivariate ideals. Viewing the problem of computing $g(a)$ modulo $f$ as a problem of reduction of $g$ modulo the ideal $I = \langle y - a(x), f(x) \rangle$, we introduce bivariate polynomials in a way different from the inverse Kronecker substitution mentioned above. Gröbner bases are commonly used for reductions modulo multivariate ideals. A division with remainder similar to that in Eq. (5) would be achieved via reduction by an appropriate Gröbner basis of $I$, provided we could compute this basis and perform the reduction in good complexity. However, already the size of the Gröbner basis itself may be $\Theta(n^{3/2})$ (see the example below), hence exceed our target complexity.

For an ideal given by two generic bivariate polynomials of degree $n$ (hence the ideal is of degree $n^2$) and the graded lexicographic order, van der Hoeven and Larrieu avoid the use of an
explicit Gröbner basis. They show that a *concise representation* of the basis of size only $\tilde{O}(n^2)$ is sufficient for reducing a polynomial modulo the ideal in time $\tilde{O}(n^2)$ [29]; the concise representation consists in particular of truncations of well chosen polynomials in the ideal. It is unclear to us whether a similar truncation strategy could be applied specifically to $\mathcal{I}$, whose degree is only $n$. Instead, the matrices of relations we compute give a set of small degree polynomials in $\mathcal{I}$ that may not generate the whole ideal (see Section 4.1), but provide a process of complexity $\tilde{O}(n^\kappa)$ for the reduction modulo $\mathcal{I}$ of Eq. (5). These polynomials generate the same ideal as the first polynomials in the Gröbner basis of $\mathcal{I}$ for the lexicographic order (see Corollary 4.3).

The concise representation of Gröbner bases has also been exploited by van der Hoeven and Lecerf for computing the minimal polynomial of the multiplication by $y$ modulo $\mathcal{J}$, when $\mathcal{J} = \langle f_1, f_2 \rangle$ is generated by two generic polynomials $f_1, f_2 \in \mathbb{K}[x, y]$ and $\mathbb{K}$ is a finite field [36, Sec. 4]. They apply the transposition principle to a bivariate modular composition map modulo $\mathcal{J}$, then compute the minimal polynomial from the resulting bivariate power projections [41, Sec. 6]. The evaluation of the composition map modulo $\mathcal{J}$ is again in $\tilde{O}(n^2)$, thanks to the concise representation [36]. In our case of $\mathcal{I} = \langle y - a(x), f(x) \rangle$ and for a generic $a$, matrices of relations allow us to compute the minimal polynomial of the multiplication by $y$ modulo $\mathcal{I}$ in complexity $\tilde{O}(n^\kappa)$ (see Section 10.1); matrices of relations are obtained via a bivariate power projection process that can be regarded, in part, as dual to Nüsken and Ziegler’s bivariate modular composition algorithm (Section 3.4.3).

*Note.* For a sufficiently large field $\mathbb{K}$, take $f = (x-1)\cdots(x-n)$, where $n = k(k+1)/2$, and $a \in \mathbb{K}[x]$ the polynomial of degree smaller than $n$ such that $a(i) = [\sqrt{2i}]$ for $1 \le i \le n$. Then the reduced Gröbner bases for the graded lexicographic order and for the lexicographic order, both with $y < x$, coincide. They contain one polynomial with leading term $x^i y^{k-i}$ for each $i \in \{0, \dots, k\}$. Counting the number of monomials of these polynomials shows that this basis has $k(k+1)(k+2)/3 + (k+1)$ monomials; this is of the order of $n^{3/2}$.

### 1.6.3 Modular composition and multipoint evaluation with precomputation.
Quasi-linear modular composition $g(a)$ rem $f$ is feasible after precomputations on $(f, a)$ only, for $a$ generic and $f$ square-free [57].

Likewise, after precomputations on the evaluation points and under genericity assumptions on them, quasi-linear multivariate multipoint evaluation is feasible [35], as well as quasi-linear bivariate interpolation [57]. Furthermore, for bivariate evaluation, genericity can be replaced by randomization [34].

In these works, the precomputation stages are at least as expensive as the fastest known corresponding modular composition or multipoint evaluation algorithms. They have a feature common with our composition algorithm: from $f, a$ (or from the evaluation points), they compute a set of polynomials which belong to $\langle y - a(x), f(x) \rangle$ (or vanish at the points), and allow for efficient degree reduction of the polynomial to compose with (or to evaluate). This set is either akin to several matrices of relations of $M_m^{(a,f)}$ for a small number of values of $m$ ranging from 1 to $n$ [57], or is a collection of well-chosen polynomials in several Gröbner bases for subsets of the points so as to build a multivariate divide and conquer evaluation tree [34, 35].

## 1.7 Outline

Section 2 introduces some notation and our computational model. Section 3 details baby steps/giant steps techniques used in our composition algorithm: known ones such as in Brent and Kung's composition, and new ones such as for computing truncated powers which give access to $(X^\top M_a^k X)_k$. Section 4 studies matrices of relations and how they are used in our composition algorithm, whereas Section 5 shows how to compute them efficiently by matrix fraction reconstruction under some
assumptions on $(f, a, m)$. Section 6 presents an algorithm for the change of basis of Section 1.3: it finds the minimal polynomial $\mu_Y$ and an inverse composition $\alpha$ such that $\alpha(\gamma) \equiv a \mod f$, under assumptions on $(f, \gamma, m)$. Section 7 studies these assumptions, and in particular gives precise generic situations where they hold. Section 8 describes our main randomized composition algorithm and proves its correctness for a separable $f$ and for a purely inseparable $f$ (generalizing $f = x^n$); then Section 9 handles the general case of composition modulo any $f$. Finally, in Section 10, we state resulting complexity improvements for several variants of modular composition and other related problems.

## 2 PRELIMINARIES

*Notation.* In this article, $\mathbb{K}$ is an arbitrary field. For bivariate polynomials in variables $x$ and $y$, $\deg_X$ and $\deg_Y$ give the degree in $x$ and in $y$. For any polynomial or power series $p = \sum_i p_i x^i$, we use the following notation for a “slice” of it: $[p]_j^k = p_j + p_{j+1}x + \cdots + p_{j+k}x^k$. The ideal generated by polynomials $f_1, \ldots, f_k$ in an ambient ring (which will be clear from the context) is denoted by $\langle f_1, \ldots, f_k \rangle$.

Vectors, such as elements of $\mathbb{K}^m$ or $\mathbb{K}[y]^m$, are seen as column vectors by default; when row vectors are considered this is explicit in our notation, e.g. $\mathbb{K}^{1 \times m}$ or $\mathbb{K}[y]^{1 \times m}$. We often identify a polynomial $g_0(y) + \cdots + g_{m-1}(y)x^{m-1}$ in $\mathbb{K}[x, y]_{<(m, \cdot)}$ with the column vector $(g_0 \cdots g_{m-1})^\top$ in $\mathbb{K}[y]^m$ of its coefficients on the basis $(1, x, \ldots, x^{m-1})$ of the $\mathbb{K}[y]$-module $\mathbb{K}[x, y]_{<(m, \cdot)}$.

For $a$ and $f$ in $\mathbb{K}[x]$, $M_a$ denotes the matrix of the linear map of multiplication by $a$ in $\mathbb{K}[x]/\langle f(x) \rangle$ with basis $(1, x, \ldots, x^{n-1})$, and $\mu_a$, resp. $\chi_a$, denotes the minimal, resp. characteristic polynomial of $a$ in $\mathbb{K}[x]/\langle f(x) \rangle$ (that is, the minimal and characteristic polynomials of $M_a$).

Whenever the context is sufficiently clear, particularly in Sections 4, 5 and 7, notation such as $M_m^{(a,f)}$, $v_m^{(a,f)}$ defined in the introduction is shortened into $M_m$, $v_m$. We keep the superscripts in important statements.

*Computational model.* Our algorithms are written in pseudo-code, using standard syntax elements (for-loops, if statements, ...). Informally, we count all arithmetic operations $\{+, -, \times, \div\}$ and zero-tests in $\mathbb{K}$ at unit cost. The underlying complexity model is the *computation tree* [17, Sec. 4.4].

A computation tree over $\mathbb{K}$ is a binary tree whose nodes are partitioned into *input nodes* that form an initial segment of the tree starting at the root, *computation nodes* with outdegree 1, *branching nodes* with outdegree 2 and *output nodes* at the leaves. To each node is associated a label. Computation nodes are labelled by constants in $\mathbb{K}$ or operations in $\{+, -, \times, \div\}$, in which case they also carry references to two previous input or computation nodes; branching nodes are labelled by zero-tests, referring to some previously computed quantity. Each leaf $v$ is labelled with a sequence of references $(u_1, \ldots, u_{\ell(v)})$ to previous input or computation nodes. The *cost* of a computation tree is simply its height $\tau$, that is, the maximum length of a path from the root to a leaf.

It then makes sense to *evaluate* a computation tree at an element of $\mathbb{K}^s$ — called *input to the tree*, where $s$ is the number of input nodes, following a path from the root to a leaf. After the input nodes, the path is constructed as follows. Each computation node is assigned a value derived from the label it carries, when it is defined. Otherwise, e.g., in case of a division by 0, the path stops. At a branching node the path branches left or right depending on whether the value it refers to vanishes or not. At a leaf $v$ with label $(u_1, \ldots, u_{\ell(v)})$, the output of the computation is the tuple of the values computed at nodes $u_1, \ldots, u_{\ell(v)}$. In that case, the computation tree is called *evaluable* at the input. Overall, the computation requires at most $\tau$ arithmetic operations in $\mathbb{K}$. An algorithm is called *quasi-linear* when the height of its computation tree is linear (up to logarithmic factors) in the number of inputs. It is called *quasi-optimal* when this height is linear (up to logarithmic factors) in the number of inputs *plus* the maximum number of values returned by the output nodes.
A computation tree takes inputs of fixed length. In order to solve a problem for inputs of arbitrary size, we need a *family* of trees, parametrized by the input size. Every algorithm we describe using pseudo-code in this article, and all algorithms that we rely on from the literature, can be described by a family of computation trees.

The translation from pseudo-code to computation tree is usually rather direct, and as is customary in the literature, we do not do it explicitly. In a nutshell, for-loops and recursive calls are “unrolled”; if statements that test whether a computed quantity vanishes yield branching nodes, etc. Some operations in our pseudo-code may not be directly available in our model (as we only allow arithmetic operations in $K$ and zero-test), but they can be rewritten in a way that complies with our requirements. This is for instance the case when we compute the degree of a polynomial (as in Euclid’s GCD algorithm): this can be achieved by scanning its coefficients, in order of decreasing degree, until a nonzero one is found.

We allow our algorithms to return flags (such as `FAIL`, or `CERT/NOCERT`). This can be done in this model, by returning constants in the vector of outputs, such as 1 for `FAIL` and 0 otherwise.

Finally, several of our algorithms rely on randomization; however, we do not want to introduce another arithmetic operation for the selection of random field elements. One reason for this is that we invoke a result by van der Hoeven and Lecerf [33] on the transformation of computation trees for *directed evaluation* in Section 9.2, and that result is explicitly written in a deterministic model. Instead, “random” field elements are given to our procedures as extra input parameters.

# 3 SIMULTANEOUS MODULAR OPERATIONS BY MATRIX MULTIPLICATION

A key ingredient in fast modular composition algorithms is to turn the problem into the simultaneous evaluation of polynomials of smaller degree, and exploit the structure brought by this simultaneity using matrix multiplication. In this section, after reviewing Brent and Kung’s original algorithm and giving a direct extension of it, we use this idea in two further contexts: Nüsken and Ziegler’s bivariate modular composition algorithm, and the computation of truncations of powers of the form $a^k \text{ rem } f$. Both arise in our algorithms, and are bottlenecks in their complexity.

## 3.1 Brent and Kung's algorithm

3.1.1 *Modular composition.* We start with a review of Brent and Kung's algorithm to compute $g(a) \text{ rem } f$, pointing out the impact of rectangular matrix multiplication [37] and how the runtime depends on the degrees of both $f$ and $g$ [65, Fact 3.1]. This can be seen as an introduction to the Nüsken-Ziegler algorithm, which generalizes this approach to a bivariate $g$.

**PROPOSITION 3.1.** Given polynomials $f \in K[x]$ of degree $n$, $a$ in $K[x]_{<n}$ and $g$ in $K[y]_{<d}$, Algorithm *MODULARCOMPOSITION-BRENTKUNG* computes $g(a) \text{ rem } f$ using $\tilde{O}((1+n/d)d^{\omega_2/2})$ operations in $K$.

*PROOF.* Correctness follows from noticing that at Step 7, $b_i \equiv g_{ir} + g_{ir+1}a + \dots + g_{ir+r-1}a^{r-1} \pmod f$ holds for all $i$, where $g_j$ is the coefficient of degree $j$ in $g$ for all $j$. The cost of the algorithm comes from $\Theta(d^{1/2})$ multiplications modulo $f$, which use $\tilde{O}(nd^{1/2})$ operations in $K$, and a matrix product in sizes $s \times r$ and $r \times n$, with both $s$ and $r$ in $\Theta(d^{1/2})$. This product can be done through $\lfloor n/d \rfloor \le n/d+1$ matrix products in sizes $s \times r$ and $r \times d$, each of which takes $O(d^{\omega_2/2})$ operations in $K$. $\square$

*Note.* In the analysis, dividing the matrix product into blocks, as we did, is sub-optimal. Using rectangular matrix multiplication directly, the runtime can be described by the finer estimate $\tilde{O}(d^{\omega_2 \log(n)/\log(d)/2})$. Here, the notation $\omega_\theta$ is a feasible exponent for rectangular matrix multiplication for any real number $\theta$: there is an algorithm which multiplies an $n \times [n^\theta]$ matrix by an $[n^\theta] \times n$ matrix using $O(n^{\omega_\theta})$ operations [52]. However, this refinement complicates notation, and would
<caption>**Algorithm 3.1 MODULARCOMPOSITION-BRENTKUNG(f, a, g)**</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[y]_{<d}$

Output: $g(a)$ rem $f$

$$1: r \leftarrow [d^{1/2}], s \leftarrow [d/r]$$

$$2: \hat{a}_0 \leftarrow 1$$

$$3: \textbf{for } i = 1, \dots, r \textbf{ do } \hat{a}_i \leftarrow a \cdot \hat{a}_{i-1} \textbf{ rem } f \qquad \triangleright \hat{a}_i = a^i \text{ rem } f$$

$$4: A \leftarrow \text{matrix } (\text{coeff}(\hat{a}_i, j))_{\substack{0 \le i < r \\ 0 \le j < n}} \text{ in } \mathbb{K}^{r \times n} \qquad \triangleright \text{ coefficient of degree } j \text{ of } \hat{a}_i$$

$$5: G \leftarrow \text{matrix } (\text{coeff}(g, ir + j))_{\substack{0 \le i < s \\ 0 \le j < r}} \text{ in } \mathbb{K}^{s \times r}$$

$$6: B = (b_{i,j})_{\substack{0 \le i < s \\ 0 \le j < n}} \leftarrow \text{GA in } \mathbb{K}^{s \times n}$$

$$7: \textbf{for } i = 0, \dots, s-1 \textbf{ do } b_i \leftarrow b_{i,0} + \dots + b_{i,n-1}x^{n-1}$$

$$8: \textbf{return } b_0 + b_1\hat{a}_r + \dots + b_{s-1}\hat{a}_r^{s-1} \text{ rem } f \qquad \triangleright \text{ Horner evaluation}$$

not be of use for our main results. The same remark holds for several other runtime estimates in this section, such as Lemmas 3.3 and 3.5.

### 3.1.2 Power projection.
The transposition principle implies the existence of an algorithm POWER-PROJECTION with the same asymptotic runtime as Algorithm MODULARCOMPOSITION-BRENTKUNG and with the following signature [65].

<caption>**Algorithm 3.2 POWERPROJECTION(f, a, d, (r_i)_{0 \le i < n})**</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a$ in $\mathbb{K}[x]_{<n}$, $d$ in $\mathbb{N}$, $(r_i)_{0 \le i < n}$ in $\mathbb{K}^n$
Output: $(\ell(1), \ell(a), \dots, \ell(a^{d-1} \bmod f))$, with $\ell(b_0 + \dots + b_{n-1}x^{n-1}) = r_0b_0 + \dots + r_{n-1}b_{n-1}$

Whereas seeing the details of Algorithm MODULARCOMPOSITION-BRENTKUNG is useful as a preamble to the Nüsken-Ziegler algorithm, Algorithm POWERPROJECTION only plays the role of a subroutine in one other algorithm given just below. Moreover, giving its pseudo-code would require us to introduce concepts such as transposed product, that would not used any further in this text. We refer the reader to [67], which gives all details but uses classical matrix arithmetic (with $\omega_2 = 4$), so the runtime of that version is $\tilde{O}(d^2 + nd)$ instead of $\tilde{O}((1 + n/d)d^{\omega_2/2})$.

### 3.1.3 Small minimal polynomial.
Modular composition can be sped up when the minimal polynomial $\mu_a$ of a modulo $f$ has degree at most $d$, for some (small) integer $d \le n$. To compute $g(a)$ rem $f$, the idea is that once $\mu_a$ is known, $\tilde{g} = g$ rem $\mu_a$ can be computed, and then $g(a) \equiv \tilde{g}(a) \mod f$. The computation of the latter by Proposition 3.1 benefits from $\tilde{g}$ having degree less than $d$.

It remains to discuss how to compute $\mu_a$. Here, we follow an algorithm of Shoup's (the deterministic version, for $f$ irreducible, is in [65, Thm. 4]; the randomized one is in [66, Sec. 4]). We take a random linear form $\ell: \mathbb{K}[x]/\langle f \rangle \to \mathbb{K}$ and compute the sequence $(\ell(a^k \mod f))_{0 \le k < 2d}$. With high probability, its minimal polynomial is $\mu_a$; the algorithm verifies whether it is the case, and returns either a correct result or FAIL. In Algorithm MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL, $\mu_a$ is computed using an Extended Euclidean scheme called MINIMALPOLYNOMIALFORSEQUENCE [22, Algo. 12.9].

The following lemma analyses the runtime of this procedure, and the probability of success. As per our convention at the end of Section 2, the “random” linear form $\ell$ is actually given as an argument, through the vector $(r_0, \dots, r_{n-1}) \in \mathbb{K}^n$ of its coefficients.

**LEMMA 3.2.** *Given $f \in \mathbb{K}[x]$ of degree $n$, $a \in \mathbb{K}[x]_{<n}$, $g \in \mathbb{K}[y]_{<n}$, $d$ in $\{1, \dots, n\}$ and $(r_i)_{0 \le i < n}$ in $\mathbb{K}^n$, Algorithm MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL uses $\tilde{O}(nd^{(\omega_2/2)^{-1}})$ operations*
**Algorithm 3.3 MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL(f, a, g, d, (rᵢ)₀≤i<n)**

Input: $f$ of degree $n$ in $K[x]$, $a$ in $K[x]_{<n}$, $g$ in $K[y]_{<n}$, $d$ in $\{1, \dots, n\}$, $(r_i)_{0\le i<n}$ in $K^n$

Output: $g(a)$ rem $f$ or FAIL

1: $(v_0, \dots, v_{2d-1}) \leftarrow \text{POWERPROJECTION}(f, a, 2d, (r_i)_{0\le i<n})$

2: $\mu \leftarrow \text{MINIMALPOLYNOMIALFORSEQUENCE}(v_0, \dots, v_{2d-1})$

3: $t \leftarrow \text{MODULARCOMPOSITION-BRENTKUNG}(f, a, \mu)$

4: **if** $t \neq 0$ **then** **return** FAIL

5: **else** **return** MODULARCOMPOSITION-BRENTKUNG($f, a, g$ rem $\mu$)

in $K$ and returns either $g(a)$ rem $f$ or FAIL. If $\mu_a$ has degree at most $d$, and the entries of $(r_i)_{0\le i<n}$ are chosen uniformly and independently from a finite subset $S$ of $K$, then with probability at least $1-n/\text{card}(S)$ the algorithm returns $g(a)$ rem $f$ and computes $\mu_a$ as a by-product. If $\mu_a$ has degree more than $d$, the algorithm returns FAIL.

PROOF. For any given $a$ in $K[x]_{<n}$ and $(r_0, \dots, r_{n-1})$, the algorithm computes a polynomial $\mu$ and tests whether $\mu(a) \equiv 0 \bmod f$; if it is the case, it reduces $g$ modulo $\mu$ before doing a modular composition. Hence, the output may only be $g(a)$ rem $f$ or FAIL, as claimed; it is FAIL if and only if the value $t$ at Step 3 does not vanish.

By the discussion in Section 3.1.2 and Proposition 3.1, the call to POWERPROJECTION takes $\tilde{O}((1+n/d)d^{\omega_2/2})$ operations in $K$; because we take $d \le n$, this is $\tilde{O}(nd^{(\omega_2/2)^{-1}})$. Step 2 then computes a nonzero annihilating polynomial of degree at most $d$ in $\tilde{O}(d)$ operations in $K$ [22, Algo. 12.9]. The remaining lines call Algorithm MODULARCOMPOSITION-BRENTKUNG with a last argument of degree at most $d$, so the cost is $\tilde{O}(nd^{(\omega_2/2)^{-1}})$ again; this establishes the claim on the runtime.

Suppose first that $\mu_a$ has degree greater than $d$. Then since $\deg(\mu) \le d$, $\mu(a)$ rem $f$ cannot be zero, so the output is FAIL, as claimed.

Finally, suppose that the minimal polynomial $\mu_a$ has degree at most $d$ and that the entries of $(r_i)_{0\le i<n}$ are chosen uniformly at random and independently from a set $S$ in $K$. With $M_a$ the multiplication matrix by $a \bmod f$ and 1 the vector $(1, 0, \dots, 0)$, the sequence $((r_i)^T M_a^k 1)_{k \ge 0}$ is $(\ell(a^k \mod f))_{k \ge 0}$ and the sequence $(M_a^k 1)_{k \ge 0}$ is $(a^k \text{rem } f)_{k \ge 0}$. Following the probabilistic analysis of Wiedemann's algorithm [42, Lem. 2; 43, Lem. 1], the probability that their minimal polynomials coincide is at least $1-n/\text{card}(S)$. When this occurs, Step 2 computes $\mu_a$; the value $t$ at Step 3 is then zero, and the output is $g(a)$ rem $f$. $\square$

The main idea in this algorithm — computing an annihilating polynomial for $a$ and using it to reduce $g$ — is actually at the core of our main algorithm as well. Key differences are that we compute several annihilating polynomials (which we call relations), and use them to reduce $g$ into a bivariate polynomial. We then apply Nüsken and Ziegler's extension of Brent and Kung's algorithm, which we present now.

## 3.2 Bivariate composition

Here we describe the Nüsken-Ziegler algorithm for modular composition [59], which computes $g(x, a)$ rem $f$ for a bivariate $g$ in $K[x, y]$.

First, however, we address the following question: given $f$ of degree $n$ in $K[x]$, $a$ in $K[x]_{<n}$ and an $s$-uple $(g_0, \dots, g_{s-1})$ in $K[x, y]^{s}_{<(m,r)}$, compute all compositions

$$ (g_0(x, a) \text{ rem } f, \dots, g_{s-1}(x, a) \text{ rem } f) \in K[x]^s. $$

The solution designed by Nüsken and Ziegler [59] boils down to a multiplication of polynomial matrices. Writing the polynomials $g_i$ as the rows of their coefficients in $y$ gives an $s \times r$ matrix $G$
whose entries are polynomials in $K[x]_{<m}$. Writing the powers of $1, a, \dots, a^{r-1}$ rem $f$ in a column vector $A$ reduces the simultaneous composition to a matrix-vector product $GA$. This is turned into a matrix-matrix product by spreading the coefficients of $A$ as follows. If

$$g_i(x, y) = \sum_{0 \le j < r} g_{i,j}(x)y^j,$$

then computing the product

$$B = \begin{pmatrix} g_{0,0}(x) & \cdots & g_{0,r-1}(x) \\ \vdots & & \vdots \\ g_{s-1,0}(x) & \cdots & g_{s-1,r-1}(x) \end{pmatrix} \begin{pmatrix} [a^0 \text{ rem } f]_0^{m-1} & \cdots & [a^0 \text{ rem } f]_{\lceil \frac{n}{m} \rceil -1m}^{m-1} \\ \vdots & & \vdots \\ [a^{r-1} \text{ rem } f]_0^{m-1} & \cdots & [a^{r-1} \text{ rem } f]_{\lceil \frac{n}{m} \rceil -1m}^{m-1} \end{pmatrix}$$

yields a matrix whose entry $B_{i,l}$ is

$$B_{i,l} = \sum_{0 \le j < r} g_{i,j}(x) [a^j \text{ rem } f]_{lm}^{m-1}.$$

Summing the $B_{i,l}x^{\ell m}$ modulo $f$, for $\ell = 0, \dots, \lceil n/m \rceil - 1$, then provides $g_i(x, a)$ rem $f$ for $i = 0, \dots, s - 1$ at low cost. This is detailed in Algorithm 3.4 and Lemma 3.3.

<caption>**Algorithm 3.4 SIMULTANEOUSBIVARIATEMODULARCOMPOSITION(f, a, g₀, ..., gₖ₋₁, m, r)**</caption>

Input: $f$ of degree $n$ in $K[x]$, $a$ in $K[x]_{<n}$, $(g₀, \dots, gₖ₋₁)$ in $K[x, y]^{s}_{<(m,r)}$

Output: $(g₀(x, a) \text{ rem } f, \dots, gₖ₋₁(x, a) \text{ rem } f)$

1: $\hat{a}_0 \leftarrow 1$

2: **for** $i = 1, \dots, r$ **do** $\hat{a}_i \leftarrow a \cdot \hat{a}_{i-1} \text{ rem } f$

3: $A \leftarrow \text{matrix} ([\hat{a}_i]_{jm}^{m-1})_{\substack{0 \le i < r \\ 0 \le j < \lceil \frac{n}{m} \rceil}} \text{ in } K[x]_m^{r \times \lceil \frac{n}{m} \rceil}$

4: $G \leftarrow \text{matrix} (g_{i,j}(x))_{\substack{1 \le i \le s \\ 0 \le j < r}} \text{ in } K[x]_m^{s \times r}$, where $g_i(x, y) = \sum_j g_{i,j}(x)y^j$

5: $B = (B_{i,j})_{\substack{0 \le i < r \\ 0 \le j < \lceil \frac{n}{m} \rceil}} \leftarrow \text{GA}$

6: **for** $i = 0, \dots, s-1$ **do** $b_i \leftarrow (\sum_{0 \le j < \lceil \frac{n}{m} \rceil} B_{i,j} x^{jm}) \text{ rem } f$

7: **return** $(b₀, \dots, bₖ₋₁)$

**LEMMA 3.3** ([59, LEM. 10(iii)]). *Algorithm SIMULTANEOUSBIVARIATEMODULARCOMPOSITION computes* $(g₀(x, a) \text{ rem } f, \dots, gₖ₋₁(x, a) \text{ rem } f)$. *Assuming* $s \in \tilde{O}(r)$, *it uses* $\tilde{O}(c(n, m, r²)) = \tilde{O}((m + n/r²)r^{\omega₂})$ *operations in K, with* $c(\cdot)$ *from Eq. (2).*

**PROOF.** Steps 1 and 2 use $\tilde{O}(rn)$ operations. Similarly, for each $i = 0, \dots, s-1$, Step 6 uses $\lceil n/m \rceil$ additions in $O(m)$ operations each and one reduction in $\tilde{O}(m + n)$ operations. The total cost of Step 6 is thus $\tilde{O}(s(n + m))$.

Steps 3 and 4 do not use any arithmetic operation. The most expensive step is Step 5, the product of an $s \times r$ matrix by an $r \times \lceil n/m \rceil$ matrix, both with entries in $K[x]_{<m}$. Using the same kind of block decomposition as in Proposition 3.1, this is done using $\lceil \lceil n/m \rceil/r² \rceil \in n/(mr²) + O(1)$ products in sizes $s \times r$ and $r \times r²$. With the assumption $s \in \tilde{O}(r)$, each of them uses $\tilde{O}(mr^{\omega₂})$ operations in $K$, for a total of $\tilde{O}(c(n, m, r²)) = \tilde{O}((m + n/r²)r^{\omega₂})$ operations in $K$.

The other steps, in $\tilde{O}((r+s)(n+m)) = \tilde{O}(rm+rn)$, are at most of the same order, since $\omega₂ \ge 3$. □

Algorithm SIMULTANEOUSBIVARIATEMODULARCOMPOSITION is the central step in bivariate composition as showed in Algorithm BIVARIATEMODULARCOMPOSITION, leading to the complexity stated in Proposition 3.4.
**Algorithm 3.5 BIVARIATEMODULARCOMPOSITION(f, a, g)** (Nüsken-Ziegler algorithm [59])

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[x, y]_{<(m,d)}$

Output: $g(x, a)$ rem $f$

1: $r \leftarrow \lfloor d^{1/2} \rfloor$, $s \leftarrow \lfloor d/r \rfloor$

2: Write $g(x, y) = g_0(x, y) + g_1(x, y)y^r + \dots + g_{s-1}(x, y)y^{r(s-1)}$ with $\deg_y(g_i) < r$ for $0 \le i < s$

3: $(b_0, \dots, b_{s-1}) \leftarrow \text{SIMULTANEOUSBIVARIATEMODULARCOMPOSITION}(\hat{f}, a, g_0, \dots, g_{s-1}, m, r)$

▷ $b_i = g_i(x, a)$ rem $f$

4: $\hat{a} \leftarrow a^r$ rem $f$

▷ is computed in the previous step

5: **return** $b_0 + b_1\hat{a} + \dots + b_{s-1}\hat{a}^{s-1}$

▷ Horner evaluation

**PROPOSITION 3.4** ([59, THM. 9]). Given $f \in \mathbb{K}[x]$ of degree $n$, $a \in \mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[x, y]_{<(m,d)}$, Algorithm BIVARIATEMODULARCOMPOSITION computes $g(x, a)$ rem $f$ using $\tilde{O}(c(n, m, d)) = \tilde{O}((m + n/d)d^{\omega_2/2})$ operations in $\mathbb{K}$, with $c(\cdot)$ from Eq. (2).

**PROOF.** The correctness of the algorithm is straightforward. For the complexity analysis, we first note that $s \sim r \sim d^{1/2}$. Lemma 3.3 then shows that the complexity of Step 3 is $\tilde{O}(c(n, m, r^2)) = \tilde{O}(c(n, m, d))$. The other task involving arithmetic operations is the final Horner evaluation which costs $\tilde{O}(rn)$. As in the proof of Lemma 3.3, this is smaller than the other part, since $\omega_2 \ge 3$. $\square$

## 3.3 Sequence of truncated modular powers

Another key ingredient in our composition algorithm also relies on polynomial matrix multiplication. To our knowledge this is a new algorithm, whose properties are summarized in the next lemma.

**LEMMA 3.5.** Given $f$ of degree $n$ in $\mathbb{K}[x]$, $(p_0, \dots, p_{r-1})$ in $\mathbb{K}[x]_{<n}$, $(q_0, \dots, q_{s-1})$ in $\mathbb{K}[x]_{<n}^s$ and $m \in \mathbb{N}_{>0}$, Algorithm SIMULTANEOUSTRUNCATEDMODULARMULTIPLICATION computes the simultaneous truncated modular multiplications

$$ \{[p_i q_j \text{ rem } f]_0^{m-1} \mid 0 \le i < r, 0 \le j < s\}. $$

If $s \in \tilde{O}(r)$, it uses $\tilde{O}(c(n, m, r^2)) = \tilde{O}((m + n/r^2)r^{\omega_2})$ operations in $\mathbb{K}$, with $c(\cdot)$ from Eq. (2).

The basic approach to this problem is to first compute all the products $p_i q_j$ modulo $f$ and then truncate the computed polynomials. However, this produces an intermediate result of size $nrs$, which is $\Theta(nr^2)$ when $s$ is in $\Theta(r)$, and is larger than our target complexity.

Hereafter, we use the reversal of a polynomial $p \in \mathbb{K}[x]$ with respect to $m \in \mathbb{N}$ defined by $\text{rev}(p, m) = x^m p(1/x)$; when $m = \text{deg}(p)$ this is the classical reciprocal of the polynomial $p$.

**PROOF.** For all $i < r$ and $j < s$, let $p_i q_j = h_{i,j} f + r_{i,j}$, with $\deg(r_{i,j}) < n$, be the Euclidean division of the product $p_i q_j$ by the polynomial $f$. The main task of the algorithm is to compute the truncated quotients $[h_{i,j}]_0^{m-1}$ (Steps 1 to 5); from there, the truncated remainders $[r_{i,j}]_0^{m-1}$ are easily obtained (Step 6) at a total cost of $\tilde{O}(mrs)$ operations.

For the efficient computation of the quotients $h_{i,j}$, we rely on the classical approach via reciprocals and power series operations. More specifically we use the identity

$$ \text{rev}(h_{i,j}, n-2) = \frac{\text{rev}(p_i, n-1) \text{rev}(q_j, n-1)}{\text{rev}(f, n)} \text{rem } x^{n-1} = \bar{p}_i \bar{q}_j \text{rem } x^{n-1}, $$

obtained by evaluating $p_i q_j = h_{i,j} f + r_{i,j}$ at $1/x$ and multiplying by $x^{n-2}/f(1/x) = x^{2n-2}/\text{rev}(f, n)$; here we have $\bar{p}_i = \text{rev}(p_i, n-1)$ and $\bar{q}_j = \text{rev}(q_j, n-1)/\text{rev}(f, n)$ rem $x^{n-1}$, as in the pseudo-code. The idea of our algorithm is to compute only the last $m$ coefficients of this expansion by means of two polynomial matrix multiplications.
<caption>**Algorithm 3.6 SIMULTANEOUS TRUNCATED MODULAR MULTIPLICATION** (*f*, (*p*<sub>*i*</sub>)<sub>*i*<sub>r</sub></sub>, (*q*<sub>*j*</sub>)<sub>*j*<sub>s</sub></sub>, *m*)</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $(p_0, \dots, p_{r-1})$ in $\mathbb{K}[x]_{<n}^r$, $(q_0, \dots, q_{s-1})$ in $\mathbb{K}[x]_{<n}^s$, $m \in \mathbb{N}_{>0}$
Output: $([p_i q_j \text{ rem } f]_0^{m-1})_{\substack{0 \le i < r \\ 0 \le j < s}}$

1. $(\ell, t) \leftarrow$ (quotient, remainder) in the Euclidean division $n - m - 1 = \ell m + t$ with $\ell = 0$ if $m \ge n$

2. **for** $i = 0, \dots, r-1$ **do** $\bar{p}_i \leftarrow \text{rev}(p_i, n-1)$

3. **for** $j = 0, \dots, s-1$ **do** $\bar{q}_j \leftarrow$ power series expansion $\text{rev}(q_j, n-1)/\text{rev}(f, n)$ rem $x^{n-1}$

4. Form the matrices

$$P_1 \leftarrow ([\bar{p}_i]_{j_{m+t}}^{m-1})_{\substack{0 \le i < r \\ 0 \le j \le l}} \in \mathbb{K}[x]_m^{r \times (l+1)} \quad P_2 \leftarrow ([\bar{p}_i]_{j_{m+t}}^{m-1})_{\substack{0 \le i < r \\ 0 \le j < l}} \in \mathbb{K}[x]_m^{r \times l}$$

and

$$Q_1 \leftarrow ([\bar{q}_j]_{(\ell-i)m}^{m-1})_{\substack{0 \le i \le l \\ 0 \le j < s}} \in \mathbb{K}[x]_m^{(\ell+1) \times s} \quad Q_2 \leftarrow ([\bar{q}_j]_{(\ell-1-i)m}^{m-1})_{\substack{0 \le i < l \\ 0 \le j < s}} \in \mathbb{K}[x]_m^{l \times s}$$

5. $H \leftarrow [P_1 Q_1]_0^{m-1} + [P_2 Q_2]_m^{m-1} + ([\bar{p}_i]_0^{t-1} [\bar{q}_j]_{\ell m+1}^{m-1+t-1})_{\substack{t-1 \\ 0 \le i < r \\ 0 \le j < s}} \triangleright H = (\bar{h}_{i,j})_{i,j}$ is in $\mathbb{K}[x]_m^{r \times s}$

6. **for** $i = 0, \dots, r-1$ and $j = 0, \dots, s-1$ **do** $r_{i,j} \leftarrow (p_i q_j - \text{rev}(\bar{h}_{i,j}, m-1)f)$ rem $x^m$

7. **return** $(r_{i,j})_{\substack{0 \le i < r \\ 0 \le j < s}}$

For any $t \in \{0, \dots, m-1\}$, for any polynomials $a, b$ written as

$$a = [a]_0^{t-1} + x^t \sum_{i \ge 0} a_i x^{im} \quad \text{with } \deg(a_i) < m, \quad b = \sum_{j \ge 0} b_j x^{jm} \quad \text{with } \deg(b_j) < m,$$

and for any positive integer $\ell$, one has

$$[ab]_{\ell m + t}^{m-1} = \left[ \sum_{i+j=\ell} a_i b_j \right]_0^{m-1} + \left[ \sum_{i+j=\ell-1} a_i b_j \right]_m^{m-1} + [[a]_0^{t-1} [b]_{\ell m+1}^{m-1+t-1}]_0^{m-1} \quad (10)$$

(The last summand is a product of small degree polynomials that is 0 when $t = 0$.) We use this formula with $\ell$ and $t$ as defined in Step 1, so that the left-hand side is $[ab]_{n-m-1}^{m-1}$; applying this to $a = \bar{p}_i$ and $b = \bar{q}_j$ gives $\bar{h}_{i,j} = [\text{rev}(h_{i,j}, n-2)]_{n-m-1}^{m-1}$, and thus $[h_{i,j}]_0^{m-1}$ by reversal.

Since $\ell \sim n/m$, using this formula for a single pair $i, j$ requires $\tilde{O}(n)$ operations in $\mathbb{K}$ and thus is as costly as computing $\bar{p}_i \bar{q}_j$ rem $x^n$. In our algorithm the gain comes from using this formula simultaneously for several products, in which case matrix multiplication helps.

The first multiplication in Step 5 is the matrix product

$$\begin{pmatrix} [\bar{p}_0]_t^{m-1} & \cdots & [\bar{p}_0]_{n-m-1}^{m-1} \\ \vdots & \ddots & \vdots \\ [\bar{p}_{r-1}]_t^{m-1} & \cdots & [\bar{p}_{r-1}]_{n-m-1}^{m-1} \end{pmatrix} \begin{pmatrix} [\bar{q}_0]_{\ell m}^{m-1} & \cdots & [\bar{q}_{s-1}]_{\ell m}^{m-1} \\ \vdots & \ddots & \vdots \\ [\bar{q}_0]_0^{m-1} & \cdots & [\bar{q}_{s-1}]_0^{m-1} \end{pmatrix}$$

Its entries are the first summand in Eq. (10) for $a = \bar{p}_i$ and $b = \bar{q}_j$, for $0 \le i < r$ and $0 \le j < s$. Similarly, the second summand in Eq. (10) is obtained from the matrix product

$$\begin{pmatrix} [\bar{p}_0]_t^{m-1} & \cdots & [\bar{p}_0]_{n-2(m-1)}^{m-1} \\ \vdots & \ddots & \vdots \\ [\bar{p}_{r-1}]_t^{m-1} & \cdots & [\bar{p}_{r-1}]_{n-2(m-1)}^{m-1} \end{pmatrix} \begin{pmatrix} [\bar{q}_0]_{(\ell-1)m}^{m-1} & \cdots & [\bar{q}_{s-1}]_{(\ell-1)m}^{m-1} \\ \vdots & \ddots & \vdots \\ [\bar{q}_0]_0^{m-1} & \cdots & [\bar{q}_{s-1}]_0^{m-1} \end{pmatrix}$$

In terms of complexity, the multiplication $P_1 Q_1$ involves $r \times (\ell + 1)$ and $(\ell + 1) \times s$ matrices, while $P_2 Q_2$ involves $r \times \ell$ and $\ell \times s$ matrices; all four operands have degree less than $m$.
Since $\ell = \lfloor(n-1)/m\rfloor - 1$, we have $\ell + 1 \le n/m$, so each matrix product can be done using at most $\lceil n/(mr^2) \rceil \le n/(mr^2) + 1$ products in sizes $r \times r^2$ and $r^2 \times s$. Since $s \in \tilde{O}(r)$, each of these take $\tilde{O}(mr^{\omega_2})$, for a total cost of $c(n, m, r^2) = \tilde{O}((m + n/r^2)r^{\omega_2})$.

The other operations performed by the algorithm are $O(r)$ power series expansions at precision $n-1$ in $\tilde{O}(n)$ operations each (precisely, one inverse and $t$ multiplications, see Step 3), and $O(r^2)$ power series expansions at precision $m$ in $\tilde{O}(m)$ operations each (precisely, at most $3rs$ multiplications and $rs$ subtractions, see Steps 5 and 6). This amounts to a total of $\tilde{O}(nr + mr^2)$ operations, and can thus be neglected, since $\omega_2 \ge 3$. $\square$

Using simultaneous truncated modular multiplication combined with a baby steps/giant steps strategy leads to Algorithm TRUNCATEDPOWERS, with the following properties.

**PROPOSITION 3.6.** Given $f$ in $K[x]$ of degree $n$, $a$ and $b$ in $K[x]_{<n}$, $m$ and $d$ in $N_{>0}$, Algorithm TRUNCATEDPOWERS computes the truncations

$$[a^k b \text{ rem } f]_0^{m-1}, \quad 0 \le k < d$$

using $\tilde{O}(c(n, m, d)) = \tilde{O}((m + n/d)d^{\omega_2/2})$ operations in $K$, with $c(\cdot)$ from Eq. (2).

**PROOF.** The algorithm computes $1, a, \dots, a^{r-1}$ rem $f$ and $b, ba^r, \dots, ba^{(s-1)r}$ rem $f$, which costs $\tilde{O}(nr)$ operations in $K$ since $r \sim s$. From these two sets of polynomials, Algorithm SIMULTANEOUS-TRUNCATEDMODULARMULTIPLICATION is then used to compute $[ba^k \text{ rem } f]_0^{m-1}$ for $0 \le k \le rs - 1$ using $\tilde{O}((m+n/r^2)r^{\omega_2})$ operations, by Lemma 3.5; since $\omega_2 \ge 3$, this is larger than $\tilde{O}(nr)$. The choice of $s$ makes $(s-1)r < d \le rs$, so the output consists of the terms $k = i + rj$ for $j < s - 1$ and $i < r$, and for $j = s - 1$ and $i < d - (s-1)r \in \{1, \dots, r\}$. $\square$

**Algorithm 3.7 TRUNCATEDPOWERS($f, a, b, m, d$)**

**Input:** $f$ of degree $n$ in $K[x]$, $a$ and $b$ in $K[x]_{<n}$, $m$ and $d$ in $N_{>0}$

**Output:** the truncated powers $[ba^k \text{ rem } f]_0^{m-1}$ for $0 \le k < d$

1. $r \leftarrow \lfloor d^{1/2} \rfloor$; $s \leftarrow \lfloor d/r \rfloor$

2. $\hat{a}_0 \leftarrow 1$; **for** $i = 1, \dots, r$ **do** $\hat{a}_i \leftarrow a \cdot \hat{a}_{i-1}$ **rem** $f$ ▶ $\hat{a}_i = a^i \text{ rem } f$

3. $\bar{a}_0 \leftarrow b$; **for** $j = 1, \dots, s-1$ **do** $\bar{a}_j \leftarrow \hat{a}_r \cdot \bar{a}_{j-1}$ **rem** $f$ ▶ $\bar{a}_j = ba^{jr}$ **rem** $f$

4. $(c_{i,j})_{\substack{0 \le i < r \\ 0 \le j < s}} \leftarrow \text{SIMULTANEOUSTRUNCATEDMODULARMULTIPLICATION}(f, \hat{a}_0, \dots, \hat{a}_{r-1}, \bar{a}_0, \dots, \bar{a}_{s-1}, m)$

5. **for** $i = 0, \dots, r-1$ and $j = 0, \dots, s-2$ **do** $r_{i+rj} \leftarrow c_{i,j}$

for $i = 0, \dots, d-1-(s-1)r$ **do** $r_{i+r(s-1)} \leftarrow c_{i,s-1}$

6. **return** $(r_k)_{0 \le k < d}$

Finally, Algorithm BLOCKTRUNCATEDPOWERS computes truncations of products of the form $x^i a^k \text{ rem } f$, which are needed in our composition algorithm; here, we assume that $f(0)$ is nonzero (see Remark 3.8).

**PROPOSITION 3.7.** Given $f$ in $K[x]$ of degree $n$ with $f(0) \ne 0$, $a$ in $K[x]_{<n}$, $m$ and $d$ in $N_{>0}$, Algorithm BLOCKTRUNCATEDPOWERS computes

$$[x^i a^k \text{ rem } f]_0^{m-1}, \quad 0 \le i < m, 0 \le k < d$$

using $\tilde{O}(c(n, m, d)) + O(m^2d) = \tilde{O}((m + n/d)d^{\omega_2/2}) + O(m^2d)$ operations in $K$, with $c(\cdot)$ from Eq. (2).
<caption>**Algorithm 3.8 BLOCKTRUNCATEDPOWERS(f, a, m, d)**</caption>

**Input:** $f$ of degree $n$ in $\mathbb{K}[x]$, with $f_0 = f(0) \neq 0$, $a$ in $\mathbb{K}[x]_{<n}$, $m$ and $d$ in $\mathbb{N}_{>0}$
**Output:** the truncated powers $[x^i a^k \text{ rem } f]_0^{m-1}$, for $0 \le i < m$ and $0 \le k < d$

1: $(r_k)_{0 \le k < d} \leftarrow \text{TRUNCATEDPOWERS}(f, a, x^{m-1}, 2m-1, d)$

2: $f_n \leftarrow \text{coeff}(f, n)$

3: **for** $k = 0, \dots, d-1$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp; $a_{m-1,k} \leftarrow r_k$

5: &nbsp;&nbsp;&nbsp;&nbsp; **for** $i = m-1, \dots, 2, 1$ **do**

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $c \leftarrow a_{i,k}(0)f_n/f_0$

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $a_{i-1,k} \leftarrow (a_{i,k} + c[f]_0^{m+i-1})/x$

8: **return** $([a_{i,k}]_0^{m-1})_{\substack{0 \le i < m \\ 0 \le k < d}}$

$$ \triangleright r_k = [x^{m-1} a^k]_0^{2m-2} $$

$\triangleright$ leading coefficient

$$ \triangleright a_{i-1,k} = [x^{i-1} a^k \text{ rem } f]_0^{m+i-2} $$

**PROOF.** Proposition 3.6 shows that the first step computes the sequence $[x^{m-1}a^k]_0^{2m-2}$ for $k = 0, \dots, d-1$ in the announced complexity. The remaining truncations are obtained from the identity

$$ [x p \text{ rem } f]_0^j = x[p \text{ rem } f]_0^{j-1} - \frac{p_{n-1}}{f_n} [f]_0^j, $$

for any integer $j$ and polynomial $p$, where $p_{n-1}$ is the coefficient of degree $n-1$ of $p \text{ rem } f$ and $f_n$ is the coefficient of degree $n$ in $f$. If we know $[x p \text{ rem } f]_0^j$, we get $p_{n-1}f_0/f_n$ as its constant coefficient, whence $p_{n-1}$ since $f_0 \neq 0$ and from there $[p \text{ rem } f]_0^{j-1}$ is easily obtained. At iteration $k$ of the loop at Step 3, the truncation $[x^{m-1}a^k \text{ rem } f]_0^{2m-2}$ computed previously is used to deduce all $[x^{m-1-i}a^k \text{ rem } f]_0^{2m-2-i}$ for $1 \le i < m$ in $O(m^2)$ operations. Thus this loop has a total cost of $O(m^2d)$ operations. $\square$

**REMARK 3.8.** The assumption $f(0) \neq 0$ is harmless in the context of modular composition: in the computation of $g(a) \text{ rem } f$, one can rather evaluate $g(y)$ at $a(x+c)$ modulo $f(x+c)$ for a randomly chosen $c \in \mathbb{K}$, and unshift the result. See Steps 3 and 12 in Algorithm MODULARCOMPOSITIONBASECASE.

## 3.4 Notes

### 3.4.1 Linear algebra interpretation.
Representing polynomials by their vector of coefficients leads to viewing the operations performed by Algorithms BIVARIATEMODULARCOMPOSITION and TRUNCATEDPOWERS as computing the product of special matrices by column vectors. Recall the notation $M_a$ for the $n \times n$ matrix of multiplication by $a$ mod $f$ in the basis $(1, x, \dots, x^{n-1})$, and $X$ for the matrix $(I_m)^T \in \mathbb{K}^{n \times m}$ with $m \in \{1, \dots, n\}$. Then Algorithms BIVARIATEMODULARCOMPOSITION and TRUNCATEDPOWERS correspond respectively to multiplication by

$$ K_{m,d}^{(a,f)} = (X \quad \dots \quad M_a^{d-1}X) \in \mathbb{K}^{n \times (md)} \quad \text{and} \quad L_{m,d}^{(a,f)} = \begin{pmatrix} X^T \\ \vdots \\ X^T M_a^{d-1} \end{pmatrix} \in \mathbb{K}^{(md) \times n}. \quad (11) $$

Indeed, $K_{m,d}^{(a,f)}$ is the matrix of the mapping $\kappa_{m,d}^{(a,f)}$ of bivariate modular composition with bounded degrees, as computed by Algorithm BIVARIATEMODULARCOMPOSITION:

$$ \kappa_{m,d}^{(a,f)} : \mathbb{K}[x, y]_{< (m,d)} \to \mathbb{K}[x]_{< n} \\ g(x, y) \mapsto g(x, a) \text{ rem } f. $$
On the other hand, $L_{m,d}^{(a,f)}$ represents the mapping $\lambda_{m,d}^{(a,f)}$ that extracts the low-degree part of multiplications by powers of $a$, as computed by Algorithm TRUNCATEDPOWERS:

$$ \lambda_{m,d}^{(a,f)} : \mathbb{K}[x]_{<n} \to \mathbb{K}[x]_{<m}^{d} \\ b \mapsto ([b \ \mathrm{rem} \ f]_{0}^{m-1}, \dots, [ba^{d-1} \ \mathrm{rem} \ f]_{0}^{m-1}). $$

These maps and matrices play an important role in the study of the generic behaviour of our algorithm starting from Section 7.

### 3.4.2 Complexity equivalence.
Proposition 3.4 (Algorithm BIVARIATEMODULARCOMPOSITION) and Proposition 3.6 (Algorithm TRUNCATEDPOWERS) give similar complexity bounds for the evaluation of $\kappa_{m,d}^{(a,f)}$ and $\lambda_{m,d}^{(a,f)}$, but the computational equivalence of these problems, possibly up to some conditions, is still unclear to us in general.

However, for $m = 1$, when $f(0) \neq 0$, these two problems are indeed equivalent. This is a consequence of the transposition principle in an indirect way, starting from the equality

$$ L_{1,n}^{(a,f)} v_b = (M_b K_{1,n}^{(a,f)})^{\mathrm{T}} \mathbf{1}, $$

where $v_b$ is the vector associated to $b$, $M_b$ is the matrix of multiplication by $b \mod f$, and $\mathbf{1}$ is the first canonical vector. First, this equality gives a way to evaluate $\lambda_{1,n}^{(a,f)}$ for the cost of one multiplication by $M_b^{\mathrm{T}}$ (i.e., $\tilde{O}(n)$ by the transposition principle), plus one multiplication by the transpose of $K_{1,n}^{(a,f)}$, which has the same asymptotic cost as that of $K_{1,n}^{(a,f)}$ itself, by the same principle. Conversely, if $v = M_b^{\mathrm{T}} \mathbf{1}$, the equality reads $(K_{1,n}^{(a,f)})^{\mathrm{T}} v = L_{1,n}^{(a,f)} v_b$, so that, again by the transposition principle, the evaluation of $\kappa_{1,n}^{(a,f)}$ reduces to that of $\lambda_{1,n}^{(a,f)}$ provided $v_b$ can be computed from $v$ in low complexity. When $f(0) \neq 0$, this can be done in $\tilde{O}(n)$ by solving a linear system of Hankel type [67, Sec. 3].

If $f(0) = 0$, it is unclear whether such a reduction holds: in the special case $f = x^n$, the map $\lambda_{1,n}^{(a,f)}$ becomes much simpler, as it simply computes the sequence $a_0^j b_0$ for $0 \le i < d$, where $a_0$ and $b_0$ are the constant coefficients of $a$ and $b$. This only requires a linear number $O(d)$ of operations. On the other hand, $\kappa_{1,d}$ is the composition of a univariate polynomial $g(y)$ of degree less than $d$ with the power series $a(x)$ and no quasi-linear complexity result is known for this operation.

### 3.4.3 Transposition of the Nüsken-Ziegler algorithm.
Finally, we discuss a different approach to Algorithm BLOCKTRUNCATEDPOWERS, that actually bypasses Algorithm TRUNCATEDPOWERS altogether, and uses the transpose of Algorithm BIVARIATEMODULARCOMPOSITION instead.

Algorithm BLOCKTRUNCATEDPOWERS computes the $m \times m$ projections $H_k = X^\mathrm{T} M_{ak} X$, for $k = 0, \dots, d-1$, using the fact that for $k < d$, $H_k$ can be deduced in $O(m^2)$ operations (for-loop at Step 5) from the column vector $\tilde{X}^\mathrm{T} M_{ak} u$ of size $2m-1$, where $\tilde{X} = (I_{2m-1} \ 0)^\mathrm{T} \in \mathbb{K}^{(2m-1) \times n}$ and $u$ is the $m$th column of $X$, i.e. the $m$th canonical vector (Step 1). (Here we have taken $m \le (n+1)/2$.) Algorithm TRUNCATEDPOWERS computes the vectors $\tilde{X}^\mathrm{T} M_a^k u$ for $0 \le k < d$ using $\tilde{O}(c(n, m, d))$ operations.

Alternatively, we can consider a recursion similar to the one in the proof of Proposition 3.7, but now for learning a new coefficient of a polynomial rather than a coefficient of a new polynomial. Assuming $f(0) \neq 0$, for a polynomial $p$ one has

$$ [xp \ \mathrm{rem} \ f]_{i+1}^{0} = [p \ \mathrm{rem} \ f]_{i}^{0} + (c/f_{0})[f]_{i+1}^{0}, $$

where $c$ is the coefficient of degree 0 of $xp \ \mathrm{rem} \ f$: we see that from the row vector $1^\mathrm{T} M_{x-m+1} a_k \tilde{X}$, one can also deduce $H_k = X^\mathrm{T} M_{ak} X$ using $O(m^2)$ operations.
Now, if we set $v = M_{x-m+1}^T 1$, computing $v^T K_{d,2m-1}^{(a,f)}$ precisely gives all vectors $1^T M_{x-m+1}^T a k \bar{X}$, for $0 \le k < d$. Since $v$ can be computed in quasi-linear time, the application of the transposition principle to Algorithm BIVARIATEMODULARCOMPOSITION shows that these vectors can be computed using $\tilde{O}(c(n, m, d))$ operations. Altogether, this gives an alternative to Algorithm BLOCKTRUNCATEDPOWERS with the same asymptotic complexity.

# 4 MATRICES OF RELATIONS FOR COMPOSITION

The heart of our algorithm for finding $g(a)$ rem $f$ is the computation of a *matrix of relations*, which gives a collection of polynomials of small degree in the ideal $\mathcal{I}$ generated by $y - a$ and $f$ in $\mathbb{K}[x, y]$. For a given positive integer $m$, these polynomials are in the $\mathbb{K}[y]$-module $M_m^{(a,f)}$ obtained by degree restriction as $\mathcal{I} \cap \mathbb{K}[x, y]_{< (m, \cdot)}$.

In Section 4.1 we show that the invariant factors of $M_m^{(a,f)}$ are the $m$ invariant factors of highest degree of the characteristic matrix $yI_n - M_a$, where $M_a$ is the matrix of multiplication by $a$ mod $f$. Once a matrix of relations has been obtained, it can be used to perform composition by reducing univariate composition to a small bivariate composition problem; this is described in Section 4.2. Finally, in Section 4.3, the results of this section are applied to the efficient computation of annihilating polynomials for a modulo $f$.

In all of Section 4, notation such as $M_m^{(a,f)}$ and $v_m^{(a,f)}$ is shortened into $M_m$ and $v_m$, except for the main definitions and statements, as there is no ambiguity as to the dependency on $a$ or $f$.

## 4.1 Structure of the module of relations

This section introduces the module of relations $M_m^{(a,f)}$ and relates it to the characteristic matrix.

### 4.1.1 Definitions.

*Relations.* We call *relations* the polynomials of the ideal $\mathcal{I} = \langle y - a(x), f(x) \rangle$ of $\mathbb{K}[x, y]$; these are the bivariate polynomials $r(x, y)$ such that $r(x, a) \equiv 0 \bmod f$, i.e., they are algebraic relations satisfied by $a \bmod f$. We are interested in those relations whose $x$-degree is bounded from above by a given positive integer $m$. They form the $\mathbb{K}[y]$-module

$$ \mathcal{M}_m^{(a,f)} = \{ r(x, y) \in \mathbb{K}[x, y]_{< (m, \cdot)} \mid r(x, a(x)) \equiv 0 \bmod f \} = \mathcal{I} \cap \mathbb{K}[x, y]_{< (m, \cdot)}, $$

which is denoted $M_m$ when $a$ and $f$ are clear from the context.

This is a $\mathbb{K}[y]$-submodule of $\mathbb{K}[x, y]_{< (m, \cdot)}$, itself a free $\mathbb{K}[y]$-module with basis $(1, x, \dots, x^{m-1})$. As stated in Section 2, we often identify a polynomial $r_0(y) + \dots + r_{m-1}(y)x^{m-1}$ in $\mathbb{K}[x, y]_{< (m, \cdot)}$ with the column vector $(r_0 \cdots r_{m-1})^T$ in $\mathbb{K}[y]^m$ of its coefficients on that basis. Since $\mathbb{K}[y]$ is a principal ideal domain, $M_m$ is free as well, and it has rank $m$ since it contains $\mu_a \mathbb{K}[y]^m$, where $\mu_a$ is the minimal polynomial of $a \bmod f$.

In terms of ideals, there is a chain of inclusions $\{0\} = \langle M_0 \rangle \subseteq \dots \subseteq \langle M_{n+1} \rangle = \mathcal{I}$; the latter identity follows from the fact that $y - a$ and $f$ have $x$-degree less than $n+1$. Furthermore $M_1 \neq \{0\}$ since $\mu_a$ belongs to $\mathcal{I} \cap \mathbb{K}[y]$. For small $m$, the module $M_m$ may not contain all the information in $\mathcal{I}$: the inclusion $\langle M_m \rangle \subset \mathcal{I}$ can be strict.

*Matrix and basis of relations, determinantal degree.* A matrix of relations of $M_m$ is any nonsingular matrix in $\mathbb{K}[y]^{m \times m}$ whose columns are elements of the module $M_m$ (represented as column vectors). Such a matrix is further called a *basis of relations* if its columns generate $M_m$; all bases of relations of $M_m$ can be obtained from any single one of them via right multiplication by a unimodular matrix in $\mathbb{K}[y]^{m \times m}$, i.e., a matrix whose determinant is in $\mathbb{K} \setminus \{0\}$. It follows that any matrix of relations of $M_m$ is a square, nonsingular right multiple of any basis of relations of $M_m$, and therefore bases of
this sequence belong to $\mathcal{M}_m$, and considering their $x$-degrees shows that they are $\mathbb{K}[y]$-linearly independent.

Any polynomial $p(x, y) \in \mathcal{M}_m$ is a $\mathbb{K}[y]$-linear combination of the first $m$ of these polynomials. Indeed, it can be divided by the Gröbner basis with a remainder equal to 0; in view of its degree in $x$, only these $m$ polynomials are involved in the division. This proves the claim on the basis of $\mathcal{M}_m$ described in Eq. (12).

The matrix $T(y) \in \mathbb{K}[y]^{m \times m}$ representing this basis (with basis elements written in columns) is upper triangular, with its first min($m, s$) diagonal entries being the first min($m, s$) polynomials in Eq. (13) in this order, and with its remaining diagonal entries being nonzero elements of $\mathbb{K}$. Furthermore each of these diagonal entries divides all other entries in the same column, hence the Smith normal form of $T(y)$ has the same diagonal entries as $T(y)$, which proves the claim on the invariant factors of $\mathcal{M}_m$. $\square$

**PROOF OF PROPOSITION 4.1.** Corollary 4.3 implies that the determinantal degree $v_m$ of $\mathcal{M}_m$ is the sum of the degrees of the elements of the first min($m, s$) elements of Eq. (13). It follows that

$$v_m \le \delta_k \deg(r_0) + \dots + \delta_1 \deg(r_{k-1}) = n,$$

where the last identity comes from considering the $\mathbb{K}$-vector space dimension of $\mathbb{K}[x, y]/\mathcal{I}$. If $s \le m$, all the nontrivial invariant factors appear and the bound is reached, while otherwise $m < s$ and $\deg \det(B)$, being the sum of the degrees of $m$ nonconstant polynomials, is at least $m$.

If $B$ is a basis of $\mathcal{M}_m$, then there exists a unimodular matrix $U \in \mathbb{K}[y]^{m \times m}$ such that $UB = T$ with $T$ as in the previous proof. If moreover $B$ is upper triangular, then so is $U$ and since $\det(U) \in \mathbb{K} \setminus \{0\}$, the diagonal entries of $U$ belong to $\mathbb{K} \setminus \{0\}$. It follows that $B$ has the same diagonal entries as $T$ up to multiplication by nonzero elements of $\mathbb{K}$.

The columns of the characteristic matrix $yI_n - M_a$ represent the polynomials $x^k(y-a(x))$ for $0 \le k < n$, making this matrix a matrix of relations of $\mathcal{M}_n$. It has determinantal degree $\deg(\chi_a) = n$, which coincides with the determinantal degree of $\mathcal{M}_n$, by the previous inequalities. Thus $yI_n - M_a$ is actually a basis of $\mathcal{M}_n$ and its invariant factors are given by the previous paragraph. $\square$

**4.1.4 Note.** For $m \in \{1, \dots, n\}$, the module of relations $\mathcal{M}_m$ is isomorphic to the module of vector generators for the matrix sequence $\{M_a^k X\}_{k \ge 0}$, where $X = (I_m \ 0)^T \in \mathbb{K}^{n \times m}$ as above (this elementary fact is established within the proof of Lemma 5.2, for instance); the bases of relations are the minimal generating polynomials for that sequence [46, 70].

The relation between Coppersmith's block Wiedemann algorithm and invariant factors of a characteristic matrix was described by Kaltofen and Villard: they show that for generic projections $V$ and $W$ in $\mathbb{K}^{n \times \ell}$ and $\mathbb{K}^{n \times m}$, with $\ell \ge m$, the invariant factors of minimal generating polynomial of the sequence $(V^T A^k W)_{k \ge 0}$ are the $m$ invariant factors of largest degree of the characteristic matrix $yI_n - A$ [46, Thm. 2.12]. In our more specific setting, Proposition 4.1 shows that this relation holds when the right projection is the structured matrix $X$ (see also Section 5.1.4).

## 4.2 Composition using matrices of relations

Matrices of relations are used to reduce the univariate problem $g(a)$ for $f$ with $g \in \mathbb{K}[y]$, to a bivariate one with better degree properties, thanks to a matrix division.

**4.2.1 Division for polynomial matrices.** If $R$ is a nonsingular matrix in $\mathbb{K}[y]^{m \times m}$ and $v_g$ is a vector in $\mathbb{K}[y]^m$, then there exist quotient and remainder vectors $w$ and $v_{\tilde{g}}$ such that

$$v_g = Rw + v_{\tilde{g}}, \quad (14)$$

and each entry of $v_{\tilde{g}}$ has degree less than that of the corresponding row of $R$ [39, Thm. 6.3-15, p. 389]. The latter reference actually states a stronger condition on $v_g$, namely that the matrix
fraction $R^{-1}v_{\tilde{g}}$ is strictly proper (see Section 5.1.1); this implies the above degree condition [39, Lem. 6.3-10, p. 383], which is sufficient for our needs.

For computing this division, it is customary to use $\mathbb{K}[y]$-linear system solving. For this, we rely on a kernel basis algorithm [75]: this returns $v$ in $\mathbb{K}[y]^m$ and $r$ in $\mathbb{K}[y]$ such that $R^{-1}v_g = v/r$, with $r$ of minimal degree. From this the remainder is obtained as $v_{\tilde{g}} = R(v \text{ rem } r)/r$, and here we do not need the quotient vector $w$.

### 4.2.2 Composition Algorithm.
In the case where $v_g = (g \ 0 \ \cdots \ 0)^T$ and $R$ is a matrix of relations of $M_m$ of degree at most $d$, the remainder in the above division is a vector $v_{\tilde{g}}$ of degree less than $d$ whose entries yield $\tilde{g} \in \mathbb{K}[x, y]_{< (m, d)}$ such that $g - \tilde{g} \in M_m$. Thus, analogously to a reduction modulo a Gröbner basis of the ideal $\mathcal{I} = \langle y - a, f \rangle$, this provides a bivariate polynomial $\tilde{g}$ with smaller degree in $y$ and controlled degree in $x$, and such that $\tilde{g} - g \in \mathcal{I}$, that is, $\tilde{g}(x, a) \equiv g(a) \mod f$.

Algorithm **BIVARIATEMODULARCOMPOSITIONWITHRELATION MATRIX** is given a matrix of relations $R$ of $M_m$ as a parameter and performs this division; then it completes the composition by evaluating $\tilde{g}(x, a) \text{ rem } f$ using Algorithm **BIVARIATEMODULARCOMPOSITION**. Algorithm **BIVARIATEMODULARCOMPOSITIONWITHRELATION MATRIX** actually accepts a slightly more general input: $g$ can be a bivariate polynomial with $x$-degree less than $m$ (however, the rest of the article focuses on the case of $g$ in $\mathbb{K}[y]$ highlighted above). The algorithm accepts $g$ of arbitrary degree in $y$, but the cost analysis is done under the assumption $\text{deg}_y(g) \in O(n)$.

<caption>**Algorithm 4.1 BIVARIATEMODULARCOMPOSITIONWITHRELATION MATRIX(f, a, g, R)**</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[x, y]_{< (m, .)}$,

$R \in \mathbb{K}[y]_{\le d}^{m \times m}$ a matrix of relations of $M_m^{(a,f)}$

Output: $q(x, a) \text{ rem } f$

1. Write $g(x, y) = g_0(y) + g_1(y)x + \cdots + g_{m-1}(y)x^{m-1}$ and set $v_g \leftarrow (g_0 \cdots g_{m-1})^T \in \mathbb{K}[y]^m$

2. ▶ Compute $v \in \mathbb{K}[y]^m$ and $r \in \mathbb{K}[y]$ using [75, Algo. 1]

$$ \binom{v}{r} \in \mathbb{K}[y]^{m+1} \leftarrow \text{MINIMALNULLSPACEBASIS}((R - v_g), (d, \dots, d, \text{deg}_y(g))) $$

3. $v_{\tilde{g}} \leftarrow R(v \text{ rem } r)/r \in \mathbb{K}[y]_{<d}^m$ ▶ $v \text{ rem } r$ is the vector of entry-wise remainders

4. $\tilde{g}(x, y) \leftarrow$ the polynomial in $\mathbb{K}[x, y]_{< (m, d)}$ corresponding to $v_{\tilde{g}}$

5. **return** BIVARIATEMODULARCOMPOSITION($f, a, \tilde{g}$) ▶ $\tilde{g}(x, a) \text{ rem } f$, Algorithm 3.5

**PROPOSITION 4.4.** Given $f$ in $\mathbb{K}[x]$ of degree $n$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[x, y]_{< (m, .)}$ with $\text{deg}_y(g) = O(n)$ and a matrix of relations $R$ in $\mathbb{K}[y]_{\le d}^{m \times m}$ of $M_m^{(a,f)}$, Algorithm **BIVARIATEMODULARCOMPOSITIONWITHRELATION MATRIX** computes $g(x, a) \text{ rem } f$ using $\tilde{O}(m^\omega(d+n/m) + c(n, m, d))$ operations in $\mathbb{K}$, with $c(\cdot)$ from Eq. (2).

**PROOF.** First, Step 2 computes $r \in \mathbb{K}[y]$ and $v = rR^{-1}v_g \in \mathbb{K}[y]^m$ with $r$ of minimal degree. Indeed, since $R$ is nonsingular, the right kernel of $(R - v_g) \in \mathbb{K}[y]^{m \times (m+1)}$ has rank 1. We use [75, Algo. 1] to compute a basis $(v^\top r)^T$ of this kernel. Thus by construction $Rv = rv_g$ holds, and the fact that $(v^\top r)^T$ generates the kernel ensures that the greatest common divisor of $r$ and all the entries of $v$ is 1, hence the minimality of $\text{deg}(v)$ and $\text{deg}(r)$.

At Step 3 one considers the vector $\tilde{v} = v \text{ rem } r \in \mathbb{K}[y]^m$ such that $\text{deg}(\tilde{v}) < \text{deg}(r)$ and $v = rw+\tilde{v}$ for some $w \in \mathbb{K}[y]^m$. It follows that $v_g = Rv/r = Rw + v_{\tilde{g}}$, where $v_{\tilde{g}} = R\tilde{v}/r$ is the vector computed at Step 3; by construction the $i$th entry of $v_{\tilde{g}}$ has degree less than that of the $i$th row of $R$. In short, Steps 2 and 3 compute a vector $v_{\tilde{g}} \in \mathbb{K}[y]^m$ which has degree less than $d$ and is a remainder of $v_g$ modulo $R$. Since $R$ is a matrix of relations, the polynomial $\tilde{g}(x, y)$ at Step 4 is such that
$\tilde{g}(x, a) \equiv g(x, a) \bmod f$. The correctness follows, since $\tilde{g}(x, a) \bmod f$ is the polynomial returned by `BIVARIATEMODULARCOMPOSITION(f, a, \tilde{g})` (see Proposition 3.4).

As required by Algorithm 1 of [75], the tuple of integers $(d, \dots, d, \deg_y(g)) \in \mathbb{Z}^{m+1}$ bounds the column degrees of $(R - v_g)$. Then, since the sum of this tuple is $md + \deg_y(g)$, with $\deg_y(g) = O(n)$, Step 2 costs $\tilde{O}(m^\omega(d + n/m))$ operations [75, Thm. 4.1]. The minimality of $\deg(r)$ implies $\deg(r) \le \deg \det(R) \le md$, and then $v$ has degree at most $\deg \det(R)R^{-1}v_g \le (m-1)d + n$ since $\det(R)R^{-1}$ is the transpose of the cofactor matrix of $R$. Thus the computation of $\tilde{v} = v \bmod r$ in Step 3 uses $\tilde{O}(m(mm+d))$ operations, which is smaller than the cost of Step 2. Next, the matrix-vector product $R\tilde{v}$ can be performed in $\tilde{O}(m^\omega d)$ operations: write the column $\tilde{v}$ of degree $< md$ as $m$ columns of degree $< d$ via $y^d$-adic expansion; use a matrix-matrix product to left-multiply these columns by $R$; finally recombine the resulting columns into a single column which gives $R\tilde{v}$. To obtain $v_g$ it remains to divide each entry of $R\tilde{v}$ by $r$, which costs $\tilde{O}(m^2d)$ since $\deg(R\tilde{v}) < (m+1)d$. By Proposition 3.4, the call at Step 5 uses $c(n, m, d)$ operations. The cost bound in the Proposition follows. $\square$

*Note.* Comparing Proposition 4.4 with Proposition 3.4, note that when $m \sim n^\eta$ and $d \sim n^{1-\eta}$ with $\eta$ from Eq. (3), then the complexity bound of Proposition 4.4 is the same as the one given by the Nüsken-Ziegler algorithm, however the $y$-degree of $g$ can now go up to the order of $n$.

## 4.3 Annihilating polynomials using matrices of relations

Our main algorithm requires an annihilating polynomial for $a$, that is, a polynomial $h$ in $\mathbb{K}[y]$ such that $h(a) \equiv 0 \bmod f$. It can readily be obtained from a matrix of relations.

**PROPOSITION 4.5.** Let $R \in \mathbb{K}[y]^{m \times m}_{\le d}$ be a matrix of relations of $M_m^{(a,f)}$. Its determinant is a nonzero annihilating polynomial for a modulo $f$. It has degree at most $md$ in $\mathbb{K}[y]$ and can be computed from $R$ using $\tilde{O}(m^\omega d)$ operations in $\mathbb{K}$.

**PROOF.** As a polynomial combination of relations in $M_m$, the entry (1, 1) of the (upper triangular) Hermite normal form of a matrix of relations is a relation in $\langle f, y - a \rangle \cap \mathbb{K}[y]$, so it is a nonzero multiple of the minimal polynomial of $a$. This implies the same property for the determinant, since it is a multiple of that entry. The bound on the degree of the determinant is straightforward, and the cost bound is from [49, Thm. 1.1]. $\square$

*Note.* For the computations of the minimal polynomial and of the characteristic polynomial of a modulo $f$, see Section 10.1.

# 5 COMPUTING MATRICES OF RELATIONS

In this section, we give an algorithm computing a matrix of relations. This study may be viewed as a specialization of the formalism developed by Kaltofen and Villard for the block Wiedemann approach (see Sections 1.1.3 and 1.1.4) in terms of manipulations of bivariate polynomials in the ideal generated by $y - a$ and $f$.

As already done in Section 4, notation such as $M_m^{(a,f)}$ and $v_m^{(a,f)}$ is shortened into $M_m$ and $v_m$ in this section, except in the main statements.

In Section 5.1, we show that for $m \in \{1, \dots, n\}$, denominators of irreducible right matrix fraction descriptions of $(yI_n - M_a)^{-1}X$ with $X = (I_m \ 0)^\top \in \mathbb{K}^{n \times m}$ yield bases of $M_m^{(a,f)}$. For efficiency reasons, a further truncation is required: this leads us to introduce modules $M_{\ell,m}^{(a,f)}$ whose bases are the denominators of irreducible right matrix fraction descriptions of $Y^T(yI_n - M_a)^{-1}X$, where $Y^T = (I_\ell \ 0) \in \mathbb{K}^{\ell \times n}$, with $\ell \in \{1, \dots, n\}$; thus we use structured left and right block projections. If
$\ell = n$, $Y$ is the identity matrix of size $n$, and we recover $M_m^{(a,f)}$, but this value is too large for our cost objectives. Instead, we focus on $\ell = m$, and thus $Y = X$.

Section 5.2 describes how a basis of $M_{\ell,m}^{(a,f)}$ can be reconstructed using so-called minimal approximant bases [4, 69], from sufficiently many terms of the power series expansion of the matrix $H = X^T(yI_n - M_a)^{-1}X$.

This strategy is turned into an algorithm for computing matrices of relations in Section 5.3: the expansion of $H$ is obtained via Algorithm `BLOCKTRUNCATEDPOWERS`, while approximant bases are computed using a matrix Padé version of the Berlekamp-Massey algorithm [4, 27]. The correctness and efficiency of this approach depends on a fundamental condition on $M_a$, i.e., on $f$ and $a$ (Proposition 5.6, first item). First, it expresses that the left projection does not prevent us from getting the right denominators of $(yI_n - M_a)^{-1}X$ from those of $H$. It also ensures the existence of matrices of relations of “small” degree, and in this way appropriately limits the number of terms of the expansion of $H$ that are required for the reconstruction. We prove in Section 7 that these properties are satisfied for generic inputs; in Section 8, we further study cases where randomization can ensure such a condition.

Verifying the condition on $M_a$, or verifying that a certain matrix is a matrix of relations, are expensive tasks: except for some restricted cases, the algorithm of Section 5.3 does not certify that its output is indeed a matrix of relations. As such, this would lead to a Monte Carlo composition algorithm. To achieve Las Vegas composition instead, in Section 5.4 we propose an algorithm which either detects that the above-mentioned output is not a matrix of relations, or uses this output to build a certified matrix of relations of slightly larger dimensions.

## 5.1 Matrices of relations as denominators of matrix fractions

This section relates denominators of some matrix fractions to bases of the module of relations $M_m$ and of a truncated version $M_{\ell,m}$ of it.

### 5.1.1 Definitions.

*Matrix Fractions.* We first recall several notions on matrix fractions that can be found in Kailath's book [39, Chap. 6]. Let $N$ be in $\mathbb{K}[y]^{\ell \times m}$, and let $D \in \mathbb{K}[y]^{m \times m}$ be nonsingular. The right fraction $ND^{-1}$ is said to be *irreducible* if $N$ and $D$ are right coprime, i.e. any right divisor common to $N$ and $D$ is unimodular, or equivalently $UN + VD = I_m$ for some $U \in \mathbb{K}[y]^{m \times \ell}$ and $V \in \mathbb{K}[y]^{m \times m}$ [39, Lem. 6.3.5 p. 379]. It is said to be *strictly proper* if for each nonzero entry of $F$, the degree of the numerator is less than the degree of the denominator. It is called a *right fraction description* of $F = ND^{-1} \in \mathbb{K}(y)^{\ell \times m}$. Similarly, $F = \hat{D}^{-1}\hat{N}$ is called a *left fraction description* of $F$. A matrix $F \in \mathbb{K}(y)^{\ell \times m}$ is said to be *describable in degree $d$* if it admits both a left and a right fraction description with denominators of degree at most $d$.

*Truncated Module of Relations.* For efficiency reasons, we consider a $\mathbb{K}[y]$-module similar to $M_m$, but where only the first $\ell$ coefficients of the polynomials are required to be 0, for some positive integer $\ell$, where $n = \deg(f)$. Explicitly, for $\ell, m \in \mathbb{N}_{>0}$ we define the $\mathbb{K}[y]$-modules

$$ M_{\ell,m}^{(a,f)} = \{r(x, y) \in \mathbb{K}[x, y]_{<(m, \ell)} \mid [a(x)^k r(x, a(x)) \text{ rem } f]_0^{\ell-1} = 0 \text{ for all } k \ge 0\}, $$

together with the usual simplified notation $M_{\ell,m}$. They satisfy the inclusions $M_{1,m} \supseteq M_{2,m} \supseteq \dots \supseteq M_{n,m} = M_m$. The determinantal degree of $M_{\ell,m}$ is denoted $v_{\ell,m}$. Of particular interest is the case when $M_{m,m} = M_m$.
### 5.1.2 Relation between bases of relations and denominators of matrix fractions.

**PROPOSITION 5.1.** For $\ell, m \in \{1, \dots, n\}$, the columns of a matrix $D \in \mathbb{K}[y]^{m \times m}$ form a basis of $\mathcal{M}_{\ell, m}^{(a,f)}$ if and only if $D$ is the denominator of an irreducible right fraction description $ND^{-1}$ of

$$(I_{\ell} \ 0)(yI_{n} - M_{a})^{-1}X \in \mathbb{K}[y]^{\ell \times m};$$

the denominator of any right fraction description of this matrix is a right multiple of any such basis $D$.

### 5.1.3 Proof of Proposition 5.1. For a matrix of rational functions $F \in \mathbb{K}(y)^{\ell \times m}$, we let

$$\mathcal{D}(F) = \{v \in \mathbb{K}[y]^m \mid Fv \in \mathbb{K}[y]^\ell\}, \quad (15)$$

which is a $\mathbb{K}[y]$-submodule of $\mathbb{K}[y]^m$ of rank $m$. Then, we can establish the relation between the module $\mathcal{M}_{\ell, m}$ and the matrix in Proposition 5.1.

**LEMMA 5.2.** For $\ell, m$ in $\{1, \dots, n\}$, one has $\mathcal{M}_{\ell, m} = \mathcal{D}((I_\ell \ 0)(yI_n - M_a)^{-1}X)$.

**PROOF.** Taking $Y^\top = (I_\ell \ 0)$, define $H(y) = Y^\top(yI_n - M_a)^{-1}X$ and $H_k = Y^\top M_a^k X \in \mathbb{K}^{\ell \times m}$, so that, by power series expansion in $y^{-1}$,

$$H(y) = \sum_{k \ge 0} H_k y^{-k-1} = \sum_{k \ge 0} Y^\top M_a^k X y^{-k-1}.$$

Let $r(x, y) = \sum_{0 \le i \le d} r_i(x) y^i \in \mathbb{K}[x, y]_{<(m, \cdot)}$ be of $y$-degree $d$, and let $v_i \in \mathbb{K}^m$ be the coefficient vector of $r_i$ for $i = 0, \dots, d$. Then, for $k \ge 0$,

$$[a^k r(x, a) \operatorname{rem} f]_0^{\ell-1} = \left[ \sum_{0 \le i \le d} a^{k+i} r_i \operatorname{rem} f \right]_0^{\ell-1} = \sum_{0 \le i \le d} [a^{k+i} r_i \operatorname{rem} f]_0^{\ell-1}$$

and $[a^{k+i} r_i \operatorname{rem} f]_0^{\ell-1}$ has coefficient vector $Y^\top M_a^{k+i} X v_i = H_{k+i} v_i$. Hence, $[a^k r(x, a) \operatorname{rem} f]_0^{\ell-1}$ has coefficient vector $H_k v_0 + \dots + H_{k+d} v_d$. Therefore $r(x, y)$ is in $\mathcal{M}_{\ell, m}$ if and only if

$$H_k v_0 + \dots + H_{k+d} v_d = 0 \quad \text{for all } k \ge 0. \quad (16)$$

On the other hand, defining $H_k = 0$ for $k < 0$, the expansion of $Hv$ at infinity reads

$$Hv = \sum_{k \ge 0} H_k y^{-k-1} \sum_{0 \le i \le d} v_i y^i = \sum_{k \ge -d} (H_k v_0 + \dots + H_{k+d} v_d) y^{-k-1}, \quad (17)$$

which implies that Eq. (16) holds if and only if $Hv$ has polynomial entries. $\square$

**Proposition 5.1** is then a direct consequence of the following general result on matrix fractions, which is a reformulation of [39, Thm. 6.5-4 and Lem. 6.5-5, p. 441].

**LEMMA 5.3.** Let $F \in \mathbb{K}(y)^{\ell \times m}$ be a matrix of rational fractions. The columns of $D \in \mathbb{K}[y]^{m \times m}$ form a basis of $\mathcal{D}(F)$ if and only if $D$ is the denominator of an irreducible right fraction description $ND^{-1}$ of $F$. Besides, the denominator of any right fraction description of $F$ is a right multiple of such a $D$.

### 5.1.4 Notes.
The role of the truncated modules $\mathcal{M}_{\ell, m}$ is to reduce the cost of computations: we decrease the dimension of the relevant matrices using a structured left projection. The more usual approach [46] uses generic projections matrices; our choice here is similar to the one used for the efficient computation of generic resultants [71].

Although not used in this work, genericity on the left is sufficient: if $V \in \mathbb{K}^{n \times \ell}$ is generic with $\ell \in \{m, \dots, n\}$, then one has $\mathcal{M}_m = \mathcal{D}((yI_n - M_a)^{-1}X) = \mathcal{D}(Y^\top(yI_n - M_a)^{-1}X)$. The latter occurs if and only if $\operatorname{rank}(V^\top P, V^\top PA, V^\top PA^2, \dots) = v_m$ for a well chosen full rank matrix $P \in \mathbb{K}^{n \times v_m}$, and a restriction $A \in \mathbb{K}^{v_m \times v_m}$ of $M_a$ to the invariant subspace generated by $X$ [70, Lem. 4.2]. The rank condition is satisfied for a generic projection [70, Cor. 6.4 and its proof].
In terms of generators of matrix sequences, Eq. (16) shows that the denominators of Proposition 5.1 are bases of modules of vector generators for the matrix sequence $\{(I_\ell 0)M_a^k X\}_{k \ge 0}$ [46, Lem. 2.8].

## 5.2 Reconstructing denominators of matrix fractions via approximant bases

Algorithm `BLOCKTRUNCATEDPOWERS` from Section 3.3 allows one to compute a truncated power series expansion of $H(y) = X^T(yI_n - M_a)^{-1}X$. When the precision of this expansion is sufficient, a basis of $M_{m,m}$ can be reconstructed.

### 5.2.1 Definitions.

**Weak Popov matrices.** Let $P \in \mathbb{K}[y]^{m \times m}$ be a matrix whose column $j$ has degree $d_j \ge 0$. The (column) leading matrix of $P$ is the matrix in $\mathbb{K}^{m \times m}$ whose entry $(i, j)$ is the coefficient of degree $d_j$ of the entry $(i, j)$ of $P$. Then $P$ is said to be (column) reduced if its leading matrix is invertible. This is the case if and only if [39, Eq. (24), p. 384]

$$ \deg \det(P) = d_1 + \cdots + d_m. \qquad (18) $$

A (column) reduced matrix is in (column) weak Popov form if its leading matrix is invertible and upper triangular. Any submodule of $\mathbb{K}[y]^m$ has at least one basis which is in weak Popov form [5, 39].

**Approximant bases.** Let $F \in \mathbb{K}[[y]]^{m \times k}$ be a matrix of power series and $\sigma \in \mathbb{N}$ be a nonnegative integer. A matrix $P \in \mathbb{K}[y]^{k \times k}$ is an *approximant basis* of $F$ at order $\sigma$ if its columns form a basis of the $\mathbb{K}[y]$-module $\{v \in \mathbb{K}[y]^k \mid Fv = 0 \text{ mod } y^\sigma\}$, which is free of rank $k$. This approximant basis is said to be *minimal* if it is reduced. Minimal approximant bases are also called $\sigma$-bases, or order bases [4, 69].

### 5.2.2 Denominators from approximant bases.
We are going to use approximant bases for solving equations of the type of Eq. (6). As pointed out in Section 1.2, we use expansions at $y = 0$ rather than infinity (see Remark 5.7).

**PROPOSITION 5.4.** Let $H \in \mathbb{K}(y)^{m \times m}$ be strictly proper, and $\delta$ be the determinantal degree of $\mathcal{D}(H)$ (notation from Eq. (15)). Suppose that $H$ has a power series expansion $H = \sum_{k \ge 0} S_k y^k$ at $y = 0$, with $S_k \in \mathbb{K}^{m \times m}$. Let

$$ F = \left( \sum_{k=0}^{2d-1} S_k y^k - I_m \right) \in \mathbb{K}[y]^{m \times (2m)}, $$

and let

$$ P = \begin{pmatrix} D & P_1 \\ N & P_2 \end{pmatrix} \in \mathbb{K}[y]^{(2m) \times (2m)} $$

be an approximant basis at order $2d$ of $F$ in weak Popov form, with each submatrix of size $m \times m$. Then the following properties hold:

(i) $D$ is weak Popov; $\deg(N) < \deg(D)$; the sum of the degrees of the diagonal entries of $D$ is $\deg \det(D)$ and satisfies $\deg \det(D) \le \delta$.

(ii) If $\deg \det(D) = \delta$ and each of the $m$ rightmost columns of $P$ has degree at least $\deg(D)$, then $ND^{-1}$ is an irreducible description of $H$.

(iii) If $H$ is describable in degree $d$, then $ND^{-1}$ is an irreducible description of $H$ such that $\deg(D) \le d$ and each of the $m$ rightmost columns of $P$ has degree at least $\deg(D)$.

The first item gives general properties of the approximant basis in weak Popov form, whereas Items (ii) and (iii) give sufficient conditions to guarantee it recovers an irreducible fraction description of $H$.
### 5.2.3 Proof of Proposition 5.4.

**LEMMA 5.5.** Let $P \in \mathbb{K}[y]^{m \times m}$. If $P$ is reduced and $B \in \mathbb{K}[y]^{m \times m}$ is a right-multiple $B = PU$ with $U$ nonsingular, then $\deg(P) \le \deg(B)$. If $P \in \mathbb{K}[y]^{m \times m}$ is weak Popov, with diagonal entries of respective degrees $d_1, \dots, d_m \in \mathbb{N}$, and $v \in \mathbb{K}[y]^m$ is a nonzero right-multiple of $P$ whose bottommost entry of largest degree is in row $i$ and has degree $d$, then $d_i \le d$.

**PROOF.** The first claim follows from the predictable degree property [39, Thm. 6.3-13, p. 387].
The second one is from [56, Lem. 1.17]. $\square$

We now prove Proposition 5.4. Consider an irreducible fraction description $QR^{-1} = H$ for some $Q \in \mathbb{K}[y]^{m \times m}$ and some weak Popov $R \in \mathbb{K}[y]^{m \times m}$. Since $H$ is strictly proper we have $\deg(Q) < \deg(R)$, and thus the $i$th column of $(\begin{smallmatrix} R \\ Q \end{smallmatrix})$ has its bottommost entry of largest degree in row $i$; let $d_i$ be this degree.

In Item (i), the first two claims follow from the definition of $P$ being weak Popov. In particular $D$ is column reduced, hence Eq. (18) shows that $\deg \det(D)$ is the sum of column degrees of $D$, which is also the sum of diagonal degrees of $D$ since $D$ is weak Popov. The identity $F(\begin{smallmatrix} R \\ Q \end{smallmatrix}) = 0$ implies the same identity modulo $y^{2d}$, and therefore $(\begin{smallmatrix} R \\ Q \end{smallmatrix})$ is a right-multiple of $P$. Hence, by Lemma 5.5, $d_i$ is at least the degree of the $i$th column of $P$, which is the degree of the $i$th column of $D$; it follows that $\deg \det(D) \le d_1 + \cdots + d_m$. On the other hand, since $R$ is reduced we have $d_1 + \cdots + d_m = \deg \det(R)$, proving the last claim of Item (i).

Concerning Item (ii), the assumption $\deg \det(D) = \delta = \deg \det(R)$ implies that the sum of column degrees of $D$ is $d_1 + \cdots + d_m$, while as showed above the $i$th column of $D$ has degree at most $d_i$. Thus $D$ has the same column degrees $(d_1, \dots, d_m)$ as $R$. In particular $\deg(D) = \deg(R) > \deg(Q)$. Then, since by assumption the $m$ rightmost columns of $P$ have bottommost entries of largest degree in rows at least $m+1$ and of degree at least $\deg(D)$, it follows from Lemma 5.5 that $(\begin{smallmatrix} R \\ Q \end{smallmatrix})$ is a right-multiple of the leftmost $m$ columns of $P$. This means $(\begin{smallmatrix} R \\ Q \end{smallmatrix}) = (\begin{smallmatrix} D \\ N \end{smallmatrix})U = (\begin{smallmatrix} DU \\ NU \end{smallmatrix})$ for some $U \in \mathbb{K}[y]^{m \times m}$, and $U$ is unimodular since $R$ and $D$ are nonsingular with $\deg \det(R) = \deg \det(D)$. Hence $H = QR^{-1} = ND^{-1}$ and the fraction $ND^{-1}$ is irreducible.

The following proof of Item (iii) reflects that of [27, Lem. 3.7]. The assumption implies first the existence of a left fraction $H = \hat{R}^{-1}\hat{Q}$ with $\deg(\hat{Q}) < \deg(\hat{R}) \le d$, and second the degree bound $\deg(R) \le d$ thanks to the degree minimality of reduced bases (see Lemma 5.5). The above paragraph shows in particular $\deg(D) \le \max_i(d_i) = \deg(R) \le d$.

Now, since $\hat{R}(\sum_{0 \le k < 2d} S_k y^k) = \hat{Q} \mod y^{2d}$, left-multiplying by $\hat{R}$ both sides of $F(\begin{smallmatrix} D \\ N \end{smallmatrix}) \equiv 0 \mod y^{2d}$ shows that $\hat{Q}D - \hat{R}N$ is a right-multiple of $y^{2d}\hat{R}$. On the other hand, $\hat{Q}D - \hat{R}N$ has degree less than $2d$. Hence it is zero, and $H = \hat{R}^{-1}\hat{Q} = ND^{-1}$. To prove that the latter fraction is irreducible, assume by contradiction that $D$ and $N$ have a nonsingular common right divisor $B \in \mathbb{K}[y]^{m \times m}$, with $\deg \det(B) > 0$. Then $H = (NB^{-1})(DB^{-1})^{-1}$ yields $F(\begin{smallmatrix} DB^{-1} \\ NB^{-1} \end{smallmatrix}) \equiv 0 \mod y^{2d}$, and $P \operatorname{diag}(B^{-1}, I_m)$ is a right-multiple of $P$ (since $P$ is a basis): this is impossible since $\deg \det(P \operatorname{diag}(B^{-1}, I_m)) < \deg \det(P)$.

It remains to prove the last degree assertion. By contradiction, assume that $P$ has a column $(\begin{smallmatrix} v_0 \\ v_1 \end{smallmatrix})$ of index larger than $m$ with $v_0$ and $v_1$ in $\mathbb{K}[y]^m$ both of degree less than $d$. Then an argument similar to the one above shows that $\hat{Q}v_0 - \hat{R}v_1 = 0$. Altogether we obtain a matrix $(\begin{smallmatrix} D \\ N \end{smallmatrix} v_0)$ of rank $m+1$ which is in the right kernel of $(\hat{Q} \ \hat{R}) \in \mathbb{K}[y]^{m \times (2m)}$ whose rank is $m$: this is not possible.

### 5.2.4 Notes.
The existence of appropriate left and right descriptions of $H$ was used before for the reconstruction of matrix fractions within the approximant framework [27, Sec. 3.2]. Our proof is similar to that of [27, Lem. 3.7], with the additional use of the weak Popov form.

Reduced forms were introduced [73] as a way to get a better control over the degrees when computing with polynomial matrices and matrix fractions, see e.g. [39, Lem. 6.3-11, p. 385] for
**Algorithm 5.1 CANDIDATEBASIS(f, a, m, d)**

**Input:** $f \in \mathbb{K}[x]$ of degree $n$, with $f(0) \neq 0$, $a \in \mathbb{K}[x]_{<n}$ with $\gcd(a, f) = 1$, $m \le n$ and $d$ in $\mathbb{N}_{>0}$

**Output:** a weak Popov matrix $R \in \mathbb{K}[y]_{\le 2d}^{m \times m}$ and a flag in $\{\text{CERT}, \text{NoCERT}\}$; $R$ is a basis of $M_m$ in either of the following cases:

* $v_{m,m} = v_m$ and $H = X^T(yI_n - M_a)^{-1}X$ is describable in degree $d$, in which case $\deg(R) \le d$

* the flag is CERT, which implies $v_{m,m} = v_m = n$

1. Truncated expansion of $H$: compute $S_k = -X^T M_a^{-k-1} X$ for $k < 2d$ using Algorithm 3.8

$$ (A_{i,k}^*)_{\substack{0 \le i < m \\ 0 \le k < 2d+2}} \leftarrow \text{BLOCKTRUNCATEDPOWERS}(f, a^{-1} \mod f, m, 2(d+1)) $$

$$ S_{i,k} \in \mathbb{K}^m \leftarrow \text{vector of coefficients of } -A_{i,k+1}^* \in \mathbb{K}[x]_{<m}, \text{ for } 0 \le i < m \text{ and } 0 \le k < 2d $$

2. Fraction reconstruction: compute approximant basis using algorithm from [27, 38]

$$ F \in \mathbb{K}[y]_{<2d}^{m \times 2m} \leftarrow (\sum_{0 \le k < 2d} S_k y^k - I_m) \text{ where } S_k = (S_{0,k} \cdots S_{m-1,k}) \in \mathbb{K}^{m \times m} $$

$$ P \in \mathbb{K}[y]_{<2d}^{2m \times 2m} \leftarrow \text{PM-BASIS}(F^T, 2d, 0)^T, \text{ with } P \text{ in weak Popov form} $$

3. Return candidate matrix and result of basic certification

$$ R \leftarrow P_{1..m,1..m} $$

if the sum of diagonal degrees of $R$ is equal to $n$

and each of the $m$ rightmost columns of $P$ has degree $\ge \deg(R)$

then return $(R, \text{CERT})$ else return $(R, \text{NoCERT})$

► Item (ii) of Proposition 5.4

proper fractions, [39, Thm. 6.3-13, p. 387] for a predictable degree property, and [39, Thm. 6.5-10, p. 458] concerning the *minimality* of the column degrees. Weak Popov forms were introduced later [5, 55] (under the name quasi-Popov and up to column permutation) and provide a refined degree control as illustrated by Lemma 5.5.

## 5.3 Candidate basis of relations

Algorithm CANDIDATEBASIS takes as input a polynomial $f \in \mathbb{K}[x]$ of degree $n$ with $f(0) \neq 0$, a polynomial $a \in \mathbb{K}[x]_{<n}$ such that $\gcd(a, f) = 1$, and two positive integers $m \le n$ and $d$. With this input, it computes an $m \times m$ matrix of degree at most $2d$.

The algorithm starts by computing a truncated expansion at order $2d$ of $H = X^T(yI_n - M_a)^{-1}X$ at $y=0$ using Algorithm BLOCKTRUNCATEDPOWERS. Then, it computes a $2m \times 2m$ minimal approximant basis as in Proposition 5.4 using the algorithm PM-BASIS of [27], and extracts a potential basis of relations. In some cases we can certify that it is a indeed basis of $M_m$, but it is not always possible to do so; a flag is returned to indicate this. This certification is actually an optimization, rather than strictly necessary; Section 5.4 discusses this question in more detail.

**PROPOSITION 5.6.** Given $f \in \mathbb{K}[x]$ of degree $n$ with $f(0) \neq 0$, $a \in \mathbb{K}[x]_{<n}$ such that $\gcd(a, f) = 1$, and two positive integers $m \le n$ and $d$, Algorithm CANDIDATEBASIS uses $\tilde{O}(m^\omega d + c(n, m, d))$ operations in $\mathbb{K}$, with $c(\cdot)$ from Eq. (2), and computes a weak Popov matrix $R \in \mathbb{K}[y]_{\le 2d}^{m \times m}$. The matrix $R$ is a basis of $M_m^{(a,f)}$ in either of the following cases:

* The determinantal degree $v_{m,m}^{(a,f)}$ is equal to $v_m^{(a,f)}$ and the fraction $H(y) = X^T(yI_n - M_a)^{-1}X$ is describable in degree $d$; in that case we further have $\deg(R) \le d$; if in addition $v_m^{(a,f)} = n$ then the flag is CERT.

* The flag is CERT, which implies $v_{m,m}^{(a,f)} = v_m^{(a,f)} = n$.

**PROOF.** Proposition 3.7 shows that Step 1 uses $\tilde{O}(m^2d + c(n, m, d))$ operations to compute the vectors $S_{i,k} \in \mathbb{K}^m$. These vectors are such that the matrices $S_k$ built in Step 2 are $S_k = -X^T M_a^{-k-1} X$; as a result, the matrix $S = \sum_{0 \le k < 2d} S_k y^k$ considered at Step 2 is the power series expansion of $H$
truncated at order 2d. Then Step 2 correctly computes a weak Popov approximant basis P for $F = (S - I_m)$ at order 2d with $\deg(P) \le 2d$ using $\tilde{O}(m^\omega d)$ operations [27, Thm. 2.4] [38, Prop. 3.2]. (Note that transposes are used at Step 2 because in [27, 38] approximant bases are considered row-wise, rather than column-wise here.) The claimed cost bound for Algorithm CANDIDATEBASIS is proved.

For the first item, assume that $H$ is describable in degree $d$. Then Item (iii) of Proposition 5.4 ensures that $R$ is the denominator of an irreducible right fraction description of $H$, that $\deg(R) \le d$, and that each of the $m$ rightmost columns of $P$ has degree at least $\deg(R)$. From Proposition 5.1 we obtain that $R$ is a basis of $M_{m,m}$, hence a basis of $M_m$ when $v_{m,m} = v_m$. This also proves the last claim of the item: if $v_{m,m} = v_m = n$, then $\deg \det(R) = n$ and this is the sum of diagonal degrees of $R$ since this matrix is in weak Popov form; hence the flag CERT is returned.

For the second item, assume that the output flag is CERT. Then the sum of diagonal degrees of $R$ is $n$; according to Item (i) of Proposition 5.4, this sum is also $\deg \det(R)$ and is at most $\delta$, the determinantal degree of bases of $\mathcal{D}(H)$. On the other hand Proposition 5.1 implies that $\delta$ is the determinantal degree $v_{m,m}$ of $M_{m,m}$. Hence $n = \deg \det(R) \le \delta = v_{m,m}$, from which we deduce $\deg \det(R) = \delta = v_{m,m} = v_m = n$, since $v_{m,m} \le v_m \le n$ always holds. Since the output flag is CERT we know in addition that each of the $m$ rightmost columns of $P$ has degree at least $\deg(R)$. Thus Item (ii) of Proposition 5.4 applies, and $R$ is the denominator of an irreducible right fraction description of $H$. We conclude as done for the first item that $R$ is a basis of $M_m$. $\square$

**REMARK 5.7.** The assumption that $f$ and $a$ are coprime is used here to ensure that $M_a$ is invertible, so that the expansion $H = \sum_{k \ge 0} S_k y^k = \sum_{k \ge 0} (-X^T M_a^{-k-1} X) y^k$ at $y=0$ can be used for fraction reconstruction. This is different from what happened in the proof of Proposition 5.1, where we used the expansion at infinity $H = \sum_{k \ge 0} H_k y^{-k-1}$.

This assumption on $\gcd(f, a)$ is harmless in our context: in the computation of $g(a)$ the $f$, one can instead evaluate $g(y-c)$ at $y=a+c$ for a randomly chosen $c \in \mathbb{K}$, ensuring $\gcd(a+c, f) = 1$ with good probability. See Step 2 in Algorithm MODULARCOMPOSITIONBASECASE.

*Notes.* For some families of approximation instances, PM-BASIS has been used to design faster minimal approximant basis algorithms [38, 74]. Yet, the instances considered here are ones where PM-BASIS is the fastest known algorithm.

A candidate matrix of relations in $K[y]_{\le 2d}^{m \times m}$ corresponds to $m$ polynomials in $K[x, y]_{\le (m, 2d)}$. Using Algorithm SIMULTANEOUSBIVARIATEMODULARCOMPOSITION to verify that the evaluations of these polynomials at $a \bmod f$ are zero uses $\tilde{O}(c(n, m, d^2))$ operations in $K$, by Lemma 3.3. For the values of $m$ and $d$ used to obtain the exponent $\kappa < 1.43$ in our main algorithm, this is $O(n^{2.55})$, and thus too costly.

## 5.4 Certified matrix of relations

In general, when Algorithm CANDIDATEBASIS does not certify its result, we do not know methods to verify that the matrix it returns is a matrix of relations within our complexity bound.

Instead, from a matrix $R$ computed by Algorithm CANDIDATEBASIS, Algorithm MATRIXOFRELATIONS either detects that it is not a matrix of relations of $M_m$, or constructs from $R$ a matrix of relations of $M_{m'}$ of degree at most $2d$, for some $m' < 2m$. This is the key towards making our modular composition algorithm Las Vegas, rather than Monte Carlo.

To achieve this, instead of evaluating all columns of $R$ at $a \bmod f$, Algorithm MATRIXOFRELATIONS evaluates only two polynomials built randomly from these columns (and only one polynomial in the special case $m=1$), which is within our target complexity using the Nüsken-Ziegler algorithm. If these evaluations are not both zero, then $R$ was not a matrix of relations. Otherwise the algorithm constructs a Sylvester matrix from these two vectors (see e.g., [22, Sec. 6.3] for the
definition and properties of the Sylvester matrix). When this matrix is nonsingular, it is a matrix of relations of a module $M_{m'}$ for $m' \le \max(1, 2(m-1))$; since $m'$ cannot be much larger than $m$, this matrix can be used for efficient composition.

<caption>**Algorithm 5.2 MATRIXOFRELATIONS(f, a, m, d, (rᵢ)₃≤i≤m)**</caption>

**Input:** $f \in \mathbb{K}[x]$ of degree $n$, with $f(0) \ne 0$, $a \in \mathbb{K}[x]_<n$ with $\gcd(a, f) = 1$, $m \le n$ and $d$ in $\mathbb{N}_{>0}$, $(r_i)_{3\le i \le m} \in \mathbb{K}^{m-2}$
**Output:** either **FAIL** or a matrix $R' \in \mathbb{K}[y]^{m' \times m'}_{\le 2d}$ of relations of $M_{m'}$ with $m' \le \max(1, 2(m-1))$
1: ▶ **Use Algorithm 5.1 to find a candidate basis of relations**
   $(R, \text{FLAG}) \in \mathbb{K}[y]^{m \times m}_{\le 2d} \times \{\text{CERT}, \text{NoCERT}\} \leftarrow \text{CANDIDATEBASIS}(f, a, m, d)$
   **if** $\text{FLAG} = \text{CERT}$ **then return** $R$
2: ▶ **Case m = 1, check that R_{1,1} \in K[y]_{<2d+1} annihilates a mod f**
   **if** $m = 1$ **then**
      **if** **MODULARCOMPOSITION-BRENTKUNG(f, a, R_{1,1}) \neq 0** **then return** **FAIL**
      **else return** $R$
3: ▶ **Build candidate relations and verify them**
   $r(x, y) \leftarrow R_{*,1}$; $s(x, y) \leftarrow R_{*,2} + r_3 R_{*,3} + \dots + r_m R_{*,m}$
   ▶ both in $\mathbb{K}[x, y]_{< (m, 2d+1)}$
   **if** **BIVARIATEMODULARCOMPOSITION(f, a, r) \neq 0**
   **or** **BIVARIATEMODULARCOMPOSITION(f, a, s) \neq 0**
   **then return** **FAIL**
   ▶ Algorithm 3.5
   **if** $m = 2$ **then return** $R$
4: ▶ **Construct and return the Sylvester matrix of f and s, if it is nonsingular**
   **if** $\gcd_x(r, s) \ne 1$ **then return** **FAIL**
   ▶ $r$ and $s$ not coprime as elements of $\mathbb{K}(y)[x]$ **return** the Sylvester matrix of $(r, s)$ as in [22, Sec. 6.3, Eq. (5)], with rows in reversed order, viewing $r$ and $s$ as polynomials in $x$ over $\mathbb{K}[y]$

**PROPOSITION 5.8.** Given $f \in \mathbb{K}[x]$ of degree $n$ with $f(0) \ne 0$, $a \in \mathbb{K}[x]_<n$ with $\gcd(a, f) = 1$, two positive integers $m \le n$ and $d$, and $(r_i)_{3\le i \le m} \in \mathbb{K}^{m-2}$, Algorithm MATRIXOFRELATIONS uses $\tilde{O}(m^\omega d + c(n, m, d))$ operations in $\mathbb{K}$, with $c(\cdot)$ from Eq. (2), and returns either **FAIL** or a matrix of relations $R' \in \mathbb{K}[y]^{m' \times m'}$ of $M_{m'}^{(a,f)}$ where $m' \le \max(1, 2(m-1))$.

If $v_{m,m}^{(a,f)} = v_m^{(a,f)}$, the fraction $H = X^T(yI_n - M_a)^{-1}X$ is describable in degree $d$, and $(r_3, \dots, r_m)$ are chosen uniformly and independently at random from a finite subset $S$ of $\mathbb{K}$, then failure happens with probability at most $(m-1)/\text{card}(S)$ and in case of success, $\deg(R') \le d$.

**PROOF.** If $\text{FLAG} = \text{CERT}$ at Step 1, then from the second item of Proposition 5.6 an appropriate matrix of relations is returned. Now assume that $\text{FLAG} = \text{NoCERT}$ and Algorithm MATRIXOFRELATIONS does not return **FAIL**. If $m=1$ then the relation has been checked at Step 2, proving the result. Otherwise, let $R' \in \mathbb{K}[y]^{m' \times m'}$ be the output matrix, which is constructed from the polynomials $r$, $s$ of $x$-degree less than $m$; in particular, $m' = \deg_x(r) + \deg_x(s) \le 2(m-1)$ [22, Sec. 6.3]. The fact that the test at Step 4 has not failed ensures that $r$ and $s$ are coprime as univariate polynomials in $\mathbb{K}(y)[x]$, and therefore $R'$ is nonsingular [22, Cor. 6.15]. Furthermore, since the tests at Step 3 have not failed, $r$ and $s$ are relations of $M_m$. It follows that the columns of $R'$, which are by construction multiples of $r$ and $s$ in $\mathbb{K}[x, y]$ represented as vectors in $\mathbb{K}[y]^{m'}$, are relations of $M_{m'}$. Besides, the construction of the Sylvester matrix does not increase the $y$-degree, hence $\deg(R') \le \deg(R) \le 2d$. We have proved the fact that if the output is not **FAIL**, then it is a matrix of relations of $M_{m'}$.

For the complexity bound, the cost for finding $R$ is given in Proposition 5.6, while the ones for checking that $R_{1,1}$, $r$ and $s$ are relations are given in Propositions 3.1 and 3.4. As for the gcd test
at Step 4, it can be done via the resultant of $r$ and $s$ with respect to $x$, computed using $\tilde{O}(m^2d)$ operations [63].

It remains to prove the third assertion and the probability bound. Since when $FLAG = \text{CERT}$ a basis is returned with no randomization, assume $FLAG = \text{NoCERT}$. The assumptions here and the first item of Proposition 5.6 ensure that $R$ is a basis of $M_m$ with $\deg(R) \le d$, hence $\deg(R') \le d$. In that case failure never occurs at Step 2 for $m=1$. It never occurs either at Step 3 for $m \ge 2$, and $r$ and $s$ are relations of $M_m$. The columns of $R$ represent bivariate polynomials $b_1, \dots, b_m \in \mathbb{K}[x, y]_{<(m,d+1)}$ and we claim that $\gcd_x(b_1, \dots, b_m) = 1$, meaning that there is a $\mathbb{K}[y]$-linear combination of $b_1, \dots, b_m$ which is in $\mathbb{K}[y] \setminus \{0\}$. Since $R$ is nonsingular, the first column of a transformation for the (upper triangular) Hermite normal form of $R$ provides such a combination. It follows that $\text{FAIL}$ is returned with probability at most $(m-1)/\text{card}(S)$ at Step 4 [22, Thm. 6.46]. $\square$

*Note.* The computation of *CERT* by Algorithm *CANDIDATEBASIS* is only an optimization. Algorithm *MATRIXOFRELATIONS* works as it is, even if *CERT* is never returned. When the candidate matrix *R* at Step 1 is not a matrix of relations, this is often detected at Step 3, but not always. Even if it is not detected, it suffices to find two coprime polynomials *r*(*x*, *y*) and *s*(*x*, *y*) that are relations to ensure that Algorithm *MATRIXOFRELATIONS* returns a matrix of relations. For example, it may happen that *R* is not a matrix of relations but some columns of it still give low-degree relations of $M_m$.

# 6 CHANGE OF BASIS

In this section we present an algorithm for performing a change of basis in $A = \mathbb{K}[x]/\langle f \rangle$. This algorithm is used in a randomized manner in Section 8, in order to handle arbitrary inputs with good complexity bounds. Our approach is based on an extension of the approximant bases used in Section 5; we start with necessary definitions.

## 6.1 Definitions

We use an extension of the forms of polynomial matrices introduced in Section 5.2.1, called *shifted* forms [5, 69]. For a given tuple $t = (t_1, \dots, t_m) \in \mathbb{Z}^m$ and a column vector $v \in \mathbb{K}[y]^m$, the *t*-shifted degree of $v$ is $\max_{1 \le i \le m}(\deg(v_i) + t_i)$. Then, for a matrix $P \in \mathbb{K}[y]^{m \times m}$ whose jth column has *t*-shifted degree $d_j \in \mathbb{Z}$, the (*column*) *t*-shifted leading matrix of $P$ is the matrix in $\mathbb{K}^{m \times m}$ whose entry $(i, j)$ is the coefficient of degree $d_j - t_i$ of the entry $(i, j)$ of $P$. Then $P$ is said to be *t*-shifted weak Popov if this *t*-shifted leading matrix is invertible and upper triangular.

We also need the corresponding normal form: $P$ is said to be *t*-shifted Popov if it is *t*-shifted weak Popov and its *row* leading matrix is the identity of $\mathbb{K}^{m \times m}$ [5, 39]. For a given $t$, any submodule of rank $m$ of $\mathbb{K}[y]^m$ admits a unique basis in *t*-shifted Popov normal form [5, Thm. 3.7]. By definition, *t*-shifted Popov matrices are also (non-shifted) *row* reduced; in particular, Hermite normal forms are *t*-shifted Popov for an appropriate choice of $t$, hence are *row* reduced.

Row reduced matrices allow for a *division with remainder* with stronger properties than the one for general nonsingular matrices presented in Section 4.2.1; namely they ensure uniqueness of the remainder. Precisely, if a matrix $P \in \mathbb{K}[y]^{m \times m}$ is *row* reduced, for any vector $v \in \mathbb{K}[y]^m$ there exists a unique vector $\tilde{v} \in \mathbb{K}[y]^m$ such that $v - \tilde{v}$ is a right multiple of $P$ and the $i$th entry of $\tilde{v}$ has degree less than the $i$th row of $P$ [39, Thm. 6.3-15, p. 389].

We also use the fact that, by definition, for any block decomposition $P = (P_{11} \ P_{12})$ of a matrix $P$, if $P$ is in Hermite (resp. *t*-shifted Popov) normal form, then:

* $P_{11}$ and $P_{22}$ are in Hermite normal form (resp. in shifted Popov normal form with respect to the corresponding subtuple of $t$);

* each column of $P_{12}$ (resp. $P_{21}$) is its own remainder in the division by $P_{11}$ (resp. $P_{22}$).
Finally, t-shifted forms induce the notion of *t-shifted approximant bases* [5, 38, 74], which are approximant bases (see Section 5.2.1) in *t-shifted Popov normal form*.

## 6.2 Inverse modular composition and change of basis via approximant bases

Let $f$ be in $\mathbb{K}[x]$ of degree $n$. A core ingredient for the randomization in our composition algorithm is an instance of *inverse modular composition*, which is used to change the basis of $\mathbb{A} = \mathbb{K}[x]/\langle f \rangle$ from $(1, x, \dots, x^{n-1})$ to $(1, \gamma, \dots, \gamma^{n-1})$ mod $f$, for some $\gamma \in \mathbb{K}[x]$ whose minimal polynomial $\mu_{\gamma}$ modulo $f$ has degree $n$. This change of basis induces the $\mathbb{K}$-algebra isomorphism

$$ \phi_{\gamma}: \mathbb{A} \to \mathbb{K}[y]/\langle\mu_{\gamma}\rangle, \qquad (19) $$

which maps any $u \in \mathbb{A}$ to $v$ such that $v(\gamma) \equiv u \mod f$. Given $a$ in $\mathbb{K}[x]_{<n}$, this section explains how to compute the unique polynomial representative $\alpha \in \mathbb{K}[y]_{<n}$ of $\phi_{\gamma}(a \mod f)$, i.e., the unique $\alpha \in \mathbb{K}[y]_{<n}$ such that $\alpha(\gamma) \text{rem } f = a$.

Reversing the path followed in our modular composition approach, we first find a bivariate $\tilde{\alpha} \in \mathbb{K}[x, y]$ such that $\tilde{\alpha} - \alpha \in \mathcal{M}_m^{(\gamma,f)}$, hence $\tilde{\alpha}(x, \gamma) \equiv \alpha(\gamma) \mod f$. Then the univariate solution $\alpha$ is recovered from $\tilde{\alpha}$ and a basis of relations $R$ of $\mathcal{M}_m^{(\gamma,f)}$ by reversing the division from Eq. (14); this corresponds to a division by the Hermite normal form of $R$.

Our algorithm for computing $\tilde{\alpha}$ can be seen as a generalization to $m \ge 1$ of Shoup's algorithm for computing $\alpha$, mentioned in Section 1.1.2. The latter algorithm deals with the case $m = 1$: from the power projections $(\ell(1), \ell(\gamma), \dots, \ell(\gamma^{2n-1}))$ and $(\ell(a), \ell(\gamma a), \dots, \ell(\gamma^{n-1}a))$, it obtains both $\alpha$ and $\mu_{\gamma}$ by solving two Padé approximation problems. In the matrix case $m \ge 1$, Algorithm CHANGEOFBASIS computes solutions to equations similar to Eqs. (6) and (8) given in the introduction. These are matrix generalizations of the Padé approximation problems; their solutions provide respectively a basis of relations $R$ of $\mathcal{M}_m^{(\gamma,f)}$ and $\tilde{\alpha}$.

In more details, Steps 2 and 3 first compute the power series expansions involved in Eqs. (6) and (8), which amounts to a type of generalized power projections. Then both approximation problems are solved at once using shifted approximant bases:

* The choice of the first $2m$ columns of $F = (S - I_m s)$ at Step 4, which are the same as in Step 2 of Algorithm CANDIDATEBASIS, and the use of a corresponding “zero shift” (first $2m$ entries of the tuple $t$ at Step 4), make this equivalent to the computation in Section 5.3 (compare Steps 3 to 5 of Algorithm CHANGEOFBASIS to Steps 1 to 3 of Algorithm CANDIDATEBASIS). This yields a basis $R$ of $\mathcal{M}_m^{(\gamma,f)}$.

* Equation (8) is solved thanks to an additional series expansion in $F$ (its last column), and the use of a sufficiently large shift (the last entry $2d$ of the tuple $t$). This yields a bivariate polynomial $\tilde{\alpha}$ which is the remainder of the requested $\alpha$ in the division by $R$.

Finally, this sought $\alpha$ can be obtained by reversing this division, using a Hermite normal form computation which also provides the minimal polynomial $\mu_{\gamma}$ (Steps 6 and 7).

The assumptions in Algorithm CHANGEOFBASIS yield a slightly stronger statement in Proposition 6.1 than in Proposition 5.6 for Algorithm CANDIDATEBASIS. Indeed, we suppose that $\gamma$ is such that $\deg(\mu_{\gamma}) = n$, whereas we make no such assumption in Algorithm CANDIDATEBASIS. From the module properties in Proposition 4.1, we deduce that $\deg(\mu_{\gamma}) = n$ implies $v_m^{(\gamma,f)} = n$, which allows us to certify the basis of relations $R$ when FAIL is not returned.

Algorithm CHANGEOFBASIS may still return FAIL; Section 8 shows that when it is called with a random $\gamma$, then with high probability, it does not fail, at least under some assumptions on $f$.

**PROPOSITION 6.1.** Given $f \in \mathbb{K}[x]$ of degree $n$ with $f(0) \neq 0$, $\gamma$ and $a$ in $\mathbb{K}[x]_{<n}$, $m \le n$ and $d$ in $\mathbb{N}_{>0}$, Algorithm CHANGEOFBASIS uses $\tilde{O}(m^{\omega}d + c(n, m, d))$ operations in $\mathbb{K}$, with $c(\cdot)$ from Eq. (2),
<caption>**Algorithm 6.1 CHANGE OF BASIS (f, γ, a, m, d)**</caption>

**Input:** $f$ of degree $n$ in $\mathbb{K}[x]$, with $f(0) \neq 0$, $\gamma \in \mathbb{K}[x]_{<n}$, $a \in \mathbb{K}[x]_{<n}$, $m \le n$ and $d \in \mathbb{N}_{>0}$

**Output:** either `FAIL` or $(R, \mu, \alpha)$ where $R \in \mathbb{K}[y]^{m \times m}_{\le 2d}$ is the Popov basis of $M_m^{(\gamma,f)}$, $\mu$ is the minimal polynomial of $\gamma$ in $\mathbb{K}[x]/\langle f \rangle$ and has degree $n$, and $\alpha \in \mathbb{K}[y]_{<n}$ with $\alpha(\gamma) \equiv a \mod f$

1. **if** $\gcd(\gamma, f) \neq 1$ **then return** `FAIL`

2. ▶ Truncated expansion of $-X^T(yI_n - M_\gamma)^{-1}v_a$ using Algorithm 3.7, $v_a \in \mathbb{K}^n$ is the coefficient vector of a
   $(r_k)_{0 \le k < 2d} \leftarrow \text{TRUNCATED POWERS}(f, \gamma^{-1} \mod f, \gamma^{-1} a \mod f, m, 2d)$
   $s \in \mathbb{K}[y]^m \leftarrow \sum_{0 \le k < 2d} s_k y^k$ where $s_k \in \mathbb{K}^m$ is the coefficient vector of $r_k$

3. ▶ Truncated expansion of $X^T(yI_n - M_\gamma)^{-1}X$ using Algorithm 3.8 (analogous to Step 1 of Algorithm 5.1)
   $(\Gamma_{i,k})_{\substack{0 \le i < m \\ 0 \le k < 2d + 2}} \leftarrow \text{BLOCK TRUNCATED POWERS}(f, \gamma^{-1} \mod f, m, 2(d+1))$
   $S_{i,k} \in \mathbb{K}^m \leftarrow$ vector of coefficients of $-\Gamma_{i,k+1} \in \mathbb{K}[x]_{<m}$, for $0 \le i < m$ and $0 \le k < 2d$
   $S \in \mathbb{K}[y]^{m \times m}_{<2d} \leftarrow \sum_{0 \le k < 2d} S_k y^k$ where $S_k = (S_{0,k} \ \dots \ S_{m-1,k}) \in \mathbb{K}^{m \times m}$

4. ▶ Fraction reconstruction using [27, 38] (analogous to Step 2 of Algorithm 5.1)

$$F \in \mathbb{K}[y]^{m \times (2m+1)}_{<2d} \leftarrow (S(y) - I_m \quad s(y))$$
$$t \in \mathbb{N}^{2m+1} \leftarrow (0, \dots, 0, 2d)$$

$$\bar{P} \in \mathbb{K}[y]^{(2m+1) \times (2m+1)}_{\le 2d} \leftarrow \text{POPOV-PM-BASIS}(F^T, 2d, t)^T, \text{ with } \bar{P} \text{ in } t\text{-shifted Popov normal form}$$

$$P \leftarrow \bar{P}_{1..2m,1..2m}$$
$$v_{\bar{\alpha}} \in \mathbb{K}[y]^{m}_{<\deg(R)} \leftarrow \bar{P}_{1..m,2m+1}$$

▷ represents $\bar{\alpha}(x, y)$, expected to satisfy $\bar{\alpha}(x, \gamma) \equiv a \mod f$

5. ▶ Ensure $R$ is a basis of $M_m^{(\gamma,f)}$, from Item (ii) of Proposition 5.4 (analogous to Step 3 of Algorithm 5.1)
   $R \leftarrow P_{1..m,1..m}$

if the sum of diagonal degrees of $R$ is less than $n$

or among the $m$ rightmost columns of $P$, one has degree $< \deg(R)$ **then return** `FAIL`

6. ▶ Compute $\mu_\gamma$, and ensure it has degree $n$

$$T \in \mathbb{K}[y]^{m \times m} \leftarrow \text{Hermite normal form of } R$$

$$\mu \in \mathbb{K}[y] \leftarrow T_{1,1}; \text{ if } \deg(\mu) < n \text{ then return } \text{FAIL}$$

▷ using [49, Algo. 1 and 3]

7. ▶ Deduce $\alpha$ and return

$$\alpha \in \mathbb{K}[y]_{<n} \leftarrow \bar{\alpha}_1 - (T_{1,2}\bar{\alpha}_2 + \dots + T_{1,m}\bar{\alpha}_m) \text{ rem } \mu, \text{ where } v_{\bar{\alpha}} = (\bar{\alpha}_1 \ \dots \ \bar{\alpha}_m)$$
**return** $(R, \mu, \alpha)$

to return either `FAIL` or $(R, \mu, \alpha)$ where $R \in \mathbb{K}[y]^{m \times m}_{\le 2d}$ is the Popov basis of $M_m^{(\gamma,f)}$, $\mu$ is the minimal polynomial $\mu_\gamma$ of $\gamma \mod f$ and has degree $n$, and $\alpha$ is the unique polynomial in $\mathbb{K}[y]_{<n}$ such that $\alpha(\gamma) \equiv a \mod f$.

If $\gcd(\gamma, f) = 1$, $v_{m,m}^{(\gamma,f)} = v_m^{(\gamma,f)}$, $\deg(\mu_\gamma) = n$ and the fraction $H = X^T(yI_n - M_\gamma)^{-1}X$ is describable in degree $d$, then the output is not `FAIL`; in that case we further have $\deg(R) \le d$.

**PROOF.** We start by showing that if the algorithm does not fail, then the truncated module $M_{m,m}^{(\gamma,f)}$ and the module of relations $M_m^{(\gamma,f)}$ are equal, and $R$ is a basis of this module.

**Steps 1 to 4:** the approximant basis $\bar{P}$. If the test at Step 1 does not fail then the specifications for Steps 2 and 3 are met. At Step 2, Algorithm `TRUNCATED POWERS` returns $r_k = [a \gamma^{-k-1} \text{ rem } f]^m_{0-1}$ for $0 \le k < 2d$ using $\tilde{O}(c(n, m, d))$ operations, according to Proposition 3.6. Thus the coefficient vector $s_k \in \mathbb{K}^m$ of $r_k$ is $X^T M_\gamma^{-k-1} v_a$, where $v_a \in \mathbb{K}^n$ is the coefficient vector of $a$, so that the polynomial vector $s(y)$ computed at Step 2 is the power series expansion of $-X^T(yI_n - M_\gamma)^{-1}v_a$ truncated at order $2d$. From Proposition 3.7, the computation of $S(y)$ at Step 3 uses $\tilde{O}(m^2 d + c(n, m, d))$ operations; $S$ is the power series expansion of $X^T(yI_n - M_\gamma)^{-1}X$ truncated at order $2d$.
Step 4 computes the t-shifted Popov approximant basis $\bar{P}$ for $F = (S - I_m s)$ at order $2d$, which uses $\tilde{O}(m^\omega d)$ operations [27, Thm. 2.4; 38, Sec. 3]. Writing

$$ \bar{P} = \begin{pmatrix} P & u \\ z & \lambda \end{pmatrix} \quad \text{for some } P \in \mathbb{K}[y]_{\le 2d}^{(2m)\times(2m)}, \lambda \in \mathbb{K}[y]_{\le 2d}, u \in \mathbb{K}[y]_{\le 2d}^{2m}, \text{ and } z \in \mathbb{K}[y]_{\le 2d}^{1\times(2m)}, $$

the fact that $\bar{P}$ is t-shifted Popov and the choice $t = (0, \dots, 0, 2d)$ ensure that $P$ is (non-shifted) Popov, that $\lambda \neq 0$, and that $\deg(z) + 2d < \deg(P) \le 2d$. The latter degree bound yields $z = 0$, hence $P$ is the Popov approximant basis of $(S - I_m)$ at order $2d$. The fact that $\bar{P}$ is t-shifted Popov also ensures that the (unique) remainder in the division of $u$ by $P$ is $u$ itself, and that the $i$th entry of $u$ has degree less than the $i$th diagonal degree of $P$.

After Step 5, $R$ is a basis of $M_{m,m}^{(y,f)} = M_m^{(y,f)}$. Let $R$ be the $m \times m$ leading principal submatrix of $P$ (and of $\bar{P}$), and let $v_{\tilde{\alpha}} \in \mathbb{K}[y]^m$ be the length-$m$ top sub-vector of $u$. Similarly to the above, $R$ is (non-shifted) Popov, and the $i$th entry of $v_{\tilde{\alpha}}$ has degree less than the $i$th diagonal degree of $R$; in particular, $\deg(v_{\tilde{\alpha}}) < \deg(R)$. Considering $H(y) = X^\top(yI_n - M_\gamma)^{-1}X$, recall from Proposition 5.1 that $D(H) = M_{m,m}^{(y,f)}$, and recall that $M_{m,m}^{(y,f)} \supseteq M_m^{(y,f)}$ with equality if and only if $v_{m,m}^{(y,f)} = v_m^{(y,f)}$. In particular, bases of $D(H)$ have determinantal degree $v_{m,m}^{(y,f)} \le v_m^{(y,f)} \le n$. Applying Item (i) of Proposition 5.4 to $H$ and the approximant basis $P$ shows that the sum of diagonal degrees of $R$ is $\deg \det(R)$, and is at most $v_{m,m}^{(y,f)}$.

As a result, if Step 5 does not return FAIL, then $n \le \deg \det(R) \le v_{m,m}^{(y,f)}$, hence $v_{m,m}^{(y,f)} = v_m^{(y,f)} = n$ and $D(H) = M_{m,m}^{(y,f)} = M_m^{(y,f)}$. Furthermore, Item (ii) of Proposition 5.4 shows that $R$ is a basis of $M_m^{(y,f)}$.

After Step 6, $\mu$ is $\mu_\gamma$ and has degree $n$. Using $\tilde{O}(m^\omega d)$ operations [49], Step 6 finds the Hermite normal form $T$ of $R$. Since $T$ is a basis of $M_m^{(y,f)}$ in upper triangular form, Proposition 4.1 states that its first diagonal entry is the minimal polynomial of $\gamma$ in $\mathbb{K}[x]/\langle f \rangle$. Hence $\mu$ computed at Step 6 is this minimal polynomial. It has degree at most $n$, and the algorithm returns FAIL at this step if and only if $\deg(\mu) < n$.

After Step 6, $v_{\tilde{\alpha}}$ represents $\tilde{\alpha}(x, y)$ such that $\tilde{\alpha}(x, \gamma) \equiv a \mod f$. Let $\tilde{\alpha} \in \mathbb{K}[x, y]_{< (m, \deg(R))}$ be the polynomial whose coefficient vector is $v_{\tilde{\alpha}}$, that is, $\tilde{\alpha} = \tilde{\alpha}_1(y) + x\tilde{\alpha}_2(y) + \dots + x^{m-1}\tilde{\alpha}_m(y)$ using notation from Step 7.

The fact that $\mu = \mu_\gamma$ has degree $n$ also ensures that there exists a unique $\alpha \in \mathbb{K}[y]_{<n}$ such that $\alpha(\gamma) \equiv a \mod f$. Then let $v_\alpha = (\alpha \ 0 \ \dots \ 0)^\top \in \mathbb{K}[y]^m$, and let $v_{\tilde{\alpha}} \in \mathbb{K}[y]^m$ be the unique remainder in the division of $v_\alpha$ by $R$. The entries of $v_{\tilde{\alpha}}$ have degree strictly less than that of the corresponding row of $R$: the degree of the $i$th entry of $v_{\tilde{\alpha}}$ is less than the $i$th diagonal degree of $R$. We also define $\tilde{\alpha} \in \mathbb{K}[x, y]_{< (m, \deg(R))}$ as the polynomial whose coefficient vector is $v_{\tilde{\alpha}}$; in particular $\tilde{\alpha}(x, \gamma) = \alpha(\gamma) \equiv a \mod f$.

We now show that $\tilde{\alpha} = \tilde{\alpha}$, which yields $\tilde{\alpha}(x, \gamma) = a \mod f$. By construction, $\tilde{\alpha}(x, y) - a(x)$ is in $M_n^{(y,f)}$, and since $yI_n - M_\gamma$ is a basis of $M_n^{(y,f)}$ (see Proposition 4.1) there is a vector $v \in \mathbb{K}[y]^n$ such that $(yI_n - M_\gamma)v = Xv_{\tilde{\alpha}} - v_\alpha$. Applying the predictable degree property [39, Thm. 6.3-13, p. 387] to the column reduced matrix $yI_n - M_\gamma$, all of whose columns have degree 1, we obtain that $\deg(v) + 1 = \deg(Xv_{\tilde{\alpha}} - v_\alpha) = \deg(v_{\tilde{\alpha}})$. Furthermore from $(yI_n - M_\gamma)v = Xv_{\tilde{\alpha}} - v_\alpha$ we get

$$ X^\top(yI_n - M_\gamma)^{-1}Xv_{\tilde{\alpha}} - X^\top v - X^\top(yI_n - M_\gamma)^{-1}v_\alpha = 0, $$
and considering truncated power series it follows that $Fq = (S - I_m s)q \equiv 0 \mod y^{2d}$, where

$$q = \begin{pmatrix} \tilde{u} \\ 1 \end{pmatrix} \in \mathbb{K}[y]^{2m+1} \quad \text{and} \quad \tilde{u} = \begin{pmatrix} v_{\tilde{\alpha}} \\ X^T v \end{pmatrix} \in \mathbb{K}[y]^{2m}.$$

Therefore $q$ is a right multiple of the approximant basis $\bar{P} = \begin{pmatrix} P \\ 0 \end{pmatrix}_{\lambda}$, which shows that $\lambda$ is a nonzero element of $\mathbb{K}$, and that $\tilde{u} - u$ is a right multiple of $P$. We check finally that the remainder of $\tilde{u}$ in the division by $P$—which is $u$ by construction—is $\tilde{u}$ itself. Indeed the degree of the $i$th entry of $v_{\tilde{\alpha}}$ is less than the $i$th diagonal degree of $R$, which is the $i$th diagonal degree of $P$; and as seen above all entries of $X^T v$ have degree at most $\deg(v_{\tilde{\alpha}}) - 1 < \deg(R) - 1$, with $\deg(R)$ being itself at most the $i$th diagonal degree of $P$ for $m+1 \le i \le 2m$. In particular $v_{\tilde{\alpha}} = v_{\tilde{\alpha}}$, hence $\tilde{\alpha} = \tilde{\alpha}$ and $\tilde{\alpha}(x, \gamma) = a \mod f$.

Step 7 computes $\alpha \in \mathbb{K}[x]_{<n}$ such that $\alpha(\gamma) = a \mod f$. Since $\deg(\mu) = n$, the Hermite normal form of $R$ has the shape $T = \begin{pmatrix} \mu \\ 0 \end{pmatrix}_{I_{m-1}}^{T_{1,*}}$, with $T_{1,*} = (T_{1,2} \cdots T_{1,m})$ and $\deg(T_{1,j}) < \deg(\mu) = n$ for $2 \le j \le m$. Then, the polynomial $\alpha = \tilde{\alpha}_1 - (T_{1,2}\tilde{\alpha}_2 + \cdots + T_{1,m}\tilde{\alpha}_m)$ rem $\mu$ constructed at Step 7 has degree less than $n$ and, by construction as well, the vector $v_{\tilde{\alpha}} = (\alpha \ 0 \ \cdots \ 0)^T \in \mathbb{K}[y]^m$ is such that $v_{\tilde{\alpha}} - v_{\tilde{\alpha}}$ is a right multiple of $T$. (In fact, $v_{\tilde{\alpha}}$ is the unique remainder in the division of $v_{\tilde{\alpha}}$ by $T$.) In particular, $v_{\tilde{\alpha}} - v_{\tilde{\alpha}}$ is a right multiple of $R$, meaning that $v_{\tilde{\alpha}}$ is equal to $v_{\tilde{\alpha}}$ modulo relations of $M_m^{(\gamma,f)}$, which implies $\alpha(\gamma) = \tilde{\alpha}(x, \gamma) \equiv a \mod f$. The computation of $\alpha$ costs $\tilde{O}(nm)$ operations in $\mathbb{K}$.

This concludes the proof of the properties of $(R, \mu, \alpha)$ in the case where the algorithm does not return FAIL. Furthermore, adding the above costs yields the cost bound claimed in the lemma, which therefore holds in general since the cost can only be smaller when the algorithm returns FAIL.

*Proof of the last claim.* The assumption $\gcd(\gamma, f) = 1$ ensures that Step 1 does not return FAIL, in which case we have seen that $P$ is a weak Popov approximant basis of $(S - I_m)$ at order $2d$.

From $\deg(\mu_\gamma) = n$ and Proposition 4.1 we know that $v_m^{(\gamma,f)} = n$, hence with the assumption $v_m^{(\gamma,f)} = v_m^{(\gamma,f)}$ we have $v_m^{(\gamma,f)} = v_m^{(\gamma,f)} = n$. Using $\mathcal{D}(H) = M_m^{(\gamma,f)} = M_m^{(\gamma,f)}$, and Item (iii) of Proposition 5.4 thanks to the assumption on $H(y) = X^T(yI_n - M_\gamma)^{-1}X$, we deduce that the $m$ rightmost columns of $P$ have degree at least $\deg(R)$ and that $R$ is a basis of $M_m^{(\gamma,f)}$ with $\deg(R) \le d$. In particular $\deg \det(R) = n$, and it follows that Step 5 does not return FAIL.

Then, the assumption on the degree of the minimal polynomial also ensures, using Proposition 4.1 as above, that the first diagonal entry of the Hermite normal form $T$ of $R$ is $\mu_\gamma$, and is the polynomial $\mu$ computed at Step 6. Therefore Step 6 does not return FAIL either: we have proved that, under the assumptions $\gcd(\gamma, f) = 1$, $v_{m,m}^{(\gamma,f)} = v_m^{(\gamma,f)}$, $\deg(\mu_\gamma) = n$, and $H = X^T(yI_n - M_\gamma)^{-1}X$ is describable in degree $d$, then the output is not FAIL and $\deg(R) \le d$. $\square$

*Notes.* Shoup's algorithm for computing $\alpha$ in the case $m = 1$ uses only $n$ terms of the sequence $(\ell(\gamma^k a))_k$, or more generally $d$ terms, where $d$ is a known bound on $\deg(\mu_a)$. Here as well, if one knows that the sought basis of relations satisfies $\deg(R) \le d$, for example under the conditions of Proposition 6.1 ensuring success, then the algorithm may be modified so as to require only $d$ terms of the expansion of $-X^T(yI_n - M_\gamma)^{-1}v_\alpha$ instead of $2d$. The vector $v_{\tilde{\alpha}}$ would appear in the approximant basis at order $d$, and from there one would consider a residual approximant problem focusing on obtaining the missing part of $R$. This is not detailed here, as this would complicate the presentation without bringing an improvement to the asymptotic complexity.

Modular composition and inverse composition are very similar. They both involve the computation of a matrix of relations and use symmetric steps with similar complexities. Indeed, the division
with remainder of Section 4.2.1 is used in both algorithms to change between univariate and bivariate representations efficiently. Also, the application of Algorithm BIVARIATEMODULARCOMPOSITION at the last step of composition in Algorithm BIVARIATEMODULARCOMPOSITIONWITHRELATION-MATRIX is reflected by TRUNCATEDPOWERS as starting step of inverse composition in Algorithm CHANGEOFBASIS. Both these steps have cost $\tilde{O}(c(n, m, d))$ from Propositions 3.4 and 3.6, respectively (see also Section 3.4.3).

# 7 THE BLOCK HANKEL MATRIX $H_k^{(a,f)}$ AND ITS GENERIC PROPERTIES

Matrices of relations are obtained either by Algorithm MATRIXOFRELATIONS directly, or by Algorithm CHANGEOFBASIS after a change of basis. In both cases, for the correctness of the computation to be granted via Propositions 5.8 and 6.1, we need $v_{m,m}^{(a,f)}$ and $v_m^{(a,f)}$ to be equal (and, equivalently, $M_{m,m}^{(a,f)} = M_m^{(a,f)}$) and the fraction $H(y) = X^T(yI_n - M_a)^{-1}X$ to be describable in degree $d$, or the same statement with $\gamma$ in place of $a$. It is thus important to understand when these properties hold.

Recall from Section 3.4.1 the matrices $K_{m,d}^{(a,f)}$ and $L_{m,d}^{(a,f)}$, that are defined for $m \in \{1, \dots, n\}$ by

$$ K_{m,d}^{(a,f)} = (X \quad \dots \quad M_a^{d-1}X) \in \mathbb{K}^{n \times (md)} \quad \text{and} \quad L_{m,d}^{(a,f)} = \begin{pmatrix} X^T \\ \vdots \\ X^T M_a^{d-1} \end{pmatrix} \in \mathbb{K}^{(md) \times n}, $$

and that correspond to Algorithms BIVARIATEMODULARCOMPOSITION and TRUNCATEDPOWERS respectively, and also to the maps $\kappa_{m,d}^{(a,f)}$ and $\lambda_{m,d}^{(a,f)}$. Their product forms the block Hankel matrix

$$ H_k^{(a,f)} = L_{m,d}^{(a,f)} K_{m,d}^{(a,f)} = \begin{pmatrix} H_0 & H_1 & \dots & H_{d-1} \\ H_1 & \ddots & \ddots & H_d \\ \vdots & \ddots & \ddots & \vdots \\ H_{d-1} & H_d & \dots & H_{2d-2} \end{pmatrix} \in \mathbb{K}^{(md) \times (md)}, \quad (20) $$

with $H_k = X^T M_a^k X$ for $k$ in $\mathbb{N}$. This matrix, and in particular its rank, is strongly related to the two properties mentioned above [70; 46, p. 97].

The outcomes of this section are the following. For any positive parameters $m \le n$ and $d$, as soon as $\text{rank}(H_k^{(a,f)}) = v_m^{(a,f)}$, then $v_{m,m}^{(a,f)} = v_m^{(a,f)}$ and $H$ is describable in degree $d$ (Section 7.1). This happens in particular when $f(0) \ne 0$, $d \ge \lceil n/m \rceil$ and either $\text{deg}(a) = m$ (Section 7.2) or for a generic choice of $a$ (Section 7.3). Also, for generic choices of the roots of $f$ and of the values of $a$ at these roots, $\text{rank}(H_k^{(a,f)}) = v_m^{(a,f)}$ as soon as $d \ge \lceil v_m^{(a,f)}/m \rceil$ (Section 7.4). As in previous sections, notation such as $v_m^{(a,f)}$, $H_k^{(a,f)}$, $\lambda_{m,d}^{(a,f)}$, etc. is often shortened into $v_m$, $H_{k,m,d}$, $\lambda_{m,d}$, etc.

These results will be used in Section 8 for the analysis of the randomized composition algorithm when $f$ is separable (Section 8.3), or when $f$ is purely inseparable, which includes the case of power series composition (Sections 8.4 and 8.5).

## 7.1 Relation between block Hankel matrix rank and fraction description degree

The key condition to control the degrees of fraction descriptions of $H(y)$ and obtain matrices of relations is the equality

$$ \text{rank}(H_k) = v_m. $$

The special case when $\text{rank}(H_k) = n$ is common, and appears naturally later on. The proof of the following result relies in an essential manner on Lemma 7.2, which we give next (the references we cite only give a sketch of proof).
**PROPOSITION 7.1.** Given $f \in \mathbb{K}[x]$ of degree $n$, $a \in \mathbb{K}[x]_{<n}$, and positive integers $m \le n$ and $d$, the rank of $\mathrm{Hk}_{m,d}^{(a,f)}$ is at most $v_m^{(a,f)}$. In case of equality, we have $v_{m,m}^{(a,f)} = v_m^{(a,f)}$ and $H(y) = X^T(yI_n - M_a)^{-1}X$ is describable in degree $d$.

In particular, if $\mathrm{Hk}_{m,d}^{(a,f)}$ has rank $n$, then $v_{m,m}^{(a,f)} = v_m^{(a,f)} = n$ and $H(y)$ is describable in degree $d$.

**PROOF.** Using Proposition 4.1 the inclusion $\mathcal{M}_m \subseteq \mathcal{M}_{m,m}$ implies $v_{m,m} \le v_m \le n$, so that by Lemma 7.2 below we have rank($\mathrm{Hk}_{m,d}$) $\le v_{m,m} \le v_m \le n$. If $\mathrm{Hk}_{m,d}$ has rank $v_m$, then $v_m = v_{m,m}$, and the claim on $H$ follows again from Lemma 7.2. The case where the rank is $n$ follows similarly. $\square$

**LEMMA 7.2** ([46, SEC. 2.1] AND [71, LEM. 2.4]). For positive integers $m \le n$ and $d$, the rank of $\mathrm{Hk}_{m,d}^{(a,f)}$ is at most $v_{m,m}^{(a,f)}$, with equality if and only if $H(y) = X^T(yI_n - M_a)^{-1}X$ is describable in degree $d$.

**PROOF.** We denote by $H_k = X^T M_a^k X \in \mathbb{K}^{m \times m}$ the coefficient in the expansion of $H$ at infinity:

$$H(y) = X^T(yI_n - M_a)^{-1}X = \sum_{k \ge 0} H_k y^{-k-1} = \sum_{k \ge 0} X^T M_a^k X y^{-k-1}.$$

To show that the rank is at most $v_{m,m}$ we first note that $\mathrm{Hk}_{m,d}$ is a submatrix of $\mathrm{Hk}_{m,d+1}$ for $d \ge 0$, the sequence $(\mathrm{rank}(\mathrm{Hk}_{m,d}))_{d \ge 0}$ is thus nondecreasing. Since $\mathcal{M}_{m,m}$ is the module of vector generators for the sequence $\{H_k\}_{k \ge 0}$ (Section 5.1), the minimal generating polynomial $F \in \mathbb{K}[y]^{m \times m}$ in Popov form for that sequence is a basis of $\mathcal{M}_{m,m}$ ([70, Def. 2.5] and [46, Def. 2.3]). It follows that $\mathrm{deg} \det(F) = v_{m,m}$, and [46, Eq. (2.6)] shows that for $d \ge n$, the rank of $\mathrm{Hk}_{m,d}$ is $v_{m,m}$. So the first claim is proved.

If the rank of $\mathrm{Hk}_{m,d}$ is equal to $v_{m,m}$, then this rank is also that of the infinite matrix corresponding to the system (see also Eq. (16))

$$H_k v_0 + \cdots + H_{k+d} v_d = 0 \quad \text{for } k \ge 0, \qquad (21)$$

thus a solution to

$$H_k v_0 + \cdots + H_{k+d} v_d = 0 \quad \text{for } 0 \le k \le d-1 \qquad (22)$$

is also a solution to Eq. (21). Since the rank of $\mathrm{Hk}_{m,d}$ is maximal we also know that the last block column of $\mathrm{Hk}_{m,d+1}$ is a linear combination of the previous ones. This provides with $m$ linearly independent $R_1, \dots, R_m \in \mathbb{K}[y]^m$, of degree $d$, whose coefficient vectors in $y$ are solutions to Eq. (22), hence to Eq. (21). Let $R$ be the matrix in $\mathbb{K}[y]^{m \times m}$ whose $j$th column is $R_j$. Using Eq. (21) we deduce that $HR = Q$ with $Q \in \mathbb{K}[y]^{m \times m}$ (see also Eq. (17)). This gives a right fraction description $H = QR^{-1}$ (which may not be irreducible) with denominator of degree $d$. The same reasoning on the left side gives a left matrix description of degree $d$, hence $H$ is describable in degree $d$.

Conversely, a right matrix description $H = QR^{-1}$ with $R$ of degree at most $d$ gives $R_j$'s whose coefficient vectors are solutions to Eq. (21). Since $F$ is a basis of the module of vector generators for $\{H_k\}_{k \ge 0}$, $R$ must be a multiple of $F$. By minimality $F$ has degree at most $d$ [39, Thm. 6.5-10, p. 458], and using [46, Eq. (2.6)] the rank of the infinite block Hankel matrix restricted to its first $d$ block columns is maximal. Starting from a left description, in an analogous way we obtain that the rank restricted to the first $d$ block rows is maximal, which yields that $\mathrm{Hk}_{m,d}$ has rank $\mathrm{deg} \det(F) = v_{m,m}$. $\square$

## 7.2 Families with $\mathrm{Hk}_{m,d}^{(a,f)}$ of rank $n$

A simple condition implies the equality $\mathrm{rank}(\mathrm{Hk}_{m,d}) = n$ of Proposition 7.1.

**PROPOSITION 7.3.** Let $f \in \mathbb{K}[x]$ have degree $n$, let $a \in \mathbb{K}[x]_{<n}$, and let $m$ be a positive integer. If $f(0) \ne 0$ and $\mathrm{deg}(a) = m$ (hence $1 \le m < n$), then the block Hankel matrix $\mathrm{Hk}_{m,d}^{(a,f)} \in \mathbb{K}^{(md) \times (md)}$ has rank $n$ for all $d \ge [n/m]$.
The rest of this subsection is devoted to the proof of this result. It is a basis for the genericity result in the next subsection.

PROOF. Given $c$ in $K$, any minor of $Hk_{m,d}^{(ca,f)}$ is equal to the corresponding minor of $Hk_{m,d} = Hk_{m,d}^{(a,f)}$ times a power of $c$. It follows that $\mathrm{rank}(Hk_{m,d}^{(ca,f)}) = \mathrm{rank}(Hk_{m,d})$ for any $c \neq 0$, and therefore in the rest of the proof we can assume that $a$ is monic of degree $m$.

By Eq. (20), it is sufficient to show that the mappings $\kappa_{m,d}$ and $\lambda_{m,d}$ associated to $K_{m,d}$ and $L_{m,d}$ are surjective and injective, respectively.

The mapping $\kappa_{m,d}$ is surjective. By assumption, $n \leq md$ so that surjectivity of $\kappa_{m,d}$ is equivalent to the matrix $K_{m,d} \in K^{n \times (md)}$ from Eq. (11) having full row rank $n$. Indeed, the first $n$ columns of $K_{m,d}$ are the coefficients of the family of polynomials $x^i a^j$ rem $f$, for $0 \leq i < m$ and $0 \leq j < d$, with $0 \leq i + jm < n$. Since $\mathrm{deg}(a) = m$, these columns form an upper triangular matrix, with 1's on the diagonal; this proves the claim.

The mapping $\lambda_{m,d}$ is injective. Equivalently, we have to show that $L_{m,d}$ has full column rank $n$. This follows from the structure of this matrix, seen at the level of polynomials.

LEMMA 7.4. With the notation and hypotheses of Proposition 7.3, let

$$p_i = [ax^{n-m+i} \operatorname{rem} f]_0^{m-1}, \quad i = 0, \dots, m-1.$$

Then,

(i) if $m \leq n/2$, the $m$ polynomials $p_0, \dots, p_{m-1}$ are linearly independent;

(ii) if $n/2 < m$, the $n-m$ polynomials $p_{2m-n}, \dots, p_{m-1}$ are linearly independent.

PROOF. The two cases require different proofs, sharing common ingredients. For $i \geq 0$, let $r_i = x^{n+i} \operatorname{rem} f$. For $b$ in $K[x]_{<n}$, we then have

$$x^i b \operatorname{rem} f = [x^i b]_0^{m-1} + \delta_{b,i}, \qquad (23)$$

for some $\delta_{b,i}$ in $\mathrm{Span}(r_0, \dots, r_{i-1})$, in particular $\delta_{b,0} = 0$. Applying this to $b = r_0 = x^n \operatorname{rem} f$ yields $r_i = [x^i r_0]_0^{m-1} + \delta_{r_0,i}$. Taking this relation modulo $x^m$ gives $[r_i]_0^{m-1} = [x^i r_0]_0^{m-1} + \mu_i$, with $\mu_i$ in $\mathrm{Span}([r_0]_0^{m-1}, \dots, [r_{i-1}]_0^{m-1})$. As a result, for $i \geq 0$, we get

$$\mathrm{Span}([r_0]_0^{m-1}, \dots, [r_i]_0^{m-1}) = \mathrm{Span}([r_0]_0^{m-1}, \dots, [x^i r_0]_0^{m-1}).$$

Writing $f = f_0 + \dots + f_{n-1}x^{n-1} + x^n$, we get $r_0 = -f_0 - f_1x - \dots - f_{n-1}x^{n-1}$. Since $f_0 \neq 0$ by assumption, $[x^i r_0]_0^{m-1}$ has valuation $i$ for $0 \leq i < m$; this implies that $\mathrm{Span}([r_0]_0^{m-1}, \dots, [r_i]_0^{m-1})$ has dimension $i+1$ for $0 \leq i < m$.

*Proof of Item (i).* Let $b = ax^{n-m} \operatorname{rem} f$ in Eq. (23). Upon reduction modulo $x^m$, for $0 \leq i < m$, we obtain the relation $p_i = [x^i b]_0^{m-1} + \mu'_i$, with $\mu'_i$ in $\mathrm{Span}([r_0]_0^{m-1}, \dots, [r_{i-1}]_0^{m-1})$.

Since $ax^{n-m}$ is monic of degree $n$ ($a$ has degree $m$), with valuation at least $n-m \geq m$ (here, $m \leq n/2$), we get $[b]_0^{m-1} = [r_0]_0^{m-1}$, and thus $[x^i b]_0^{m-1} = [x^i r_0]_0^{m-1}$ for $0 \leq i < m$. This gives $p_i = [r_i]_0^{m-1} + \mu'_i - \mu_i$, with $\mu'_i - \mu_i$ in $\mathrm{Span}([r_0]_0^{m-1}, \dots, [r_{i-1}]_0^{m-1})$. In particular, taking all $i$ up to $m-1$, we get the equality $\mathrm{Span}(p_0, \dots, p_{m-1}) = \mathrm{Span}([r_0]_0^{m-1}, \dots, [r_{m-1}]_0^{m-1})$, and we saw that the latter has dimension $m$. Item (i) is proved.

*Proof of Item (ii).* Assume that $q = c_m x^m + \dots + c_{n-1} x^{n-1}$ is such that $[aq \operatorname{rem} f]_0^{m-1} = 0$. We prove that all $c_i$'s vanish.

We can rewrite $aq \operatorname{rem} f$ as $x^m b \operatorname{rem} f$, with $b = a(q/x^m)$; since $a$ has degree $m$, $b$ is in $K[x]_{<n}$. Applying Eq. (23) to $b$ and $i = m$, our assumption that $[x^m b \operatorname{rem} f]_0^{m-1} = 0$ implies that $[\delta_{b,m}]_0^{m-1} =$
0. The linear independence of $[r_0]_0^{m-1}, \dots, [r_{m-1}]_0^{m-1}$ then shows that $\delta_{b,m}$ itself is zero, so that $x^m b \operatorname{rem} f = [x^m b]_0^{n-1}$.

This shows that the remainder $[x^m b]_n^{m-1}$ rem $f$ vanishes. Because $x$ is invertible modulo $f$ (since $f_0 \neq 0$), it follows that $[b]_{n-m}^{m-1}$ vanishes modulo $f$, or equivalently that $[b]_{n-m}^{m-1} = 0$. Since $a$ is monic of degree $m$, and since $n-m < m$, the definition of $b$ then implies that all coefficients $c_i$'s vanish. Hence, Item (ii) is proved. $\square$

Let now $v \in K[x]_{<n}$ be such that

$$[v]_0^{m-1} = [av \operatorname{rem} f]_0^{m-1} = \dots = [a^{d-1}v \operatorname{rem} f]_0^{m-1} = 0.$$

We prove that $\deg(v) < n - mi$ for $i = 0, \dots, d-1$. For $d = \lceil n/m \rceil$, this gives $\deg(v) < m$; together with the assumption $[v]_0^{m-1} = 0$, this proves that $v = 0$.

The proof is by induction. For $i = 0$, there is nothing to prove. If the claim holds for some index $i < d-1$, since $a$ has degree $m$, for any $w$ in $K[x]_{<n}$, the polynomial $[aw \operatorname{rem} f]_0^{m-1}$ splits into two parts:

$$[aw \operatorname{rem} f]_0^{m-1} = [a[w]_0^{n-m-1}]_0^{m-1} + [a[w]_{n-m}^{m-1} \operatorname{rem} f]_0^{m-1}.$$

Apply this identity with $w = a^i v \operatorname{rem} f$. Then, both the left-hand side and the first summand vanish: the former because $[a^{i+1}v \operatorname{rem} f]_0^{m-1} = 0$, the latter because $[a^i v \operatorname{rem} f]_0^{m-1} = 0$, i.e., $w = a^i v \operatorname{rem} f$ has valuation at least $m$. We deduce that $[a[w]_{n-m}^{m-1} \operatorname{rem} f]_0^{m-1} = 0$, with $w = a^i v \operatorname{rem} f$.

* If $m \le n/2$, the linear independence of the polynomials $p_j = [ax^{n-m+j} \operatorname{rem} f]_0^{m-1}$, for $j = 0, \dots, m-1$, then shows that $[w]_{n-m}^{n-1} = [a^i v \operatorname{rem} f]_{n-m}^{n-1}$ vanishes.

* If $m > n/2$, then the assumption that $w$ has valuation at least $m$, with thus $m > n-m$, shows that $[w]_{n-m}^{n-1} = [w]_m^{n-1}$. In this case, the linear independence of the polynomials $p_j$ for $j = 2m-n, \dots, m-1$ shows that $[w]_{n-m}^{n-1} = [w]_m^{n-1} = 0$.

In other words, in both cases, we have proved that $w = a^i v \operatorname{rem} f$ has degree less than $n-m$.

On the other hand, the induction assumption that $\deg(v) < n - mi$ implies that $a^i v \operatorname{rem} f = a^i v$, so the latter has degree less than $n-m$. Since $a^i$ has degree $mi$, this shows that $\deg(v) < n-m(i+1)$, as claimed. $\square$

## 7.3 Generic regularity in $a$ and $f$

In all this document, genericity is understood in the Zariski sense:

**Definition 7.5.** A property $\mathcal{P}$ of certain parameters $(u_1, \dots, u_s)$ holds for a generic choice of $(u_1, \dots, u_s)$ in $K^s$ if there exists a nonzero polynomial $\Delta$ in $K[\bar{u}_1, \dots, \bar{u}_s]$ (where the $\bar{u}_i$'s are new indeterminates) such that $\Delta(u_1, \dots, u_s) \neq 0$ implies that $\mathcal{P}(u_1, \dots, u_s)$ holds.

Note that if $K$ is finite, there may be no choice of the $u_i$'s in $K$ for which $\Delta$ does not vanish, but such points exist in a finite extension of $K$ of sufficiently large degree (such as $O(\log(n))$) when the degree of $\Delta$ is polynomial in $n$, as is the case below).

**PROPOSITION 7.6.** Let $f$ in $K[x]$ be of degree $n$ and such that $f(0) \neq 0$. For any $m \in \{1, \dots, n\}$ there exists a nonzero polynomial $\Delta_{f,m}$ in $K[\bar{a}_0, \dots, \bar{a}_{n-1}]$ of degree at most $2n^2/m$ such that for $a = a_0 + \dots + a_{n-1}x^{n-1}$ in $K[x]_{<n}$, if $\Delta_{f,m}(a_0, \dots, a_{n-1}) \neq 0$ then $\mathrm{Hk}_{m,d}^{(a,f)} \in K^{(md) \times (md)}$ has rank $n$ for any $d \ge \lceil n/m \rceil$.

### 7.3.1 Proof of Proposition 7.6.

**LEMMA 7.7.** Let $m, n$ be positive integers, with $m \in \{1, \dots, n\}$, and let $\bar{f} = \bar{f}_0 + \dots + \bar{f}_{n-1}x^{n-1} + x^n$ and $\bar{a} = \bar{a}_0 + \dots + \bar{a}_{n-1}x^{n-1}$ be polynomials in $\mathbb{Z}[\bar{a}_0, \dots, \bar{a}_{n-1}, \bar{f}_0, \dots, \bar{f}_{n-1}][x]$. Then any $n$-minor of $\mathrm{Hk}_{m, \lceil n/m \rceil}^{(\bar{a}, \bar{f})}$ has degree at most $2n^2/m$ in $\bar{a}_0, \dots, \bar{a}_{n-1}$ and $2n^2(n-1)/m$ in $\bar{f}_0, \dots, \bar{f}_{n-1}$.
PROOF. The entries of the multiplication matrix $M_{\tilde{a}}$ are the coefficients of $x^k \tilde{a}$ rem $\tilde{f}$ for $k = 0, \dots, n-1$, and are therefore polynomials of degree 1 in the coefficients $\tilde{a}_0, \dots, \tilde{a}_{n-1}$ and at most $n-1$ in the coefficients $\tilde{f}_0, \dots, \tilde{f}_{n-1}$. In turn, the coefficients of $M_{\tilde{a}}^j$ have degree at most $j$ in $\tilde{a}_0, \dots, \tilde{a}_{n-1}$ and $j(n-1)$ in $\tilde{f}_0, \dots, \tilde{f}_{n-1}$. For $0 \le i, j < \lfloor n/m \rfloor$, the $m \times m$ block of coordinates $(i, j)$ in $\text{Hk}_{m, \lfloor n/m \rfloor}^{(\tilde{a}, \tilde{f})}$ is a submatrix of $M_{\tilde{a}}^{i+j}$; it has degree at most $i+j$ in $\tilde{a}_0, \dots, \tilde{a}_{n-1}$ and $(i+j)(n-1)$ in $\tilde{f}_0, \dots, \tilde{f}_{n-1}$. As a result, any $n$-minor of this matrix has degree at most $m \lfloor n/m \rfloor (\lfloor n/m \rfloor - 1) \le 2n^2/m$ in $\tilde{a}_0, \dots, \tilde{a}_{n-1}$ and $m \lfloor n/m \rfloor (\lfloor n/m \rfloor - 1)(n-1) \le 2n^2(n-1)/m$ in $\tilde{f}_0, \dots, \tilde{f}_{n-1}$. $\square$

Take $f$ of degree $n$ with $f(0) \neq 0$. Proposition 7.3 with $a = x^m$ shows that at least one $n$-minor of $\text{Hk}_{m, \lfloor n/m \rfloor}^{(x^m, f)}$ is nonzero, so the corresponding $n$-minor of $\text{Hk}_{m, \lfloor n/m \rfloor}^{(\tilde{a}, f)}$ is not identically zero. We take this minor for $\Delta_{f,m}$, and its degree is then bounded by Lemma 7.7.

### 7.3.2 Note: basis of relations for a generic $a$
For any $f$ in $K[x]$ with $f(0) \neq 0$, and for a generic $a$ in $K[x]_{<n}$, Proposition 7.6 shows that the rank of $\text{Hk}_{m,d}^{(a,f)}$ is $n$, with $d = \lfloor n/m \rfloor$. From Proposition 7.1 we then obtain $v_m^{(a,f)} = v_{m,m}^{(a,f)} = n$ and the describability of $H$ in degree $d$. Therefore, by Proposition 5.6, Algorithm CANDIDATEBASIS returns a basis of $M_m^{(a,f)}$ and the flag CERT.

## 7.4 Generic rank for a separable $f$

We now study the rank of $\text{Hk}_{m,d}$, for a generic choice of the roots of $f$, and for a generic choice of the values of $a$ at these roots, subject to certain combinatorial conditions.

### 7.4.1 Definitions.
Consider pairwise distinct $\xi_1, \dots, \xi_n$ in an algebraic closure $\bar{K}$ of $K$. To such points, we associate the polynomial $f = (x - \xi_1) \cdots (x - \xi_n)$. We also consider $a \in K[x]_{<n}$, and we say that a takes values $\lambda_1, \dots, \lambda_r$ at $\xi_1, \dots, \xi_n$ with multiplicities $\ell_1, \dots, \ell_r$ if the following holds:

- $\lambda_1, \dots, \lambda_r$ are pairwise distinct elements in $\bar{K}$;

- $\ell_1 + \cdots + \ell_r = n$, with all $\ell_i$ positive integers;

- for $i = 1, \dots, r$, $a(\xi_{\sigma_i+1}) = \cdots = a(\xi_{\sigma_i+\ell_i}) = \lambda_i$, where we write $\sigma_i = \ell_1 + \cdots + \ell_{i-1}$ (the empty sum for $i=1$ is zero).

In view of our application, we also assume that the $\xi_i$'s are such that $f$ is in $K[x]$.

### 7.4.2 Generic rank.

**PROPOSITION 7.8.** Fix positive integers $m \in \{1, \dots, n\}$ and $\ell = (\ell_1, \dots, \ell_r)$ such that $\ell_1 + \cdots + \ell_r = n$. Then, there exists a nonzero polynomial $\Gamma_{\ell,m} \in \mathbb{Z}[\tilde{\xi}_1, \dots, \tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r]$ such that the following holds. For pairwise distinct nonzero $\xi_1, \dots, \xi_n$ in $\bar{K}$ such that $f = c(x - \xi_1) \cdots (x - \xi_n)$ with $c \in K \setminus \{0\}$ is in $K[x]$ and for $a \in K[x]$ that takes values $\lambda_1, \dots, \lambda_r$ at $\xi_1, \dots, \xi_n$ with multiplicities $\ell_1, \dots, \ell_r$, if $\Gamma_{\ell,m}(\tilde{\xi}_1, \dots, \tilde{\xi}_n, \lambda_1, \dots, \lambda_r)$ is nonzero, then

$$ \text{rank}(\text{Hk}_{m,d}^{(a,f)}) = v_m^{(a,f)} \quad \text{for any } d \ge \lfloor v_m^{(a,f)}/m \rfloor, $$

with in addition the equality

$$ v_m^{(a,f)} = \sum_{i=1}^{r} \min(\ell_i, m). $$

Finally, for any pairwise distinct $\lambda_1, \dots, \lambda_r$, the polynomial $\Gamma_{\ell,m}(\tilde{\xi}_1, \dots, \tilde{\xi}_n, \lambda_1, \dots, \lambda_r)$ is nonzero and has degree at most $2n^2$.
### 7.4.3 Proof of Proposition 7.8.
The rather long proof is decomposed as follows. First, the expression for the determinantal degree $v_m$ is established. For the proof of the rest of the proposition we exploit the factorization $Hk_{m,d} = L_{m,d}K_{m,d}$, that is analyzed through a series of lemmas.

The ranks of the matrices $K_{m,d}$ and $L_{m,d}$ are related to that of a simple matrix $\mathcal{P}_{\ell,m,d}$ (see Eq. (30)). This leads to the proof that for $d = \lceil v_m/m \rceil$, the rank of $K_{m,d}$ and $L_{m,d}$ is $v_m$ generically. Then we prove that generically, taking any $d_0 \ge \lceil v_m/m \rceil$ is sufficient for studying the rank of $Hk_{m,d}$. The proof is concluded by establishing that the rank is $v_m$ when $d_0$ is $r$, the number of distinct values $a(\xi_k)$'s: for this value of $d_0$, we establish that the intersection of the image of $K_{m,d}$ with the kernel of $L_{m,d}$ is reduced to 0. The polynomial $\Gamma_{\ell,m}$ and the degree bounds are derived from the proof.

Determinantal degree $v_m^{(a,f)}$. As in the proposition, let $\xi_1, \dots, \xi_n$ be pairwise distinct in $\overline{\mathbb{K}}$ and let $f = c(x - \xi_1) \cdots (x - \xi_n)$. The Lagrange interpolation polynomials

$$ \mathcal{L}_k(x) = \frac{1}{f'(\xi_k)} \frac{f(x)}{x - \xi_k} = \prod_{\ell \neq k} \frac{x - \xi_\ell}{\xi_k - \xi_\ell}, \quad k = 1, \dots, n. \qquad (24) $$

form a basis of $\overline{A} := \overline{\mathbb{K}}[x]/\langle f \rangle$. For any $a \in \overline{A}$, the matrix of multiplication by $a$ is diagonalizable, its eigenvalues are the values of $a$ at the $\xi_j$'s, and the Lagrange polynomials are eigenvectors. The characteristic polynomial $\chi_a$ of a modulo $f$ is therefore given by

$$ \chi_a = \prod_{k=1}^{n} (y - a(\xi_k)) \in \overline{\mathbb{K}}[y]. $$

For $1 \le i \le r$, we define $S_i = \{k \in \{1, \dots, n\} \mid a(\xi_k) = \lambda_i\}$ and use that

$$ S_i = \{\sigma_i + 1, \dots, \sigma_i + \ell_i\}. \qquad (25) $$

With these conventions we have the factorization

$$ \chi_a = \prod_{i=1}^{r} (y - \lambda_i)^{\ell_i}, $$

where the factors $(y - \lambda_i)$ are pairwise coprime. The Smith normal form of $yI_n - M_a$ is then known and an explicit expression for the determinantal degree $v_m$ can be given: $yI_n - M_a$ has $\max(\ell_i)$ nontrivial invariant factors; for $1 \le k \le \max(\ell_i)$, the $k$th one is $\prod_{1 \le i \le r} (y - \lambda_i)^{\epsilon_{i,k}}$, where $\epsilon_{i,k} = 1$ if $k \le \ell_i$ and 0 otherwise. From there, recalling from Eq. (7) that for $m$ in $\{1, \dots, n\}$, $v_m$ is the sum of the degrees of the first $m$ such invariant factors, we have:

$$ v_m = \sum_{k=1}^{\min(m, \max(\ell_i))} \mathrm{card}(\{i \mid \ell_i \le k\}) = \sum_{i=1}^{r} \min(\ell_i, m). \qquad (26) $$

This proves the claim regarding $v_m$ in the proposition (this claim thus holds without further assumption on the $\xi_i$'s and $\lambda_i$'s).

Maximal rank of $Hk_{m,d}^{(a,f)}$.

### LEMMA 7.9.
Let $A \in \mathbb{K}^{n \times n}$ and $m \in \mathbb{N}_{>0}$, and let $v$ be the sum of the degrees of the $\min(m, n)$ highest degree invariant factors of $yI_n - A$. Then for any collection of $m$ vectors $v_1, \dots, v_m \in \mathbb{K}^n$, one has $\dim(\mathrm{Span}(A^i v_j, 0 \le i, 1 \le j \le m)) \le v$.

### PROOF.
We let $\tilde{v} = \dim(\mathrm{Span}(A^i v_j, 0 \le i, 1 \le j \le m))$. For $1 \le j \le m$, let $d_j \ge 0$ be the first index such that $A^{d_j} v_j \in \mathrm{Span}(v_j, A v_j, \dots, A^{d_j-1} v_j, \{A^i v_k \mid 0 \le i, 0 \le k < j\})$; if $l \ge d_j$ then $A^l v_j$ also
belongs to the latter subspace of $\mathbb{K}^n$, which is therefore stable under left-multiplication by $A$. This
holds for any $1 \le j \le m$, hence $d_1 + \dots + d_m = \tilde{v}$, and the matrix

$$
P_1 = (v_1 \ Av_1 \ \dots \ Av_1^{d_1-1} \ \dots \ v_m \ Av_m \ \dots \ Av_m^{d_m-1}) \in \mathbb{K}^{n \times \tilde{v}}
$$

has rank $\tilde{v}$ and can be completed into a nonsingular matrix $P = (P_1 \ P_2) \in \mathbb{K}^{n \times n}$. By applying the
change of basis $P^{-1}AP$ we obtain

$$
P^{-1}(yI_n - A)P = \begin{pmatrix} yI_{\tilde{v}} - C & B_1 \\ 0 & yI_{n-\tilde{v}} - B_2 \end{pmatrix} \in \mathbb{K}[y]^{n \times n}, \quad (27)
$$

where $C \in \mathbb{K}^{\tilde{v} \times \tilde{v}}$, $B_1 \in \mathbb{K}^{\tilde{v} \times (n-\tilde{v})}$, $B_2 \in \mathbb{K}^{(n-\tilde{v}) \times (n-\tilde{v})}$. Thanks to the form of $P_1$, the matrix $C \in \mathbb{K}^{\tilde{v} \times \tilde{v}}$ is block upper triangular with at most $m$ companion blocks $C_j$ of dimensions $d_j$ on the diagonal (there is no block for $d_j = 0$, and at most $n$ of the $d_j$'s are nonzero). By a unimodular row transformation $U_j \in \mathbb{K}[y]^{d_j \times d_j}$, a matrix $yI_{d_j} - C_j$ can be brought into an upper triangular form $T_j(y) = U_j(y)(yI_{d_j} - C_j)$, which has diagonal entries 1 except for the last entry which is the characteristic polynomial $\chi^{(j)} = y^{d_j} - \chi_{d_j-1}^{(j)}y^{d_j-1} - \dots - \chi_0^{(j)}$ of $C_j$:

$$
\begin{pmatrix} -1 & & & \\ & \ddots & & \\ & & -1 & \\ 1 & y & \dots & y^{d_j-1} \end{pmatrix} \begin{pmatrix} yI_{d_j} - \begin{pmatrix} \chi_0^{(j)} & & & \\ & \ddots & & \\ & & \ddots & \\ & & & \chi_{d_j-1}^{(j)} \end{pmatrix} \end{pmatrix} = \begin{pmatrix} 1 & & & \\ & \ddots & & \\ & & 1 & \\ & & & \chi^{(j)}(y) \end{pmatrix} \in \mathbb{K}[y]^{d_j \times d_j},
$$

Therefore Eq. (27) can be rewritten as

$$
U(y)P^{-1}(yI_n - A)P = \begin{pmatrix} T(y) & \bar{B}_1(y) \\ 0 & yI_{n-\bar{v}} - B_2 \end{pmatrix} = \begin{pmatrix} I_{\bar{v}} & \bar{B}_1(y) \\ 0 & yI_{\bar{n}-\bar{v}} - B_2 \end{pmatrix} \begin{pmatrix} T(y) & 0 \\ 0 & I_{n-\bar{v}} \end{pmatrix}, \quad (28)
$$

where $U = \operatorname{diag}(U_1, \dots, U_m, I_{n-\tilde{v}})$ is unimodular (with no $U_j$ if $d_j = 0$), and $T \in \mathbb{K}[y]^{\tilde{v} \times \tilde{v}}$ is block upper triangular with diagonal blocks the $T_j$'s. The matrix $T$ is triangular with 1's on the diagonal except for at most $m$ entries. We deduce that the gcd of the minors of dimension $k$ of $T$ is a unit for $1 \le k \le \tilde{v} - m$, and that $T$ has at most $m$ nontrivial invariant factors [58, Ch. II, Eq. (13)]. The product of these invariant factors is $\det(T) = \prod_j \chi^{(j)}$, whose degree is $d_1 + \dots + d_m = \tilde{v}$. From the matrix product on the right-hand side of Eq. (28), these latter invariant factors divide the $m$ highest degree invariant factors of $yI_n - A$ [58, Thm. II.14]. From the definition of $\nu$ we obtain $\tilde{v} \le \nu$. $\square$

$$
\begin{equation}
\begin{gathered}
\text{With } A = M_a \text{ or } M_a^T, \text{ and } K_{m,d}, L_{m,d} \text{ from Eq. (11), for any positive integer } d, \text{ Lemma 7.9 gives} \\
\text{rank}(K_{m,d}) \le \nu_m, \text{ rank}(L_{m,d}) \le \nu_m \text{ and } \text{rank}(H_{k,m,d}) \le \nu_m.
\end{gathered}
\tag{29}
\end{equation}
$$

Next, we show that the ranks of both $K_{m,\lceil v_m/m \rceil}$ and $L_{m,\lceil v_m/m \rceil}$ are $v_m$ generically.

The relation of $K_{m,d}^{(a,f)}$ and $L_{m,d}^{(a,f)}$ to the matrix $\mathcal{P}_{\ell,m,d}$. For $\ell = (\ell_1, \dots, \ell_r)$, $m$ in $\{1, \dots, n\}$, and a positive integer $d$, we define the matrix

$$
\mathcal{P}_{\ell,m,d} = \begin{pmatrix} 1 & \bar{\xi}_1 & \dots & \bar{\xi}_1^{m-1} & \bar{\lambda}_1 & \bar{\xi}_1 \bar{\lambda}_1 & \dots & \bar{\xi}_1^{m-1} \bar{\lambda}_1^{d-1} \\ & \vdots & & & \vdots & & & \vdots \\ 1 & \bar{\xi}_n & \dots & \bar{\xi}_n^{m-1} & \bar{\lambda}_r & \bar{\xi}_n \bar{\lambda}_r & \dots & \bar{\xi}_n^{m-1} \bar{\lambda}_r^{d-1} \end{pmatrix} \in \mathbb{Z}[\bar{\xi}_1, \dots, \bar{\xi}_n, \bar{\lambda}_1, \dots, \bar{\lambda}_r]^{n \times md}, \quad (30)
$$

where the rows are indexed by the variables $\bar{\xi}_1, \dots, \bar{\xi}_n$, and where each $\bar{\lambda}_i$ occurs $\ell_i$ consecutive
times, for $i = 1, \dots, r$. The following lemma summarizes the key properties of this matrix in relation
with the rank of $K_{m,d}$ and $L_{m,d}$.
LEMMA 7.10. Let $\ell, \xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r, f, a$ and $m$ be as in Proposition 7.8, and let $d$ be a positive integer. The following holds:

* the rank of $K_{m,d}$ is equal to the rank of $\mathcal{P}_{\ell,m,d}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$;

* if all $\xi_i$'s are nonzero, the rank of $L_{m,d}$ is equal to the rank of $\mathcal{P}_{\ell,m,d}(1/\xi_1, \dots, 1/\xi_n, \lambda_1, \dots, \lambda_r)$.

PROOF. We use the same notation

$$ \kappa_{m,d} : \overline{\mathbb{K}}[x, y]_{< (m,d)} \to \overline{\mathbb{K}}[x] / \langle f \rangle \quad \text{and} \quad \lambda_{m,d} : \overline{\mathbb{K}}[x]_{< n} \to \overline{\mathbb{K}}[x]_{< m}^d $$

for the mappings induced by scalar extension from $\kappa_{m,d}$ and $\lambda_{m,d}$ from Section 3.4.1.

Taking $(x^i y^j)_{0 \le i < m, 0 \le j < d}$ for basis of $\overline{\mathbb{K}}[x, y]_{< (m,d)}$ and the Lagrange basis $\mathcal{L}_1, \dots, \mathcal{L}_n$ for $\overline{\mathbb{K}}[x]_{< n}$, the matrix of $\kappa_{m,d}$ is $\mathcal{P}_{\ell,m,d}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$. This proves the first point.

To prove the second point, take $k$ in $\{1, \dots, n\}$, and let $i$ in $\{1, \dots, r\}$ be such that $a(\xi_k) = \lambda_i$. The image of the Lagrange polynomial $\mathcal{L}_k$ by $\lambda_{m,d}$ is the polynomial vector

$$ \lambda_{m,d}(\mathcal{L}_k) = ([\mathcal{L}_k]_0^{m-1}, [a\mathcal{L}_k \operatorname{rem} f]_0^{m-1}, \dots, [a^{d-1}\mathcal{L}_k \operatorname{rem} f]_0^{m-1}) \in \overline{\mathbb{K}}[x]_{<m}^d, $$

and since the Lagrange polynomials are eigenvectors of multiplication by $a$, we get

$$ \begin{aligned} \lambda_{m,d}(\mathcal{L}_k) &= ([\mathcal{L}_k]_0^{m-1}, [\lambda_i \mathcal{L}_k]_0^{m-1}, \dots, [\lambda_i^{d-1} \mathcal{L}_k]_0^{m-1}) \\ &= ([\mathcal{L}_k]_0^{m-1}, \lambda_i [\mathcal{L}_k]_0^{m-1}, \dots, \lambda_i^{d-1} [\mathcal{L}_k]_0^{m-1}). \end{aligned} $$

Let $L' \in \overline{\mathbb{K}}^{(md) \times n}$ be the matrix whose $k$-th column (for $k = 1, \dots, n$) contains the md coefficients of the entries of $\lambda_{m,d}(\mathcal{L}_k)$. This is the matrix of $\lambda_{m,d}$, if we take the Lagrange basis for the domain $\overline{\mathbb{K}}[x]_{< n}$.

Since all $\xi_i$'s are nonzero, we get $f(0) \ne 0$, so that $f$ is invertible as a power series. Because the $\mathbb{K}$-linear transformation $b \in \mathbb{K}[x]_{< m} \mapsto [b/f]_0^{m-1}$ is invertible, $L'$ has the same rank as the matrix whose columns are the coefficients of the vectors

$$ ([[\mathcal{L}_k]_0^{m-1}/f]_0^{m-1}, \lambda_i [[\mathcal{L}_k]_0^{m-1}/f]_0^{m-1}, \dots, \lambda_i^{d-1} [[\mathcal{L}_k]_0^{m-1}/f]_0^{m-1}), $$

for $i$ and $k$ as above. On the other hand, we have $[[\mathcal{L}_k]_0^{m-1}/f]_0^{m-1} = [\mathcal{L}_k/f]_0^{m-1}$ and

$$ \frac{\mathcal{L}_k}{f} = \frac{1}{f'(\xi_k)} \frac{1}{x - \xi_k}. $$

This shows that to determine the rank of $L'$, we may as well consider the vectors

$$ \left( \left[ \frac{1}{x - \xi_k} \right]_0^{m-1}, \lambda_i \left[ \frac{1}{x - \xi_k} \right]_0^{m-1}, \dots, \lambda_i^{d-1} \left[ \frac{1}{x - \xi_k} \right]_0^{m-1} \right). $$

Now, note that

$$ \left[ \frac{1}{x - \xi_k} \right]_0^{m-1} = -\xi_k \left( 1 + \frac{1}{\xi_k} x + \dots + \frac{1}{\xi_k^{m-1}} x^{m-1} \right). $$

Thus, up to the factors $-\xi_k$, taking the md coefficients of these vectors and putting them in columns gives us the transpose of $\mathcal{P}_{\ell,m,d}(1/\xi_1, \dots, 1/\xi_n, \lambda_1, \dots, \lambda_r)$. This proves the rank equality claimed in the second item. $\square$
The rank of $K_{m,d}^{(a,f)}$ and $L_{m,d}^{(a,f)}$ for $d = \lceil v_m^{(a,f)} \rceil/m$. Together with Lemma 7.10, the next lemma establishes that the generic rank of $K_{m,\lceil v_m/m \rceil}$ and $L_{m,\lceil v_m/m \rceil}$ is $v_m$. Let $\mathcal{R}_{\ell,m}$ be the $v_m \times v_m$ submatrix of $\mathcal{P}_{\ell,m,\lceil v_m/m \rceil}$ obtained by extracting the first $\min(\ell_i, m)$ rows containing $\bar{\lambda}_i$, for $i = 1, \dots, r$ (see Eq. (26)), and the first $v_m$ columns (note that $\mathcal{P}_{\ell,m,\lceil v_m/m \rceil}$ has $m \lceil v_m/m \rceil \ge v_m$ columns).

LEMMA 7.11. For $\ell = (\ell_1, \dots, \ell_r)$, $n = \ell_1 + \dots + \ell_r$ and $m$ in $\{1, \dots, n\}$, and for any pairwise distinct $\lambda_1, \dots, \lambda_r$ in $\overline{\mathbb{K}}^r$, the determinant $w_{\ell,m}(\bar{\xi}_1, \dots, \bar{\xi}_n, \lambda_1, \dots, \lambda_r)$ of the $v_m \times v_m$ matrix $\mathcal{R}_{\ell,m}$ at $\lambda_1, \dots, \lambda_r$ is nonzero.

PROOF. We prove the non-vanishing property by exhibiting a vector $(\xi_1, \dots, \xi_n) \in \overline{\mathbb{K}}^n$ for which the evaluation $w_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ is not zero. In what follows, for $i = 1, \dots, r$, recall that we write $\sigma_i = \ell_1 + \dots + \ell_{i-1}$, so that the rows involving $\bar{\lambda}_i$ in $\mathcal{P}_{\ell,m,\lceil v_m/m \rceil}$ have indices $\sigma_i + 1, \dots, \sigma_i + \ell_i$ (see Eq. (25)).

Assume first that $m$ is invertible in $\mathbb{K}$, and choose $\delta$ in $\overline{\mathbb{K}}$ such that $\delta + \lambda_i \neq 0$ for $i = 1, \dots, r$. Then, for all $i$, the polynomial $x^m - (\delta + \lambda_i)$ is separable, since its discriminant is $m^m(\delta + \lambda_i)^{m-1}$, and we choose $\xi_{\sigma_i+1}, \dots, \xi_{\sigma_i+\min(\ell_i, m)}$ to be pairwise distinct roots of this polynomial in $\overline{\mathbb{K}}$. If $m < \ell_i$, we further take $\xi_{\sigma_i+m+1}, \dots, \xi_{\sigma_i+\ell_i}$ arbitrary in $\overline{\mathbb{K}}$ (note that $w_{\ell,m}$ does not depend on these quantities). Now, for any $\xi, \lambda$ such that $\xi^m = \delta + \lambda$, and for $j \ge 1$, we have $\lambda^j = \xi^{jm} + \sum_{k=0}^{j-1} \binom{j}{k} (-\delta)^j \xi^{(j-k)m}$. Up to invertible linear combinations of its columns, $\mathcal{R}_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ is thus the Vandermonde matrix at the roots $\xi_{\sigma_i+1}, \dots, \xi_{\sigma_i+\min(\ell_i, m)}$, $i = 1, \dots, r$. Since the $\lambda_i$'s are pairwise distinct, all these roots are pairwise distinct too, so the determinant $w_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ is nonzero.

If $m$ is 0 in $\mathbb{K}$, then for all $i$, $x^m + x - \lambda_i$ is separable, since its discriminant is $(m-1)^{m-1} \neq 0$. Again, choosing distinct roots of these polynomials and performing linear combinations of the columns of $\mathcal{R}_{\ell,m}$ leads to a nonzero Vandermonde determinant. $\square$

If the rank of $\mathrm{Hk}_{m,d}^{(a,f)}$ is $v_m^{(a,f)}$ for some $d \ge 0$, then it is $v_m^{(a,f)}$ for all $d \ge \lceil v_m^{(a,f)} \rceil/m$.

LEMMA 7.12. Let $\ell, \xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r, a, f$ and $m$ be as in Proposition 7.8. If $\mathrm{Hk}_{m,d_0}$ has rank $v_m$ for some $d_0 \ge \lceil v_m/m \rceil$, and if $w_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ and $w_{\ell,m}(1/\xi_1, \dots, 1/\xi_n, \lambda_1, \dots, \lambda_r)$ from Lemma 7.11 are nonzero, then $\mathrm{Hk}_{m,d}$ has rank $v_m$ for all $d \ge \lceil v_m/m \rceil$.

PROOF. Since $w_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ is nonzero, $\mathcal{P}_{\ell,m,\lceil v_m/m \rceil}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ has rank at least $v_m$, and so does $K_{m,\lceil v_m/m \rceil}$ (Lemma 7.10).

As a result, for $d \ge \lceil v_m/m \rceil$, $K_{m,d}$ still has rank exactly $v_m$ (recall that this rank cannot exceed $v_m$, by Eq. (29)). Thus, for such $d$, there exists a nonsingular $P \in \mathbb{K}^{(md)\times(md)}$ such that $K_{m,d}P = [K_{m,\lceil v_m/m \rceil} \ 0]$, where the zero matrix is $n \times (m(d - \lceil v_m/m \rceil))$. In the same way, since $w_{\ell,m}(1/\xi_1, \dots, 1/\xi_n, \lambda_1, \dots, \lambda_r)$ is nonzero, Lemma 7.10 also implies that $L_{m,\lceil v_m/m \rceil}$ has rank $v_m$, therefore there exists a nonsingular $Q \in \mathbb{K}^{(md)\times(md)}$ such that $QL_{m,d} = [(L_{m,\lceil v_m/m \rceil})^\mathrm{T} \ 0]^\mathrm{T}$. We obtain

$$ Q \mathrm{Hk}_{m,d} P = Q \mathrm{L}_{m,d} K_{m,d} P \begin{pmatrix} \mathrm{Hk}_{m, \lceil v_m/m \rceil} & 0 \\ 0 & 0 \end{pmatrix} \in \mathbb{K}^{(md)\times(md)}, $$

which shows that for $d \ge \lceil v_m/m \rceil$ we have $\mathrm{rank}(\mathrm{Hk}_{m,d}) = \mathrm{rank}(\mathrm{Hk}_{m,\lceil v_m/m \rceil})$. $\square$

The rank of $\mathrm{Hk}_{m,r}^{(a,f)}$ is $v_m^{(a,f)}$ generically. To establish that the rank of $\mathrm{Hk}_{m,r}$ is $v_m$ for generic choices of $\xi_1, \dots, \xi_n$, we introduce a decomposition into vector spaces associated to the $\lambda_i$'s. We then study these spaces separately; their dimensions are $\min(\ell_i, m)$, respectively, leading as expected to a total dimension $\sum_{i=1}^r \min(\ell_i, m) = v_m$.