This is achieved through a description of the images of the mappings $\kappa_{m,d}$ and $\lambda_{m,d}$ in terms of polynomials. Given positive integers $\ell = (\ell_1, \dots, \ell_r)$ and $\xi_1, \dots, \xi_n$ in $\overline{\mathbb{K}}^n$, define

$$P_{i,j} = \sum_{k \in S_i} \xi_k^j \mathcal{L}_k \in \overline{\mathbb{K}}[x], \quad i = 1, \dots, r, \quad j \ge 0, \qquad (31)$$

with the Lagrange polynomials $\mathcal{L}_1, \dots, \mathcal{L}_n$ and the sets $S_1, \dots, S_r$ from Eq. (25).

**LEMMA 7.13.** Let $\ell, \xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r, a, f$ and $m$ be as in Proposition 7.8 and let $d$ be a positive integer. The image of $\kappa_{m,d}$ lies in the linear span of the $v_m$ linearly independent polynomials $P_{i,j}$ from Eq. (31), for $1 \le i \le r$ and $0 \le j < \min(\ell_i, m)$.

**PROOF.** Let $V(x, y) = \sum_{j=0}^{m-1} c_j(y)x^j$ belong to $\overline{\mathbb{K}}[x, y]_{<(m,d)}$. Lagrange interpolation gives

$$\kappa_{m,d}(V) = V(x, a) \operatorname{rem} f = \sum_{k=1}^{n} V(\xi_k, a(\xi_k)) \mathcal{L}_k.$$

Since $V(\xi_k, a(\xi_k)) = \sum_{j=0}^{m-1} c_j(a(\xi_k)) \xi_k^j$, we deduce

$$\kappa_{m,d}(V) = \sum_{i=1}^{r} \sum_{j=0}^{m-1} c_j(\lambda_i) P_{i,j}.$$

For $i = 1, \dots, r$, at most $\ell_i$ of the polynomials $P_{i,j}$, $j = 0, \dots, m-1$, can be linearly independent, since they are all linear combinations of $\ell_i$ linearly independent $\mathcal{L}_k$. On the other hand, the polynomials $P_{i,j}$ for $j = 0, \dots, \ell_i - 1$ are linearly independent, due to the linear independence of the polynomials $\mathcal{L}_k$, and the invertibility of the Vandermonde matrix $[\xi_k^j]_{0 \le j < \ell_i} \in \overline{\mathbb{K}}^{\ell_i \times \ell_i}$. This proves that the image of $\kappa_{m,d}$ is included in the span of the polynomials $P_{i,j}$, for $i = 1, \dots, r$ and $j = 0, \dots, \min(\ell_i, m) - 1$, as claimed. $\square$

This polynomial-based interpretation then allows us to use the following decomposition.

**LEMMA 7.14.** Let $\ell, \xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r, a, f$ and $m$ be as in Proposition 7.8. The rank of $H_{k,m,r}$ is the sum of the dimensions of the vector spaces

$$V_i = \operatorname{Span}([P_{i,j}]_0^{m-1}, j = 0, \dots, \min(\ell_i, m) - 1) \qquad (32)$$

with the polynomials $P_{i,j}$ from Eq. (31) for $i = 1, \dots, r$.

**PROOF.** We first claim that for $d = r, K_{m,r}$ has rank $v_m$, or equivalently (Lemma 7.10) that $\mathcal{P}_{\ell,m,r}$ has rank $v_m$ at $(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$. Indeed, we can extract from $\mathcal{P}_{\ell,m,r}$ a $v_m \times v_m$ submatrix by keeping the first $\min(\ell_i, m)$ rows indexed by $\tilde{\lambda}_i$, for $i = 1, \dots, r$, and the columns containing the monomials $\tilde{\lambda}^{j-1}, \dots, \tilde{\lambda}^{j-1} \xi^{\min(\ell_i, m)-1}$, for $j = 1, \dots, r$. That this submatrix is nonsingular at $(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ follows from an explicit factorization of its determinant [16], [21, Eq. (2.18)]; therefore $K_{m,r}$ has rank at least $v_m$. Using Eq. (29) we deduce that $K_{m,r}$ has rank exactly $v_m$ as announced, and from Lemma 7.13, we know that the image of $\kappa_{m,r}$ is the span of the polynomials $P_{i,j}$ defined in that lemma.

It follows that the rank of $H_{k,m,r}$ is the dimension of the span of the image $\lambda_{m,r}(P_{i,j})$. For $i = 1, \dots, r$, and $j = 0, \dots, \min(\ell_i, m) - 1$,

$$\lambda_{m,r}(P_{i,j}) = ([P_{i,j} \operatorname{rem} f]_0^{m-1}, [aP_{i,j} \operatorname{rem} f]_0^{m-1}, \dots, [a^{r-1}P_{i,j} \operatorname{rem} f]_0^{m-1}),$$

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

$$\lambda_{m,r}(P_{i,j}) = ([P_{i,j}]_0^{m-1}, [\lambda_i P_{i,j}]_0^{m-1}, \dots, [\lambda_i^{r-1} P_{i,j}]_0^{m-1}) \\ = ([P_{i,j}]_0^{m-1}, \lambda_i [P_{i,j}]_0^{m-1}, \dots, \lambda_i^{r-1} [P_{i,j}]_0^{m-1}).$$
The span of these vectors is unchanged if, seeing them as row vectors, one multiplies them all on the right by the inverse of the Vandermonde matrix associated to the $\lambda_i$'s. This yields a block-diagonal matrix, with blocks that span the spaces $V_i$ of the lemma. The result on the rank follows. $\square$

The dimensions of the vector spaces from Eq. (32) can now be analyzed separately.

**LEMMA 7.15.** Fix positive integers $\ell = (\ell_1, \dots, \ell_r)$ such that $\ell_1 + \cdots + \ell_r = n$, and $m$ in $\{1, \dots, n\}$. There exists a nonzero polynomial $z_{\ell,m} \in \mathbb{Z}[\tilde{\xi}_1, \dots, \tilde{\xi}_n]$ of degree at most $(n-1)(n-\nu_m)$ such that if pairwise distinct nonzero $\xi_1, \dots, \xi_n$ do not form a zero of $z_{\ell,m}$, then $V_i$ from Eq. (32) has dimension $\min(\ell_i, m)$ for all $i$.

**PROOF.** Take $i$ in $\{1, \dots, r\}$, consider the set of indices $S_i = \{\ell_1 + \cdots + \ell_{i-1} + 1, \dots, \ell_1 + \cdots + \ell_i\}$ from Eq. (25). Let then $A_i = \prod_{k \in S_i} (x - \xi_k)$, $B_i = \prod_{k \notin S_i} (x - \xi_k)$ and $C_i = 1/B_i$ mod $A_i$. Note that $A_i$ and $B_i$ have respective degrees $\ell_i$ and $n - \ell_i$, and that $C_i$ is well-defined, since $B_i$ and $A_i$ have no common root.

For $k$ in $S_i$, by construction, $B_i$ divides the Lagrange polynomial $\mathcal{L}_k$, with a quotient of degree $n - 1 - \deg(B_i) = \ell_i - 1$. In view of Eq. (31), $B_i$ divides $P_{i,j} = \sum_{k \in S_i} \xi_k^j \mathcal{L}_k$, for all $j \ge 0$, and the quotient has degree less than $\ell_i$. We now prove that it is actually equal to $x^j C_i$ rem $A_i$. Since $f = A_i B_i$, for $k$ in $S_i$, the Lagrange polynomial $\mathcal{L}_k = f/(f'(\xi_k)(x - \xi_k))$ satisfies

$$ \frac{\mathcal{L}_k}{B_i} = \frac{1}{f'(\xi_k)} \frac{f}{B_i(x - \xi_k)} = \frac{A_i}{f'(\xi_k)} \frac{x - \xi_k}{x - \xi_k} = \frac{C_i(\xi_k)}{A'_i(\xi_k)} \frac{A_i}{x - \xi_k}. $$

In particular, for $j \ge 0$, $\xi_k^j \mathcal{L}_k / B_i$ takes the value $\xi_k^j C_i(\xi_k)$ at $\xi_k$, and 0 at all other roots of $A_i$. Taking the sum over all $k$ in $S_i$ then proves our claim that $P_{i,j}/B_i = x^j C_i$ rem $A_i$. Since $\xi_1, \dots, \xi_n$ are nonzero, $B_i(0)$ as well is nonzero, so $B_i$ is invertible as a power series and $[P_{i,j}/B_i]_{0}^{m-1} = [[P_{i,j}]_{0}^{m-1}/B_i]_{0}^{m-1}$. Thus the truncated polynomials $[P_{i,j}]_{0}^{m-1}$, for $0 \le j < \min(\ell_i, m)$, are linearly independent if and only if the truncated polynomials $[x^j C_i$ rem $A_i]$_{0}^{m-1}$ are.

When $\ell_i \le m$, the polynomials $x^j C_i$ rem $A_i$ have degree less than $m$ and their linear independence follows from that of the polynomials $x^j$, $j = 0, \dots, \ell_i - 1$, since $C_i$ is invertible modulo $A_i$. Thus in this case, we always have $\dim(V_i) = \ell_i = \min(\ell_i, m)$.

When $\ell_i > m$, we are going to prove that the polynomials $x^j C_i$ rem $A_i$, $j = 0, \dots, m-1$, are linearly independent for a generic choice of $\xi_1, \dots, \xi_n$. To achieve this, define the matrix $M_{C_i}$ whose entry $(j, \ell)$ is the coefficient of $x^{\ell-1}$ in $x^{j-1} C_i$ rem $A_i$ for $j = 1, \dots, \ell_i$ and $\ell = 1, \dots, \ell_i$; this is the multiplication matrix by $C_i$ modulo $A_i$. We also consider its inverse, the multiplication matrix $M_{B_i}$ by $B_i$ modulo $A_i$.

For our claim to hold, it is enough to guarantee that the $m \times m$ leading principal minor $K_i$ of $M_{C_i}$ be nonzero. We view this minor as a rational function in $\tilde{\xi}_1, \dots, \tilde{\xi}_n$: this is done by introducing the polynomials $\tilde{A}_i = \prod_{k \in S_i} (x - \tilde{\xi}_k)$, $\tilde{B}_i = \prod_{k \notin S_i} (x - \tilde{\xi}_k)$ and $\tilde{C}_i = 1/\tilde{B}_i$ mod $\tilde{A}_i$, all of which are in $\mathbb{Q}(\tilde{\xi}_1, \dots, \tilde{\xi}_n)[x]$. We can then define the matrices $M_{\tilde{C}_i}$ and $M_{\tilde{B}_i}$ of multiplication by respectively $\tilde{C}_i$ and $\tilde{B}_i$ modulo $\tilde{A}_i$, and the $m \times m$ leading principal minor $\tilde{K}_i$ of $M_{\tilde{C}_i}$. This is a rational function of $\tilde{\xi}_1, \dots, \tilde{\xi}_n$, whose evaluation at $\tilde{\xi}_r, \dots, \tilde{\xi}_n$ gives the scalar $K_i \in \mathbb{K}$.

Note first that $\tilde{K}_i$ is not identically zero: if we evaluate all $\tilde{\xi}_g$ at 0, for $g$ in $S_i$, $\tilde{A}_i$ becomes $x^{\ell_i}$, and the matrix $M_{\tilde{C}_i}$ becomes upper triangular, with $1/B_i(0) \ne 0$. It then remains to estimate the degree of a numerator of $\tilde{K}_i$. The Schur complement formula gives $\tilde{K}_i = \det(M_{\tilde{C}_i})\tilde{L}_i$, where $\tilde{L}_i$ is the $(\ell_i - m) \times (\ell_i - m)$ lower right minor of the inverse $M_{\tilde{B}_i}$ of $M_{\tilde{C}_i}$. The determinant of $M_{\tilde{C}_i}$ is the resultant of $\tilde{C}_i$ and $\tilde{A}_i$, that is, $1/\prod_{g \in S_i, h \notin S_i} (\tilde{\xi}_g - \tilde{\xi}_h)$. On the other hand, $\tilde{L}_i$ is a polynomial in $\mathbb{Z}[\tilde{\xi}_1, \dots, \tilde{\xi}_n]$ (since $\tilde{B}_i$ and $\tilde{A}_i$ have coefficients in $\mathbb{Z}[\tilde{\xi}_1, \dots, \tilde{\xi}_n]$, and $\tilde{A}_i$ is monic in $x$).

For $s \ge 0$, write $x^s$ rem $\tilde{A}_i = c_{s,0} + \cdots + c_{s, \ell_i - 1} x^{\ell_i - 1}$, for $c_{s,t} \in \mathbb{Z}[\tilde{\xi}_1, \dots, \tilde{\xi}_n]$. By induction on $s$, we obtain the bound $\deg(c_{s,t}) \le s - t$. From this, it follows that all entries of $M_{\tilde{B}_i}$ have degree at
most $n-1$, and that $\tilde{L}_i$ has degree at most $(n-1)(\ell_i - m) \le n\ell_i$. To conclude the proof, we let $z_{\ell,m}$ be the product of the polynomials $\tilde{L}_i$, for $i$ such that $\ell_i > m$. The degree bound follows from remarking that $\sum_{\ell_i > m}(\ell_i - m) = n - v_m$. $\square$

**Genericity polynomials and degree bounds.** Until here, the conditions we have seen are the non-vanishing of $w_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$, $w_{\ell,m}(1/\xi_1, \dots, 1/\xi_n, \lambda_1, \dots, \lambda_r)$, and $z_{\ell,m}(\xi_1, \dots, \xi_n)$. When nonzero, the first two quantities allow us to apply Lemma 7.12 and obtain the rank of $Hk_{m,\lfloor n/m \rfloor}$ from any $Hk_{m,d_0}$ with $d_0 \ge \lfloor v_m/m \rfloor$; the third condition $z_{\ell,m}(\xi_1, \dots, \xi_n) \ne 0$ allows us to take $d_0 = r$ thanks to Lemmas 7.14 and 7.15

The bound on the degree of $w_{\ell,m}$ in $\tilde{\xi}_1, \dots, \tilde{\xi}_n$ follows from summing the degrees of the columns in $\mathcal{P}_{\ell,m,\lfloor v_m/m \rfloor}$. Each block of $m$ columns involves degrees $1 + \dots + (m-1) = m(m-1)/2$, and we consider $\lfloor v_m/m \rfloor$ such blocks (the last one may not be complete), for a total of at most $(v_m + m)(m-1)/2$. Next, consider the term $w_{\ell,m}(1/\tilde{\xi}_1, \dots, 1/\tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$, which is not a polynomial in the $\tilde{\xi}_i$'s. To estimate the degree of its numerator, observe that it is a $v_m \times v_m$-minor of the matrix

$$ \begin{pmatrix}
1 & \frac{1}{\tilde{\xi}_1} & \dots & \frac{1}{\tilde{\xi}_1^{m-1}} & \bar{\lambda}_1 & \frac{1}{\tilde{\xi}_1} & \dots & \frac{1}{\tilde{\xi}_1^{m-1}} & \bar{\lambda}_1^{\lfloor v_m/m \rfloor - 1} \\
& \vdots & & \vdots & & \vdots & & & \vdots \\
1 & \frac{1}{\tilde{\xi}_n} & \dots & \frac{1}{\tilde{\xi}_n^{m-1}} & \bar{\lambda}_r & \frac{1}{\tilde{\xi}_n} & \dots & \frac{1}{\tilde{\xi}_n^{m-1}} & \bar{\lambda}_r^{\lfloor v_m/m \rfloor - 1}
\end{pmatrix} $$

Factoring out (on the right) the diagonal matrix with diagonal $(1/\tilde{\xi}_i^{m-1})_{1 \le i \le n}$, we see that the non-vanishing of $w_{\ell,m}(1/\tilde{\xi}_1, \dots, 1/\tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$ is equivalent to the non-vanishing of the corresponding $v_m \times v_m$-minor $\tilde{w}_{\ell,m}$ in

$$ \begin{pmatrix}
\tilde{\xi}_1^{m-1} & \tilde{\xi}_1^{m-2} & \dots & 1 & \tilde{\xi}_1^{m-1} \bar{\lambda}_1 & \tilde{\xi}_1^{m-2} \bar{\lambda}_1 & \dots & \bar{\lambda}_1^{\lfloor v_m/m \rfloor - 1} \\
& \vdots & & \vdots & \vdots & \vdots & & \vdots \\
\tilde{\xi}_n^{m-1} & \tilde{\xi}_n^{m-2} & \dots & 1 & \tilde{\xi}_n^{m-1} \bar{\lambda}_r & \tilde{\xi}_n^{m-2} \bar{\lambda}_r & \dots & \bar{\lambda}_r^{\lfloor v_m/m \rfloor - 1}
\end{pmatrix} $$

The degree upper bound for $\tilde{w}_{\ell,m}$ is $(v_m + m)(m-1)/2$, as for $w_{\ell,m}$.

We then take $\Gamma_{\ell,m} = w_{\ell,m}\tilde{w}_{\ell,m}z_{\ell,m}$ to prove Proposition 7.8. For the degree estimate, note that $(v_m + m)(m-1) + (n-1)(n - v_m) \le 2n^2$. For correctness, take pairwise distinct nonzero $\tilde{\xi}_1, \dots, \tilde{\xi}_n$ in $\mathbb{K}$ and let $a \in \mathbb{K}[x]$ take distinct values $\tilde{\lambda}_1, \dots, \tilde{\lambda}_r$ at $\tilde{\xi}_1, \dots, \tilde{\xi}_n$, with multiplicities $\ell_1, \dots, \ell_r$. As before, we write $f = (x - \tilde{\xi}_1) \cdots (x - \tilde{\xi}_n)$, and we assume that $f$ is in $\mathbb{K}[x]$. Finally, we suppose that $\Gamma_{\ell,m}(\tilde{\xi}_1, \dots, \tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$ is nonzero. Lemmas 7.14 and 7.15 show that for $d_0 = r$, we have $\text{rank}(Hk_{m,r}^{(a,f)}) = v_m$. Since $r \ge \lfloor v_m/m \rfloor$, by Lemma 7.12, it is then also the case for $Hk_{m,d}^{(a,f)}$ for all $d \ge \lfloor v_m/m \rfloor$, as claimed.

The only remaining claim is that for any pairwise distinct $\tilde{\lambda}_1, \dots, \tilde{\lambda}_r$, $\Gamma_{\ell,m}(\tilde{\xi}_1, \dots, \tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$ is a nonzero polynomial in $\tilde{\xi}_1, \dots, \tilde{\xi}_n$. That $z_{\ell,m}$ is nonzero is in Lemma 7.15 (this polynomial does not depend on $\tilde{\lambda}_1, \dots, \tilde{\lambda}_r$); Lemma 7.11 proves that $w_{\ell,m}(\tilde{\xi}_1, \dots, \tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$ is nonzero. That lemma also implies that $w_{\ell,m}(1/\tilde{\xi}_1, \dots, 1/\tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$ is nonzero (as a rational function), and as a consequence, this is also the case for $\tilde{w}_{\ell,m}(\tilde{\xi}_1, \dots, \tilde{\xi}_n, \tilde{\lambda}_1, \dots, \tilde{\lambda}_r)$. The claim for $\Gamma_{\ell,m}$ is thus proved.

# 8 A RANDOMIZED COMPOSITION ALGORITHM THROUGH CHANGE OF BASIS

In this section we give the base case of our modular composition algorithm that is used when $f$ is either separable or purely inseparable (which includes the case of power series). The core Algorithm `MODULARCOMPOSITIONBASECASE` is studied in Section 8.1, and a variation for computing annihilating polynomials is given in Section 8.2.

The algorithm of Section 4.2 performs bivariate modular composition within our target complexity bound, assuming the knowledge of a matrix of relations with appropriate dimension and degree.
Since such a matrix of relations of $M_m^{(\alpha,f)}$ may not exist for general $a$ and $f$, Algorithm `MODULARCOMPOSITIONBASECASE` transports the computation of $g(a)$ in $\mathbb{A} = \mathbb{K}[x]/\langle f \rangle$ to an isomorphic algebra which is expected to be more favorable to the computation.

More precisely, we pick a random $\gamma \in \mathbb{K}[x]_{<n}$; generically, its minimal polynomial $\mu_\gamma \in \mathbb{K}[y]$ has degree $n$ and is also its characteristic polynomial $\chi_\gamma$, so that the powers of $\gamma$ generate $\mathbb{A}$. This induces the $\mathbb{K}$-algebra isomorphism $\phi_\gamma$ of Eq. (19); Step 5 of Algorithm `MODULARCOMPOSITIONBASECASE` then computes a polynomial representative $\alpha$ of $\phi_\gamma(a \mod f)$ using the change of basis algorithm of Section 6. Note that a matrix of relations $R_m^{(\gamma,f)}$ is also obtained at Step 5 in preparation for the final stage. Then, with good probability, the conditions for the efficient computation of a certified matrix of relations $R_m^{(\alpha,\mu_\gamma)}$ of $M_m^{(\alpha,\mu_\gamma)}$ via the approach of Section 5.4 are fulfilled. Step 8 of Algorithm `MODULARCOMPOSITIONBASECASE` computes this matrix of relations, which then allows us to obtain the polynomial $\beta = g(\alpha)$ rem $\mu_\gamma$ at Step 9 as seen in Section 4.2. The solution $b = g(a)$ rem $f$ to the initial problem is finally recovered by applying $\phi_\gamma^{-1}$ to $\beta \mod \mu_\gamma$, which amounts to computing $b = \beta(\gamma)$ rem $f$. Since we already have $R_m^{(\gamma,f)}$ at our disposal, $b$ is obtained with the algorithm of Section 4.2 as well.

**Proposition 8.1** in Section 8.1 shows the correctness of this strategy and bounds its complexity. We then study the probability of success for $f$ separable and $f$ purely inseparable. The main point is to ensure that appropriate matrices of relations $R_m^{(\gamma,f)}$ and $R_m^{(\alpha,\mu_\gamma)}$ are actually available. For Steps 5 and 11 where a random $\gamma$ is involved, we directly rely on the generic properties of the associated block Hankel matrix $H_k^{(\gamma,f)}_{m,\lfloor n/m \rfloor}$ (Proposition 7.6). For the computation of $R_m^{(\alpha,\mu_\gamma)}$ we use the fact that $\alpha$ and $\mu_\gamma$ are sufficiently generic, hence also give access to good properties for the associated block Hankel matrix after the change of basis.

The probability of failure for a general separable $f$ is bounded in Section 8.3. The power series case and, more generally, the case of purely inseparable $f$ are treated in Sections 8.4 and 8.5. For such $f$, the success of Algorithm `MODULARCOMPOSITIONBASECASE` is proven in Section 8.4 under some assumptions on the valuation of the input polynomial $a$ and the characteristic of $\mathbb{K}$. Still in the case of $f$ purely inseparable, a complete algorithm is then given in Section 8.5: when the valuation is large (with respect to the target value $m \sim n^\eta$ with $\eta$ from Eq. (3)), then the minimal polynomial of $a$ modulo $f$ has small degree and we use the extension of Shoup's algorithm seen in Section 3.1.3. For fields $\mathbb{K}$ of small characteristic, we adapt Bernstein's composition algorithm for power series [6] to our general context.

## 8.1 Randomized composition

The procedure is detailed in Algorithm `MODULARCOMPOSITIONBASECASE`. It uses $n+m$ parameters from $\mathbb{K}$ that are available as a sequence $r$ of length $n+m$. The coefficients of the random polynomial $\gamma$ are given as part of the input as $r_i$, for $3 \le i \le n+2$; we require further parameters in order to reduce to the case where $f(0) \ne 0$ and $\gcd(a, f) = 1$ (Remarks 3.8 and 5.7), and for the random column combination performed by Algorithm `MATRIXOFRELATIONS`.

The parameter $m$ could be taken arbitrarily in $\{1, \dots, n\}$, but we choose the specific value $m = \lfloor n^\eta \rfloor$, with $\eta$ from Eq. (3), as this choice minimizes the overall cost. The following proposition describes the output of the procedure; the probability of failure is bounded in Sections 8.3 and 8.4.

**PROPOSITION 8.1.** Given $f \in \mathbb{K}[x]$ of degree $n$, $a \in \mathbb{K}[x]_{<n}$, $g \in \mathbb{K}[y]$ with $\deg(g) = O(n)$ and $r \in \mathbb{K}^{n+m}$ with $m = \lfloor n^\eta \rfloor$ and $\eta$ from Eq. (3), Algorithm `MODULARCOMPOSITIONBASECASE` returns either $g(a)$ rem $f$ or `FAIL`; it uses $\tilde{O}(n^\kappa)$ operations in $\mathbb{K}$, with $\kappa < 1.43$ as in Eq. (1).

**PROOF.** If $n=1$ then as $a$ has degree 0, the result is $g(a) \in \mathbb{K}$ and the algorithm is correct. The rest of the proof assumes $n > 1$.
<caption>**Algorithm 8.1 MODULARCOMPOSITIONBASECASE(f, a, g, r)**</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a \in \mathbb{K}[x]_{<n}$, $g \in \mathbb{K}[y]$, $r \in \mathbb{K}^{n+\lfloor n' \rfloor}$
Output: $b = g(a)$ rem $f$ or FAIL
1: **if** $n = 1$ **then** return $g(a)$
2: $g \leftarrow g(y - r_1)$, $a \leftarrow a(x) + r_1$; **if** $\gcd(a, f) \neq 1$ **then return** FAIL
3: $f \leftarrow f(x + r_2)$; $a \leftarrow a(x + r_2)$; **if** $f(0) = 0$ **then return** FAIL
4: $m \leftarrow \lfloor n' \rfloor$
5: ▶ **Change of basis:** compute a polynomial $\alpha$ such that $\alpha \equiv \phi_Y(a \bmod f) \bmod \mu_Y$
6: ▶ **Getting a basis of relations** $R^{(Y,f)}$ and the minimal polynomial $\mu_Y$ of $\gamma$ $\bmod f$
7: $\gamma \leftarrow r_3 + r_4x + \dots + r_{n+2}x^{n-1}$
8: $(R^{(Y,f)}, \mu_Y, \alpha) \leftarrow \text{CHANGE OF BASIS}(f, \gamma, a, m, \lfloor n/m \rfloor)$
9: ▶ **Algorithm 6.1**
10: **if** this call returned FAIL **then return** FAIL
11: **if** $\mu_Y(0) = 0$ **then return** FAIL
12: **if** $\mu_Y(0) \neq 0$ **then return** FAIL
13: substitute “$y$” by “$x$” in $\mu_Y$ and $\alpha$, which are then in $\mathbb{K}[x]$
14: ▶ **Compute a matrix of relations for** $(\alpha, \mu_Y)$
15: $R^{(\alpha, \mu_Y)} \leftarrow \text{MATRIX OF RELATIONS}(\mu_Y, \alpha, m, \lfloor n/m \rfloor, (r_{n+i})_{3 \le i \le m})$
16: ▶ **Algorithm 5.2**
17: **if** this call returned FAIL **then return** FAIL
18: ▶ **Bivariate modular composition in the new basis:** $\beta \equiv g(\alpha) \bmod \mu_Y$
19: ▶ **Algorithm 4.1**
20: $\beta \leftarrow \text{BIVARIATE MODULAR COMPOSITION WITH RELATION MATRIX}(\mu_Y, \alpha, g, R^{(\alpha, \mu_Y)})$
21: substitute “$x$” by “$y$” in $\beta$, which is then in $\mathbb{K}[y]$
22: ▶ **Inverse change of basis:** $b \equiv \phi^{-1}(\beta \bmod \mu_Y) \bmod f$
23: $b \leftarrow \text{BIVARIATE MODULAR COMPOSITION WITH RELATION MATRIX}(f, \gamma, \beta, R^{(Y,f)})$ ▶ **Algorithm 4.1**
24: **return** $b(x - r_2)$

The first two steps ensure that $\gcd(a, f) = 1$ and $f(0) \neq 0$. This does not impact the complexity, as shifting a polynomial of degree $O(n)$ can be achieved in $\tilde{O}(n)$ arithmetic operations [7, Chap. 1, Pb. 3.5]. The same observation applies to the last step.

At Step 5, if Algorithm CHANGEOFBASIS does not return FAIL then by Proposition 6.1 the matrix $R^{(Y,f)}$ is a basis of relations of $\mathcal{M}_m^{(Y,f)}$, $\mu_Y = \chi_Y$, and $\alpha(\gamma) \equiv a \bmod f$. It follows that $\mu_\alpha = \mu_a$ since the quotient algebras are isomorphic, and $\mu_\alpha(0) = \mu_a(0)$ implies $\gcd(\alpha, \mu_Y) = \gcd(a, f) = 1$. If the test at Step 6 does not fail then the specifications for the call to Algorithm MATRIXOFRELATIONS are met; from Proposition 5.8, if Step 8 does not return FAIL then the matrix $R^{(\alpha, \mu_Y)}$ is a matrix of relations in $(\mu_Y, y - \alpha)$. Both these matrices of relations have dimension at most $2(m-1)$, and degree at most $2\lfloor n/m \rfloor$; they are obtained in $\tilde{O}(m^\omega d + c(n, m, d)) = \tilde{O}(m^{\omega_2/2})$ operations, with $d = \lfloor n/m \rfloor$. This is $\tilde{O}(n^\kappa)$ arithmetic operations, according to Eqs. (1) and (3) and the choice of $m$ at Step 4.

The variable substitutions at Steps 7 and 10 are harmless; they make notation match with that in Algorithms MATRIXOFRELATIONS and BIVARIATEMODULARCOMPOSITIONWITHRELATIONMATRIX.

At Step 9, within the same complexity bound as above by Proposition 4.4, $\beta$ is computed such that $\beta \equiv g(\alpha) \bmod \mu_Y$ (these polynomials are temporarily in $x$). After the substitution of Step 10 the latter relation implies the existence of a polynomial $h \in \mathbb{K}[y]$ such that

$$ \beta(y) = g(\alpha(y)) + h(y)\mu_Y(y). $$

Since $\mu_Y(\gamma) \equiv 0 \bmod f$, evaluating this identity at $y = \gamma$ results in $b = \beta(\gamma) = g(a)$ rem $f$ at Step 11. □
## 8.2 Randomized annihilating polynomial

If the choice of $\gamma$ ensures that the isomorphism $\phi_{\gamma}$ is well defined (the powers of $\gamma$ generate $\mathbb{A}$), then a univariate polynomial $\mu$ over $\mathbb{K}$ is such that $\mu(a) \equiv 0 \bmod f$ if and only $\mu(\alpha) \equiv 0 \bmod \mu_{\gamma}$. Since Algorithm MODULARCOMPOSITIONBASECASE computes a matrix of relations in $\langle\mu_{\gamma}, y - \alpha\rangle$ at Step 8, an algorithm for computing such a $\mu$ follows from the results of Section 4.3.

<caption>**Algorithm 8.2 ANNIHILATINGPOLYNOMIAL(f, a, r)**</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a \in \mathbb{K}[x]_{<n}$, $r \in \mathbb{K}^{n+\lfloor n\eta \rfloor}$
Output: $\mu$ nonzero in $\mathbb{K}[y]_{\le 4n}$ such that $\mu(a) \equiv 0 \bmod f$ or FAIL
1: **if** $n = 1$ **then return** $y - a$
2: ▶ Compute a matrix of relations $R^{(\alpha, \mu_{\gamma})}$ for $(\alpha, \mu_{\gamma})$ with $\alpha = \phi_{\gamma}(a)$ execute Steps 3 to 8 of Algorithm 8.1
   **if** FAIL has been returned by one of these steps **then return** FAIL
3: $\mu \leftarrow \det(R^{(\alpha, \mu_{\gamma})})$
4: **return** $\mu$

$\triangleright a \in \mathbb{K}$

**COROLLARY 8.2.** Given $f \in \mathbb{K}[x]$ of degree $n$, $a \in \mathbb{K}[x]_{<n}$ and $r \in \mathbb{K}^{n+m}$ with $m = \lfloor n\eta \rfloor$ and $\eta$ from Eq. (3), Algorithm ANNIHILATINGPOLYNOMIAL returns either FAIL or a nonzero $\mu \in \mathbb{K}[y]_{\le 4n}$ such that $\mu(a) \equiv 0 \bmod f$; it uses $\tilde{O}(n^\kappa)$ operations in $\mathbb{K}$, with $\kappa < 1.43$ as in Eq. (1).

**PROOF.** If $n=1$ then as $a \in \mathbb{K}$, $\mu = y - a$ is such that $\mu(a) = 0$ and the algorithm is correct. Now assume that $n > 1$. The annihilating polynomials are left unchanged by the substitution $x \leftarrow x + r_2$. As in the proof of Proposition 8.1, if failure does not occur then Step 8 computes a matrix of relations of $M_{m'}^{(\alpha, \mu_{\gamma})}$, for some $m' \le 2(m-1)$, within the claimed complexity bound; this matrix has degree at most $2\lfloor n/m \rfloor$. Then Proposition 4.5 shows that $\mu$ annihilates $\alpha \bmod \mu_{\gamma}$ and thus $a \bmod f$, and that it has degree $\deg(\mu) \le 4(m-1)\lfloor n/m \rfloor$. This is at most $4n$ when $m \le \sqrt{n}$, which is the case when $m = \lfloor n\eta \rfloor$ with $\eta$ as in Eq. (3). The complexity then follows from the proof of Proposition 8.1 and Proposition 4.5 again. $\square$

## 8.3 Success of randomization for separable $f$

The probabilistic properties of the previous algorithms in the separable case are summarized in the following.

**PROPOSITION 8.3.** Let $a, f$ be polynomials in $\mathbb{K}[x]$ and $g$ be in $\mathbb{K}[y]$, with $f$ separable of degree $n$ and $\deg(a) < n$. If $r_1, \dots, r_{n+\lfloor n\eta \rfloor} \in \mathbb{K}$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then Algorithms MODULARCOMPOSITIONBASECASE and ANNIHILATINGPOLYNOMIAL return FAIL with probability at most $6n^2/\text{card}(S)$.

**PROOF.** The success of modular composition in Algorithm MODULARCOMPOSITIONBASECASE and of the computation of an annihilating polynomial in Algorithm ANNIHILATINGPOLYNOMIAL relies on: finding good shifts $r_1$ and $r_2$ in the first two steps; a choice of $\gamma$ such that $\mu_{\gamma}(0) \ne 0$ and $\mu_{\gamma}$ has degree $n$; the availability of matrices of relations $R^{(\gamma, f)}$ and $R^{(\alpha, \mu_{\gamma})}$. The probability estimate is obtained by showing the existence of polynomials whose zero sets contain the values of the parameters $r_i$ where these properties do not hold. The probability of avoiding these zero sets is then handled by the Schwartz-Zippel lemma. In what follows, as in the algorithm, we write $m = \lfloor n\eta \rfloor$.

(1) A value of $r_1$ such that $\gcd(a + r_1, f) \ne 1$. The resultant of $a(x) + r_1$ and $f(x)$ is nonzero of degree $n$ in $r_1$. Bad choices thus occur with probability at most $n/\text{card}(S)$.
(2) A value of $r_2$ such that $f(0) \neq 0$ after the shift "$x \leftarrow x + r_2$". The same reasoning as above applies to the coefficient of degree zero of $f(x + r_2)$.

The next properties all concern the same parameters ($r_3, \dots, r_{n+m}$), so their failures are not independent events, and their joint probability is bounded using a product of polynomials encoding each of them. Below, we write $\bar{\gamma} = \bar{\gamma}_0 + \dots + \bar{\gamma}_{n-1}x^{n-1}$, with the $\bar{\gamma}_i$'s new indeterminates, and consider polynomials in $\overline{\mathbb{K}}[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$ to quantify probabilities of failure.

(3) The constant coefficient $\mu_\gamma(0)$ is not 0. Write $f = c(x - \varphi_1) \cdots (x - \varphi_n)$, for pairwise distinct $\varphi_i$ in $\overline{\mathbb{K}}$ and $c \in \mathbb{K} \setminus \{0\}$. The roots of $\mu_\gamma$ are the values $\bar{\gamma}(\varphi_i)$, so $\mu_\gamma(0)$ being nonzero is equivalent to $\gcd(\gamma, f)$ being trivial. Thus, we let $\Delta_0 \in \mathbb{K}[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$ be the resultant of $\bar{\gamma}$ and $f$. This polynomial has degree $n$, and choosing $\gamma = 1$ shows that it is not identically zero.

(4) The minimal polynomial $\mu_\gamma$ has degree $n$. For any $\gamma = \gamma_0 + \dots + \gamma_{n-1}x^{n-1}$ in $\mathbb{K}[x]_{<n}$, the characteristic polynomial $\chi_\gamma \in \mathbb{K}[y]$ of $\gamma$ mod $f$ factors over $\overline{\mathbb{K}}[y]$ as $\chi_\gamma = \prod_{i=1}^n (y - \xi_i)$, where $\xi_i = \gamma(\varphi_i)$ for all $i$. We can thus let $\Delta_1 \in \overline{\mathbb{K}}[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$ be the product $\prod_{1 \le i < j \le n} (\bar{\gamma}(\varphi_i) - \bar{\gamma}(\varphi_j))$. This is a polynomial of degree $n(n-1)/2$, and the previous discussion shows that $\Delta_1(\gamma_0, \dots, \gamma_{n-1}) \neq 0$ implies that $\chi_\gamma$ is separable. In that case, $\chi_\gamma = \mu_\gamma$ by degree considerations. Finally, the polynomial $\Delta_1$ itself is nonzero since its value at $(0, 1, 0, \dots, 0)$, i.e. at $\gamma = x$, is not zero.

(5) The computation of $R^{(\gamma, f)}$ does not fail. Since $f(0) \neq 0$, Proposition 7.6 shows that the associated block Hankel matrix $\mathrm{Hk}_{m, [n/m]}^{(\gamma, f)}$ has rank $n$ as soon as the coefficients of $\gamma$ avoid the zero set of a polynomial $\Delta_{f,m}$ of degree at most $2n^2/m$.

When this condition holds, Proposition 7.1 shows that the matrix fraction $X^T(yI_n - M_\gamma)^{-1}X$ is describable in degree $\lceil n/m \rceil$, and that $v_m^{(\gamma, f)} = v_{m,m}^{(\gamma, f)}$. Since we also have $\gcd(\gamma, f) = 1$ by the item above, and since the minimal polynomial $\mu_\gamma$ of $\gamma$ mod $f$ has degree $n$, Proposition 6.1 concludes that the computation of $R^{(\gamma, f)}$ is successful.

(6) The rank of $\mathrm{Hk}_{m,d}^{(\alpha, \mu_\gamma)}$ is equal to $v_m^{(\alpha, \mu_\gamma)}$ for $d \ge \lceil v_m^{(\alpha, \mu_\gamma)} \rceil/m$. When the previous properties are all satisfied, there exists a $\mathbb{K}$-algebra isomorphism $\phi_\gamma : \overline{\mathbb{K}}[x]/\langle f \rangle \to \overline{\mathbb{K}}[y]/\langle \mu_\gamma \rangle$ which maps $a$ to $\alpha$ such that $\alpha(\gamma) \equiv a \pmod f$. Up to changing the indices of the roots $\varphi_i$, we can assume that $a$ takes values $\lambda_1, \dots, \lambda_r$ at $\varphi_1, \dots, \varphi_n$ with multiplicities $\ell_1, \dots, \ell_r$, for some positive integers $\ell_1, \dots, \ell_r$, and pairwise distinct $\lambda_1, \dots, \lambda_r$ in $\overline{\mathbb{K}}$ (as in Section 7.4, the $\varphi_i$'s are assumed to be ordered such that $a(\varphi_1) = \dots = a(\varphi_{\ell_1}) = \lambda_1$, etc). Then, since $\xi_i = \gamma(\varphi_i)$ for all $i$, the relation $\alpha(\gamma) \equiv a \pmod f$ implies that $\alpha(\xi_i) = a(\varphi_i)$ for all $i$, so that $\alpha$ takes the values $\lambda_1, \dots, \lambda_r$ at $\xi_1, \dots, \xi_n$ with multiplicities $\ell_1, \dots, \ell_r$.

The assumptions of Proposition 7.8 are satisfied. If $\Gamma_{\ell,m} \in \mathbb{Z}[\bar{\xi}_1, \dots, \bar{\xi}_n, \bar{\lambda}_1, \dots, \bar{\lambda}_r]$ is the polynomial defined in that proposition, then when $\Gamma_{\ell,m}(\xi_1, \dots, \xi_n, \lambda_1, \dots, \lambda_r)$ is nonzero, the rank of $\mathrm{Hk}_{m,d}^{(\alpha, \mu_\gamma)}$ is $v_m^{(\alpha, \mu_\gamma)}$ for $d \ge \lceil v_m^{(\alpha, \mu_\gamma)} \rceil/m$.

The relevant polynomial is thus $\Delta_3 = \Gamma_{\ell,m}(\bar{\gamma}(\varphi_1), \dots, \bar{\gamma}(\varphi_n), \lambda_1, \dots, \lambda_r) \in \overline{\mathbb{K}}[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$. Proposition 7.8 states that $\Gamma_{\ell,m}(\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}, \lambda_1, \dots, \lambda_r)$ is nonzero of degree at most $2n^2$; this is thus also the case for $\Delta_3$, since the transformation $(\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}) \mapsto (\bar{\gamma}(\varphi_1), \dots, \bar{\gamma}(\varphi_n))$ is linear and invertible (its matrix is the Vandermonde matrix at $\varphi_1, \dots, \varphi_n$).

(7) The computation of $R^{(\alpha, \mu_\gamma)}$ does not fail. When the previous properties are all satisfied, Proposition 7.1 applies with $a = \alpha$ and $f = \mu_\gamma$ and shows that $v_{m,m}^{(\alpha, \mu_\gamma)} = v_m^{(\alpha, \mu_\gamma)}$ and $X^T(yI_n - M_\alpha)^{-1}X$ is describable in degree $\lceil v_m^{(\alpha, \mu_\gamma)} \rceil/m$, where $M_\alpha$ is the multiplication matrix of $\alpha$ modulo $\mu_\gamma$. Since $\mu_\gamma(0) \neq 0$ and $\gcd(\alpha, \mu_\gamma) = \gcd(a, f) = 1$, the assumptions of Proposition 5.8 are satisfied for the
successful computation of $R^{(\alpha, \mu_Y)}$ (Algorithm `MATRIXOFRELATIONS`) with a probability of failure depending on the choices of $(r_{n+3}, \dots, r_{n+m})$ and bounded by $(m-1)/\text{card}(S)$.

Case $n=1$. In that situation steps, (5)–(7) above simplify. Since its top left corner is the identity matrix, the rank of the block-Hankel matrix is at least 1, which is equal to $n$, and thus (5)–(7) succeed with probability 1 in that case.

*Probability bounds.* The polynomial $\Delta_0\Delta_1\Delta_{f,m}\Delta_3 \in \overline{\mathbb{K}}[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$ is nonzero and has degree at most

$$d_{n,m} = \frac{n(n-1)}{2} + n + \frac{2n^2}{m} + 2n^2.$$

A choice of $(r_3, \dots, r_{n+2})$ that avoids its zero set ensures that the properties (3)–(6) hold. The other probabilities have been discussed in steps (1), (2) and (7) above. In summary, the probability of success is at least

$$\begin{cases} \left(1 - \frac{1}{\text{card}(S)}\right)^3 \ge 1 - \frac{3}{\text{card}(S)} & \text{if } n=1, \\ \left(1 - \frac{n}{\text{card}(S)}\right)^2 \left(1 - \frac{d_{n,m}}{\text{card}(S)}\right) \left(1 - \frac{m-1}{\text{card}(S)}\right) \ge 1 - \frac{2n+d_{n,m}+m-1}{\text{card}(S)} & \text{otherwise.} \end{cases}$$

In the second expression, dividing the numerator of the fraction for $n \ge 2$ by $n^2$ gives

$$\frac{5}{2} + \frac{5}{2n} + \frac{2}{m} + \frac{m-1}{n^2},$$

which decreases as a function of $n$ for $n \ge 0$ and, for fixed $n$, decreases as a function of $m$ for $m \le n$. Thus it reaches its maximum at $m=1$, $n=2$, where its value is $23/4 < 6$, proving the probability bound for $n \ge 2$, while $3 < 6$ deals with the case $n=1$.

The assertion for Algorithm `ANNIHILATINGPOLYNOMIAL` follows: Step 2 apart, it fails in the same cases as Algorithm `MODULARCOMPOSITIONBASECASE`. $\square$

## 8.4 Success of randomization for $f$ purely inseparable: small valuation

*Definition 8.4.* A degree $n$ polynomial $f$ in $\mathbb{K}[x]$ is purely inseparable if it has only one root in an algebraic closure $\overline{\mathbb{K}}$, so that it factors as $f = (x - \xi)^n$ in $\overline{\mathbb{K}}[x]$; if $n$ is a unit in $\mathbb{K}$, $\xi$ itself is in $\mathbb{K}$.

In this section, we study the probabilistic aspects of Algorithm `MODULARCOMPOSITIONBASECASE` for such polynomials, under two additional assumptions on the valuation $v = \text{val}_\xi(a - a(\xi))$: it is not 0 in $\mathbb{K}$, and it is at most the value chosen for $m$ (which is $\lceil n^\eta \rceil$ in the algorithm, for the target complexity bound). The other cases are discussed in the next section.

**PROPOSITION 8.5.** Let $a, f$ be polynomials in $\mathbb{K}[x]$ and $g$ be in $\mathbb{K}[y]$, with $f = (x - \xi)^n \in \mathbb{K}[x]$ where $\xi \in \overline{\mathbb{K}}$, and $\text{deg}(a) < n$. Let $p$ be the characteristic of $\mathbb{K}$. Suppose that $v = \text{val}_\xi(a - a(\xi))$ satisfies the following inequalities, with $\eta$ as in Eq. (3):

$$v \le \lceil n^\eta \rceil, \quad p = 0 \text{ or } v < p.$$

Take $r_1 = 0$ if $\text{gcd}(a, f) = 1$ and $r_1 = 1$ otherwise, $r_2 = 0$ if $f \ne x^n$ and $r_2 = 1$ otherwise. If $r_3, \dots, r_{n+\lceil n^\eta \rceil}$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then Algorithms `MODULARCOMPOSITIONBASECASE` and `ANNIHILATINGPOLYNOMIAL` return `FAIL` with probability at most $2n^4/\text{card}(S)$

*PROOF.* The proof follows the same steps as in Section 8.3. As before, we write $m = \lceil n^\eta \rceil$.

(1), (2) Values of $r_1$ and $r_2$. The choice of $r_1$ gives $\text{gcd}(a + r_1, f) = 1$, and $r_2$ modifies the constant coefficient of $f$ if necessary. The first two steps of Algorithm `MODULARCOMPOSITIONBASECASE` therefore provide polynomials that satisfy $\text{gcd}(a, f) = 1$, $f(0) \ne 0$, and $v = \text{val}_\xi(a - a(\xi)) \le m$.
(3) The constant coefficient $\mu_Y(0)$ is not 0. For any $\gamma = \gamma_0 + \cdots + \gamma_{n-1}x^{n-1}$, the roots of the characteristic polynomial $\chi_\gamma$ of $\gamma$ modulo $f$ are the values taken by $\gamma$ at the roots of $f$, counted with multiplicities. Since $f = (x - \xi)^n$ over $\mathbb{K}[x]$, this implies that $\chi_\gamma = \prod(y - \gamma(\xi))^n$. The minimal polynomial $\mu_\gamma$ then admits a similar factorization as $\prod(y - \gamma(\xi))^c$, for some positive $c$.

Set $\Delta_0(\bar{\gamma}_0, \dots, \bar{\gamma}_n) = \sum_{i=0}^{n-1} \bar{\gamma}_i \xi^i$; this is a (nonzero) polynomial of degree 1 which is such that $\Delta_0(\gamma_0, \dots, \gamma_{n-1}) = \gamma(\xi)$, so the non-vanishing of this quantity gives the same property for $\mu_\gamma(0)$.

(4) The minimal polynomial $\mu_\gamma$ has degree $n$. Consider now $\Delta_1(\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}) = \sum_{i=1}^{n-1} i \bar{\gamma}_i \xi^{i-1}$, which is also a nonzero polynomial of degree 1. It is such that $\Delta_1(\gamma_0, \dots, \gamma_{n-1}) = \gamma'(\xi)$, so the non-vanishing of this quantity implies that $\text{val}_\xi(\gamma - \gamma(\xi)) = 1$. This implies that the powers $1, \gamma - \gamma(\xi), (\gamma - \gamma(\xi))^2, \dots, (\gamma - \gamma(\xi))^{n-1}$ rem $(x - \xi)^n$ have respective valuations $0, 1, \dots, n-1$ at $\xi$, and thus are linearly independent. It follows that the minimal polynomial of $\gamma - \gamma(\xi)$ has degree $n$, and the same then holds for $\gamma$ itself.

(5) The computation of $R^{(\gamma, f)}$ does not fail. Here the argument of the previous section applies verbatim and relies on a polynomial $\Delta_{f,m}$ of degree at most $2n^2/m$.

(6) The rank of $\text{Hk}_{m,d}^{(\alpha, \mu_\gamma)}$ is equal to $n$ for $d \ge \lfloor n/m \rfloor$. This step is the difficult one in the proof; note that the statement slightly deviates from the one in the separable case in the definition of the threshold degree $\lfloor n/m \rfloor$.

The result is obtained by bounding the degree of the numerator of a nonzero $n \times n$ minor of $\text{Hk}_{m, \lfloor n/m \rfloor}^{(\alpha, \mu_\gamma)}$, seen as a polynomial in $\gamma_0, \dots, \gamma_{n-1}$. We first show the existence of $\gamma \in \mathbb{K}[x]_{<n}$ and $\alpha \in \mathbb{K}[y]_{<n}$ such that the block Hankel matrix $\text{Hk}_{m, \lfloor n/m \rfloor}^{(\alpha, \mu_\gamma)}$ has rank $n$. This implies the existence of a nonzero $n \times n$ minor of this matrix; the degree of this minor as a polynomial in the coefficients of $\alpha$ and $\mu_\gamma$ is controlled by Lemma 7.7. These in turn are related to the coefficients of $\gamma$, using its explicit form for $\mu_\gamma$ and a linear system for the coefficients of $\alpha$.

(6a) Generic behaviour. We start by proving the existence of $\alpha$ of degree $m$ in $\overline{\mathbb{K}}[y]$ and $\gamma$ in $\overline{\mathbb{K}}[x]_{<n}$ such that we have $\alpha(\gamma) \equiv a \bmod f$, $\gamma(\xi) \ne 0$ and $\gamma'(\xi) \ne 0$.

Write $a = a_0 + a_\nu(x - \xi)^\nu + \cdots + a_{n-1}(x - \xi)^{n-1}$, with $a_0 = a(\xi)$ and, by definition of $\nu = \text{val}_\xi(a - a(\xi))$, $a_\nu \ne 0$ and $\nu > 0$. Since we also assume that the characteristic $p$ of $\mathbb{K}$ is either zero, or greater than $\nu$, this means in particular that $\nu$ is a unit in $\mathbb{K}$. Let

$$ \tilde{a}(x) = \frac{a - a_0}{a_\nu(x - \xi)^\nu} = 1 + \sum_{1 \le i < n - \nu} \tilde{a}_i (x - \xi)^i, $$

with coefficients $\tilde{a}_i = a_{i+\nu}/a_\nu$.

- If $\nu = m$, we define $\tilde{\alpha}(y) = y^\nu = y^m$. Since $\nu \ne 0$ in $\mathbb{K}$, $\tilde{\alpha}'(1) \ne 0$ and Newton iteration guarantees the existence of a unique $\tilde{\gamma} = 1 + \sum_{1 \le i < n} \tilde{\gamma}_i (x - \xi)^i$ such that $\tilde{\gamma}^\nu \equiv \tilde{a} \bmod f$.

- If $\nu < m$, we define $\tilde{\alpha}(y) = y^\nu + y^m$. This time, we let $\tilde{\gamma}$ be the unique polynomial of the form $\tilde{\gamma} = 1 + \sum_{1 \le i < n} \tilde{\gamma}_i (x - \xi)^i$ such that $\tilde{\gamma}^\nu + (x - \xi)^{m-\nu} \tilde{\gamma}^m \equiv \tilde{a} \bmod f$. As previously, existence follows from Newton iteration, using the assumption $\nu \ne 0$ in $\mathbb{K}$.

In both cases, we set $\gamma = 1 + (x - \xi)\tilde{\gamma}$ rem $f \in \overline{\mathbb{K}}[x]$ and $\alpha = a_0 + a_\nu \tilde{\alpha}(y - 1) \in \overline{\mathbb{K}}[y]$. We can then verify that all requirements $\alpha(\gamma) \equiv a \bmod f$, $\gamma(\xi) \ne 0$ and $\gamma'(\xi) \ne 0$ are satisfied.

Since $\mu_\gamma$ then has degree $n$, and since $\mu_\gamma(0) = (-\gamma(\xi))^n$ is nonzero, Proposition 7.3 shows that $\text{Hk}_{m,d}^{(\alpha, \mu_\gamma)}$ has rank $n$ for $d \ge \lfloor n/m \rfloor$.

(6b) A polynomial in $\mathbb{K}[\bar{a}_0, \dots, \bar{a}_{n-1}, \bar{f}_0, \dots, \bar{f}_{n-1}]$. The existence of $\gamma$ and $\alpha$ implies that of a nonzero $n \times n$ minor $\delta$ of $\text{Hk}_{m, \lfloor n/m \rfloor}^{(\alpha, \chi_\gamma)}$. Let then $\Delta \in \mathbb{K}[\bar{a}_0, \dots, \bar{a}_{n-1}, \bar{f}_0, \dots, \bar{f}_{n-1}]$ be the corresponding
minor of $Hk_{m, \lfloor n/m \rfloor}^{(\bar{a}, \bar{f})}$, where $\bar{a} = \bar{a}_0 + \cdots + \bar{a}_{n-1}x^{n-1}$ and $\bar{f} = \bar{f}_0 + \cdots + \bar{f}_{n-1}x^{n-1} + x^n$ are polynomials whose coefficients are indeterminates. Lemma 7.7 shows that this is a polynomial of 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}$.

(6c) The rational functions $\bar{a}_0, \dots, \bar{a}_{n-1}$. Next, with $\bar{\gamma} = \bar{\gamma}_0 + \cdots + \bar{\gamma}_{n-1}x^{n-1}$ a polynomial whose coefficients are indeterminates, we consider $\bar{a}$ such that $\bar{a}(\bar{\gamma}) \equiv a \bmod f$. The coefficients of $\bar{a}$ are given as solutions of the linear system $\bar{a}(\bar{\gamma}) \equiv a \bmod f$, thus they are rational functions $\bar{a}_0, \dots, \bar{a}_{n-1}$ in $K(\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1})$. In this paragraph, we bound the degrees of their numerators and denominators in $K[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$, using power series inversion and composition.

We first consider the solution $u$ to $u(\bar{\gamma}) = x - \xi \bmod (x - \xi)^n$, or equivalently $u(\bar{\phi}) \equiv x \bmod x^n$, with $\bar{\phi} = \bar{\gamma}(x + \xi)$. We write $\bar{\phi} = \bar{\phi}_0 + \bar{\phi}_1x + \cdots + \bar{\phi}_{n-1}x^{n-1}$, where the coefficients $\bar{\phi}_0, \dots, \bar{\phi}_{n-1}$ are linear in $\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}$, with in particular $\bar{\phi}_0 = \Delta_0$ and $\bar{\phi}_1 = \Delta_1$. We can then write $u = \sum_{j=1}^{n-1} u_j(y - \Delta_0)^j$, where for $j \ge 1$, the coefficient $u_j$ is a rational function in $\bar{\gamma}_1, \dots, \bar{\gamma}_{n-1}$, with numerator of degree $j-1$ in $\bar{\gamma}_1, \dots, \bar{\gamma}_{n-1}$ and denominator $\Delta_1^{2j-1}$. More generally, for $i \ge 1$, the power $u^i$ has valuation $i$, and for $j \ge i$, the coefficient of $(y - \Delta_0)^j$ in it is a rational function with numerator of degree $j-i$ in $\bar{\gamma}_1, \dots, \bar{\gamma}_{n-1}$ and denominator $\Delta_1^{2j-i}$.

It follows that if we write $a = a_0 + a_\nu(x - \xi)^\nu + \cdots + a_{n-1}(x - \xi)^{n-1}$, then the solution $\bar{a}$ to the equation $\bar{a}(\bar{\gamma}) \equiv a \bmod f$ is given by $\bar{a} = a_0 + a_\nu u^\nu + \cdots + a_{n-1}u^{n-1} \bmod (y - \Delta_0)^n$. Once we rewrite $\bar{a}$ as $\bar{a}_0 + \cdots + \bar{a}_{n-1}y^{n-1}$, we see that the coefficients $\bar{a}_0, \dots, \bar{a}_{n-1}$ are rational functions with numerator of degree at most $2n-3$ in $\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}$, and denominator $\Delta_1^{2n-3}$.

(6d) The polynomial $\Delta_2$. We now evaluate the indeterminates $\bar{a}_i$ and $\bar{f}_i$ in the minor $\Delta$ of (6b) at the coefficients of $\bar{a}$ and $\chi_{\bar{\gamma}} = (y - \Delta_0)^n$, respectively. Write $(y - \Delta_0)^n$ as $\bar{q}_0 + \cdots + \bar{q}_{n-1}y^{n-1} + y^n$, so that $\bar{q}_i = \binom{n}{i}(-\Delta_0)^{n-i}$ for all $i$. It follows that $\Delta(\bar{a}_0, \dots, \bar{a}_{n-1}, \bar{q}_0, \dots, \bar{q}_{n-1})$ is a rational function in the indeterminates $\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}$, which can be written as

$$ \Delta(\bar{a}_0, \dots, \bar{a}_{n-1}, \bar{q}_0, \dots, \bar{q}_{n-1}) = \frac{\Delta_2(\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1})}{\Delta_1(\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1})^\epsilon}, \quad (33) $$

for some polynomial $\Delta_2$ of degree at most

$$ \frac{2n^2}{m}(2n-3) + \frac{2n^2(n-1)}{m}(n-1) = \frac{2n^2(n^2-2)}{m}, $$

and for some integer exponent $\epsilon \le 2n^2(2n-3)/m$.

Consider again the polynomials $\gamma$ and $\alpha$ in (6a), and their coefficients $\gamma_0, \dots, \gamma_{n-1}$ and $\alpha_0, \dots, \alpha_{n-1}$ (with actually $\alpha_{m+1} = \cdots = \alpha_{n-1} = 0$). We saw that $\gamma$ satisfies $\Delta_1(\gamma_0, \dots, \gamma_{n-1}) = \gamma'(\xi) \ne 0$, which implies that the rational functions $\bar{a}_0, \dots, \bar{a}_{n-1}$ are well-defined at $\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}$ and take $\alpha_0, \dots, \alpha_{n-1}$ for values there. This implies that the nonzero minor $\delta$ is $\delta = \Delta_2(\gamma_0, \dots, \gamma_{n-1})/\Delta_1(\gamma_0, \dots, \gamma_{n-1})^\epsilon$, and in particular that $\Delta_2$ is a nonzero polynomial.

*Probability bounds.* The end of the proof is as in the previous section. The polynomial $\Delta_0\Delta_1\Delta_f, m\Delta_2$ in $\overline{K}[\bar{\gamma}_0, \dots, \bar{\gamma}_{n-1}]$ has degree at most

$$ 1 + 1 + \frac{2n^2}{m} + \frac{2n^2(n^2-2)}{m} = \frac{2(n^4 - n^2 + m)}{m}; $$

we can now readily verify that a choice of $(r_3, \dots, r_{n+2})$ that avoids its zeros ensures that properties (3)-(6) hold. For (3)-(5), this follows immediately from the definitions.

To see that (6) holds, that is, that $Hk_{m,d}^{(\alpha, \mu_\gamma)}$ has rank $n$ for $d \ge \lfloor n/m \rfloor$, recall that the algorithm constructs $\gamma = r_3 + r_4x + \cdots + r_{n+2}x^{n+1}$. Properties (3)-(4) show that $\mu_\gamma$ has degree $n$, and that its constant coefficient is nonzero. Since in particular $\Delta_1(r_3, \dots, r_{n+2}) \ne 0$, we deduce that the rational
functions $\tilde{\alpha}_0, \dots, \tilde{\alpha}_{n-1}$ of (6c) are well-defined at $(r_3, \dots, r_{n+2})$, and that they give the coefficients of the unique polynomial $\alpha$ such that $\alpha(\gamma) \equiv a \mod f$. Since $\Delta_2(r_3, \dots, r_{n+2}) \neq 0$, it follows from Eq. (33) that $\mathrm{Hk}_{m, \lfloor n/m \rfloor}^{(\alpha, \mu_Y)}$ has rank $n$ (and thus similarly for $\mathrm{Hk}_{m, d}^{(\alpha, \mu_Y)}$, for $d \ge \lfloor n/m \rfloor$).

The other probabilities have been discussed in step (1)-(2) above and in step (7) of the previous section. Altogether, this gives a probability of success at least

$$ \left(1 - \frac{2(n^4 - n^2 + m)/m}{\mathrm{card}(S)}\right) \left(1 - \frac{m-1}{\mathrm{card}(S)}\right) \ge 1 - \frac{2(n^4 - n^2 + m)/m + m - 1}{\mathrm{card}(S)}. $$

Dividing the numerator of the last fraction by $n^4/m$ gives

$$ 2 - \frac{2n^2 - m^2 - m}{n^4} \le 2, $$

where the last inequality comes from $m \le n$. $\square$

**8.4.1 Note.** In Proposition 8.5, the role of the condition on the valuation being nonzero in $\mathbb{K}$ is shown by the following example. Take a field $\mathbb{K}$ of characteristic 2, $n = 6$, $m = 3$, $f = (x-1)^6$ and $a = (x-1)^2$. Then for any $\gamma \in \mathbb{K}[x]_{<6}$, the four polynomials $(1, a, \gamma^2, a\gamma^2)$ remain belong to the vector space generated by $(1, x^2, x^4)$ and are therefore linearly dependent. Using the expression of $M_\alpha$ from Eq. (9), we see that this implies that the block Krylov matrix $K_{m,n/m}^{(\alpha, \mu_Y)}$ of Eq. (11) is singular, and thus so is $\mathrm{Hk}_{m,n/m}^{(\alpha, \mu_Y)}$ regardless of the choice of $\gamma$.

A more general version of this counterexample when $\mathbb{K}$ has characteristic $p > 0$ is obtained with $m = p + 1$, $d = p$, $n = md$, and $\mathrm{val}_\xi(a) = p$.

## 8.5 Complete algorithm for $f$ purely inseparable

We now extend Proposition 8.5 in order to cover all cases of composition modulo a purely inseparable polynomial $f$.

If $p$ is the characteristic of $\mathbb{K}$, any purely inseparable $f$ can be written as $f(x) = (x^{p^e} - c)^\ell$ with $c$ in $\mathbb{K}$ and $\ell$ in $\mathbb{N}$ such that $p$ does not divide $\ell$, and $e = 0$ if $p = 0$ [25]; in particular the degree $n$ of $f$ is equal to $p^e\ell$. We assume that the parameters $e, \ell$ and $c$ are known, since this is the case when our algorithms have to handle this situation; indeed in the next section we introduce separable factorization techniques which allow us to compute them.

**8.5.1 Large valuation.** If $v = \mathrm{val}_\xi(a(x) - a(\xi))$ satisfies $v > \lfloor n^\eta \rfloor$, the minimal polynomial of $a$ in $\mathbb{K}[x]/\langle f \rangle$ factors over $\overline{\mathbb{K}}$ as $\mu_a(y) = (y - a(\xi))^\delta$, with $\delta = \lfloor n/v \rfloor \le \lfloor n^{1-\eta} \rfloor$. Since the latter degree is small compared to $n$, this case is handled efficiently by Algorithm MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL from Section 3.1.

**8.5.2 Small characteristic.** In the case $0 < p \le \lfloor n^\eta \rfloor$, our algorithm is based on Bernstein's composition algorithm for power series [6], which we adapt to work modulo $f(x) = (x^{p^e} - c)^\ell$. See also [30, Algorithm 3.1] for another extension of Bernstein's result, which is however not sufficient to reach our target cost for the specific kind of modulus we work with.

If $e = 0, p^e = 1$ and we are working modulo $f = (x - c)^\ell$, with $\ell = n$. In this case, to compute $b = g(a)$ remain $f$, we write $\tilde{a}(x) = a(x + c)$, we compute $\tilde{b} = g(\tilde{a})$ remain $x^\ell$, then we obtain $b$ as $\tilde{b}(x - c)$. The bottleneck is the computation of $g(\tilde{a})$ remain $x^\ell$, which can be done in $\tilde{O}(p\ell)$ operations in $\mathbb{K}$ using Bernstein's algorithm (in Algorithm MODULARCOMPOSITION-SMALLCHARACTERISTIC, that algorithm is called POWERSERIESCOMPOSITION-SMALLCHARACTERISTIC).
Suppose now that $e \ge 1$. Write $g = \sum_{i=0}^{p-1} g_i(y^p)y^i$, with $g_i \in \mathbb{K}[y]$ of degree less than $p^{e-1}\ell$. Write also $a(x) = \sum_{i=0}^{n-1} a_i x^i$, and let $\bar{a}(x) = \sum_{i=0}^{n-1} a_i^p x^i$, so that $a^p(x) = \bar{a}(x^p)$. It follows that

$$g(a) \operatorname{rem}(x^{p^e} - c)^\ell = \sum_{i=0}^{p-1} \bar{g}_i a^i \operatorname{rem}(x^{p^e} - c)^\ell$$

where, for all $0 \le i \le p-1$,

$$\bar{g}_i(x) = g_i(a^p(x)) \operatorname{rem}(x^{p^e} - c)^\ell = g_i(\bar{a}(x^p)) \operatorname{rem}(x^{p^e} - c)^\ell.$$

If we define $h_i = g_i(\bar{a}) \operatorname{rem}(x^{p^{e-1}} - c)^\ell$, it follows that $\bar{g}_i = h_i(x^p)$, so that

$$g(a) \operatorname{rem}(x^{p^e} - c)^\ell = \sum_{i=0}^{p-1} h_i(x^p) a^i \operatorname{rem}(x^{p^e} - c)^\ell.$$

The following lemma summarizes the cost of this procedure.

**Algorithm 8.3 MODULARCOMPOSITION-SMALLCHARACTERISTIC(c, e, ℓ, a, g)**

Input: $\mathbb{K}$ has characteristic $p > 0$,
$c \in \mathbb{K}$, $e \in \mathbb{N}$ and $\ell \in \mathbb{N}_{>0}$ such that $f = (x^{p^e} - c)^\ell$ has degree $n = \ell p^e$,
$a \in \mathbb{K}[x]_{<n}$, $g \in \mathbb{K}[y]_{<n}$

Output: $g(a) \operatorname{rem} f$

1: **if** $e = 0$ **then**

2: &nbsp;&nbsp;&nbsp;&nbsp;$a \leftarrow a(x + c)$

3: &nbsp;&nbsp;&nbsp;&nbsp;$b \leftarrow \text{POWERSERIESCOMPOSITION-SMALLCHARACTERISTIC}(x^n, a, g)$

▷ [6, Sec. 2]

4: **return** $b(x - c)$

5: **else**

6: &nbsp;&nbsp;&nbsp;&nbsp;Write $g = g_0(y^p) + \cdots + g_{p-1}(y^p)y^{p-1}$

7: &nbsp;&nbsp;&nbsp;&nbsp;Write $a = a_0 + \cdots + a_{n-1}x^{n-1}$

8: &nbsp;&nbsp;&nbsp;&nbsp;$\bar{a} \leftarrow a_0^p + \cdots + a_{n-1}^p x^{n-1}$

9: **for** $i = 0, \dots, p-1$ **do**

10: &nbsp;&nbsp;&nbsp;&nbsp;$h_i \leftarrow \text{MODULARCOMPOSITION-SMALLCHARACTERISTIC}(c, e-1, \ell, \bar{a}, g_i)$

11: **return** $h_0(x^p) + \cdots + h_{p-1}(x^p)a^{p-1} \operatorname{rem} f$

**LEMMA 8.6.** For a field $\mathbb{K}$ of characteristic $p > 0$, given a purely inseparable polynomial $f = (x^{p^e} - c)^\ell$ of degree $n = \ell p^e$, $a \in \mathbb{K}[x]_{<n}$ and $g \in \mathbb{K}[y]_{<n}$, Algorithm **MODULARCOMPOSITION-SMALLCHARACTERISTIC** returns $g(a) \operatorname{rem} f$ and uses $\tilde{O}(pn)$ operations in $\mathbb{K}$.

**PROOF.** Correctness follows from the previous description. For the runtime analysis when $e = 0$, the result is Bernstein's. For $e > 0$, apart from the $p$ recursive calls, Step 8 takes $\tilde{O}(n)$ operations (we raise all coefficients of $a$ to the power $p \le n$) and Step 11 takes $\tilde{O}(pn)$ operations, using Horner's rule. Remembering that $n = \ell p^e$, we deduce that the runtime $T(e, p, \ell)$ satisfies $T(e, p, \ell) = pT(e - 1, p, \ell) + \tilde{O}(p^{e+1}\ell)$ and $T(0, p, \ell) \in \tilde{O}(p\ell)$. This resolves to $T(e, p, \ell) \in \tilde{O}(p^{e+1}\ell)$, which is $\tilde{O}(pn)$. □

**8.5.3 Main algorithm.** Combining the previous results gives Algorithm **COMPOSITIONMODULOINSEPARABLE**. It first tests whether the characteristic of $\mathbb{K}$ is small enough for Algorithm
**MODULARCOMPOSITION-SMALLCHARACTERISTIC** to run within our prescribed runtime. Otherwise, rather than computing the valuation $v$, it simply calls Algorithm **MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL**; in case of failure, it falls back on Algorithm **MODULARCOMPOSITION-BASECASE**. As previously, the algorithm takes as input a vector $r$ that plays the role of random parameters.

<caption>**Algorithm 8.4 COMPOSITIONMODULOINSEPARABLE(c, e, ℓ, a, g, r)**</caption>

Input: $c$ in $\mathbb{K}$, $e$ in $\mathbb{N}$ and $\ell$ in $\mathbb{N}_{>0}$ such that $f = (x^{p^e} - c)^\ell$ has degree $n = \ell p^e$, where $p$ is the characteristic of $\mathbb{K}$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[y]_{<n}$, $r \in \mathbb{K}^{n+\lfloor n^\eta \rfloor}$ with $\eta$ from Eq. (3)

Output: $b = g(a)$ rem $f$ or **FAIL**

1: $n \leftarrow \ell p^e$

2: **if** $0 < p \le \lfloor n^\eta \rfloor$ **then return** **MODULARCOMPOSITION-SMALLCHARACTERISTIC**($c$, $e$, $\ell$, $a$, $g$)

3: $f \leftarrow (x^{p^e} - c)^\ell$

4: $b \leftarrow \text{MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL}(f, a, g, \lfloor n^{1-\eta} \rfloor, (r_i)_{0 \le i < n})$
    **if** $b \ne \text{FAIL}$ **then return** $b$

5: **if** $\gcd(a, f) = 1$ **then** $r_1 = 0$ **else** $r_1 = 1$; **if** $c \ne 0$ **then** $r_2 = 0$ **else** $r_2 = 1$

6: **return** **MODULARCOMPOSITIONBASECASE**($f, a, g, r$)

► **Proposition 8.5**

**PROPOSITION 8.7.** For a field $\mathbb{K}$ of characteristic $p$, given $c, e, \ell$ such that $f = (x^{p^e} - c)^\ell$ is purely inseparable of degree $n = \ell p^e$ ($e = 0$ if $p = 0$), $a \in \mathbb{K}[x]_{<n}$, $g \in \mathbb{K}[y]_{<n}$ and $r \in \mathbb{K}^{n+m}$ with $m = \lfloor n^\eta \rfloor$ and $\eta$ from Eq. (3), Algorithm **COMPOSITIONMODULOINSEPARABLE** uses $\tilde{O}(n^\kappa)$ operations in $\mathbb{K}$, with $\kappa < 1.43$ as in Eq. (1), and returns either $g(a)$ rem $f$, or **FAIL**.

If the entries of $r$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then the algorithm returns $g(a)$ rem $f$ with probability at least $1 - 2n^4/\text{card}(S)$.

**PROOF.** We first analyze the runtime. If $0 < p \le \lfloor n^\eta \rfloor$, then Algorithm **MODULARCOMPOSITION-SMALLCHARACTERISTIC** has cost $\tilde{O}(pn)$ by Lemma 8.6, which is thus $\tilde{O}(n^{1+\eta}) \in \tilde{O}(n^{1+(\omega-1)\eta}) = \tilde{O}(n^\kappa)$ from Eq. (3). Computing $f$ takes time $\tilde{O}(n)$ by repeated squaring. By Lemma 3.2, the call to Algorithm **MODULARCOMPOSITION-SMALLMINIMALPOLYNOMIAL** uses

$$ \tilde{O}\left(n^{1+(1-\eta)(\omega_2/2-1)}\right) = \tilde{O}\left(n^{\eta+(1-\eta)(\omega_2/2)}\right) = \tilde{O}(n^\kappa) $$

operations in $\mathbb{K}$, and by Proposition 8.1, it is also the case for Algorithm **MODULARCOMPOSITION-BASECASE**. The specifications of the subroutines imply that the output can be either $g(a)$ rem $f$ or **FAIL**, so only the probability analysis remains.

If $0 < p \le \lfloor n^\eta \rfloor$, Lemma 8.6 shows that the output is $g(a)$ rem $f$; hence, we may now assume that $p > \lfloor n^\eta \rfloor$, or $p = 0$. Let $\xi = c^{1/p^e} \in \overline{\mathbb{K}}$, so that $f = (x - \xi)^n$ in $\overline{\mathbb{K}}[x]$; let further $v$ be the valuation of $a - a(\xi)$ at $\xi$. The minimal polynomial of $a$ modulo $f$ has degree $\delta = \lfloor n/v \rfloor$.

Suppose first that $v \le \lfloor n^\eta \rfloor$, so that $\delta \ge \lfloor n^1/\lfloor n^\eta \rfloor \rfloor$. The value $b$ computed at Step 4 is either $g(a)$ rem $f$, or **FAIL**; let $\pi$ be the probability of the former (for instance, by Lemma 3.2, $\pi = 0$ if $\delta > \lfloor n^{1-\eta} \rfloor$). If **FAIL** is returned at Step 4, then we enter Step 6. At this stage, we have inequalities $v \le \lfloor n^\eta \rfloor < p$, or $v \le \lfloor n^\eta \rfloor$ and $p = 0$, so by Proposition 8.5 the call to Algorithm **MODULARCOMPOSITIONBASECASE** returns $g(a)$ rem $f$ with probability at least $1 - 2n^4/\text{card}(S)$. Overall, the probability of returning $g(a)$ rem $f$ in this case is at least $\pi + (1 - \pi)(1 - 2n^4/\text{card}(S))$, which is at least $1 - 2n^4/\text{card}(S)$.

Suppose on the other hand that $v > \lfloor n^\eta \rfloor$, so that we have in particular $v \ge n^\eta$, and thus $\delta = \lfloor n/v \rfloor \le \lfloor n^{1-\eta} \rfloor$. By Lemma 3.2, $b$ computed at Step 4 is $g(a)$ rem $f$ with probability at least
1 − n / card(S). If it is not the case, the algorithm enters Algorithm MODULARCOMPOSITIONBASECASE, which computes g ( a ) rem f with a certain probability π ′ ≥ 0. Overall, we return g ( a ) rem f with probability at least 1 − n / card(S) + π ′ ≥ 1 − n / card(S). □

9 ALGORITHM FOR GENERAL f

We now present our Las Vegas Algorithm MODULARCOMPOSITION that computes g ( a ) rem f for arbitrary input g , a , f . The analysis of this algorithm in Section 9.5 proves Theorem 1.1.

The starting point is the separable decomposition of f (Section 9.1), a generalization of square-
free decomposition from fields of characteristic zero to arbitrary base fields. This yields a partial
factorization f = f₁ · · · fₛ into pairwise coprime factors. The algorithm then proceeds by computing
g (a) modulo each of these factors and the final result is obtained by Chinese remaindering in quasi-
linear complexity [22, §10.3].

If p is the characteristic of $\mathbb{K}$ then the factors $f_i$ of the separable decomposition of $f$ are the
form $h_i(x^{p^{e_i}})^{\ell_i}$ (or more simply $h_i(x)^{\ell_i}$ when $p=0$), with integers $e_i$, $\ell_i$ and separable $h_i \in \mathbb{K}[x]$.
Composition modulo such an $f_i$ is achieved via a $\mathbb{K}$-algebra isomorphism

$$
\Psi_i : \mathbb{A}_i = \mathbb{K}[x] / \langle f_i(x) \rangle \to \mathbb{B}_i = \mathbb{K}[\theta, z] / \langle h_i(\theta), (z^{p^{e_i}} - \theta)^{\ell_i} \rangle
$$

that maps $x$ to $z$ (Proposition 9.6). If $\mathbb{L}_i$ denotes $\mathbb{K}[\theta]/\langle h_i(\theta) \rangle$, then, as a $\mathbb{K}$-vector space, $\mathbb{B}_i \simeq \mathbb{L}_i[z]/\langle(z^{p^{e_i}} - \bar{\theta}_i)^{\ell_i} \rangle$ with $\bar{\theta}_i$ the class of $\theta$ in $\mathbb{L}_i$. The computation of $g(a)$ rem $f_i$ over $\mathbb{K}$ is thus mapped to the composition

$$
g(A_i) \bmod (z^{p^{e_i}} - \bar{\theta}_i)^{\ell_i}
$$

over $\mathbb{L}_i$, with $A_i = \Psi_i(a \bmod f_i)$ and modulo the purely inseparable $(z^{p^{e_i}} - \bar{\theta}_i)^{\ell_i}$. In order to perform this last composition efficiently, it is also necessary to decrease the degree of $g$ by first reducing $g$ modulo the characteristic polynomial of $A_i$ in $\mathbb{L}_i[z]/\langle(z^{p^{e_i}} - \bar{\theta}_i)^{\ell_i} \rangle$. We call reduction of $g$ that step of the process (Proposition 9.8). It produces a representative of $G_i \in \mathbb{L}_i[y]$ such that $B_i = g(A_i) \in \mathbb{B}_i$ is obtained through the univariate modular composition

$$
G_i(A_i) \bmod (z^{p^{e_i}} - \bar{\theta}_i)^{\ell_i},
$$

which is computed with coefficients in $\mathbb{L}_i$. Finally, the class $g(a) \bmod f_i \in \mathbb{A}$ is recovered as
$\Psi_i^{-1}(B_i)$. In practice, the algorithms working with elements of $\mathbb{L}_i$ use polynomial representatives in
$\mathbb{K}[\theta]_{<\deg(f_i)}$, that are the canonical lifts of their class.

The idea of using these homomorphisms was introduced by van der Hoeven and Lecerf in the
case $e_i = 0$ [30]; it is extended to the general case in Sections 9.3 and 9.4. We keep their terminology,
calling *untangling* an algorithm that computes the map $\Psi_i$ and *tangling*, one which computes the
reverse map. Both these operations can be performed efficiently (Section 9.3).

The univariate modular composition in $\mathbb{L}_i[z]$ modulo the purely inseparable polynomial $(z^{p^{e_i}} - \bar{\theta})^{\ell_i}$ can be achieved by Algorithm COMPOSITIONMODULOINSEPARABLE of Section 8.5 when $\mathbb{L}_i$ is a field. In general however, $\mathbb{L}_i$ is a product of fields. In Section 9.2, the extension of the scope of our algorithms to this setting is obtained using a paradigm also due to van der Hoeven and Lecerf called *directed evaluation* [33].

Conventions. For $h$ of degree $d$ in $\mathbb{K}[\theta]$ and $f$ in $\mathbb{K}[\theta, z]$, monic of degree $n$ in $z$, and for any $P$
in $\mathbb{K}[\theta, z]$, we denote by $P \rem \langle h, f \rangle \in \mathbb{K}[\theta, z]_{<(d,n)}$ the polynomial obtained by reducing $P$ first
by $f$, then by $h$ (this is the normal form of $P$ modulo $(h, f)$, if we see the latter as a Gröbner basis
for the lexicographic order induced by $\theta < z$). Thus $P \rem \langle h, f \rangle$ is a canonical lift of the class of $P$
in $\mathbb{K}[\theta, z]/\langle h, f \rangle$.

If $P \in \mathbb{K}[\theta, z]$, we use the notation $\bar{P}(z)$ to denote the class (projection) of $P$ in $\mathbb{L}[z]$, where $\mathbb{L}$ will
be clear from the context.
## 9.1 Separable decomposition

Let $p$ be the characteristic of the field $\mathbb{K}$ and let $f$ in $\mathbb{K}[x]$ be of degree $n$. The separable decomposition of $f$ is the set

$$ S = \{(h_1, e_1, \ell_1), \dots, (h_s, e_s, \ell_s)\}, \quad \text{with } h_i \in \mathbb{K}[x] \text{ and } e_i, \ell_i \in \mathbb{N} \text{ for all } i, $$

that satisfies the following properties, where we write $f_i = h_i (x^{p^{e_i}})^{\ell_i}$:

(1) $f = c f_1 \cdots f_s$ with $c \in \mathbb{K} \setminus \{0\}$;

(2) for all $i \neq j$ in $\{1, \dots, s\}$, $f_i$ and $f_j$ are coprime;

(3) for all $i$ in $\{1, \dots, s\}$, $h_i \in \mathbb{K}[x]$ is separable, monic and of positive degree $d_i$;

(4) for all $i$ in $\{1, \dots, s\}$, $e_i = 0$ (if $p=0$) or $e_i$ is in $\mathbb{N}$ (if $p > 0$);

(5) for all $i$ in $\{1, \dots, s\}$, $\ell_i$ is not divisible by $p$;

(6) for all $i \neq j$ in $\{1, \dots, s\}$, $(e_i, \ell_i) \neq (e_j, \ell_j)$.

The separable decomposition of $f$ can be computed in $\tilde{O}(n)$ operations in $\mathbb{K}$ using an algorithm due to Lecerf [53]. The special case when $p=0$ recovers the more classical square-free factorization.

## 9.2 Composition over products of fields, modulo purely inseparable $f$

Let $h$ be separable of degree $d$ in $\mathbb{K}[\theta]$, and consider $f$ of the form $f = (z^{p^e} - c(\theta))^{\ell} \in \mathbb{K}[\theta, z]$, for integers $e \in \mathbb{N}$ and $\ell \in \mathbb{N}_{>0}$, where $p$ is the characteristic of $\mathbb{K}$. Given $A$ in $\mathbb{K}[\theta, z]_{<(d,n)}$ and $G$ in $\mathbb{K}[\theta, y]_{<(d,n)}$, with $n = \deg_G(f) = \ell p^e$, we consider here the computation of $B = G(\theta, A)$ rem $\langle h, f \rangle$. This question is mapped to a univariate composition problem with coefficients in $\mathbb{L} = \mathbb{K}[\theta]/\langle h \rangle$: if we let $\tilde{A}, \tilde{G}, \tilde{B}$ and $\tilde{c}$ be the projections of respectively $A, G, B$ and $c$ in $\mathbb{L}[z], \mathbb{L}[y], \mathbb{L}[z]$ and $\mathbb{L}$ (the degree constraints show that $\tilde{A}, \tilde{G}, \tilde{B}$ can be obtained without any calculation from $A, G, B$, and conversely), then $\tilde{B} = \tilde{G}(\tilde{A}) \text{ rem } (z^{p^e} - \tilde{c})^{\ell}$ as an equality in $\mathbb{L}[z]$.

When $h$ is irreducible, so that $\mathbb{L}$ is a field, the algorithm of Section 8.5 applies over $\mathbb{L}$; as reported in Proposition 8.7, if $n = \deg(f) = \ell p^e$, the runtime is $\tilde{O}(dn^\kappa)$ operations in $\mathbb{K}$, coming from $\tilde{O}((\ell p^e)^\kappa) = \tilde{O}(n^\kappa)$ times a factor in $\tilde{O}(d)$ for the cost of arithmetic operations in $\mathbb{L}$. However, we only assume $h$ separable, so that $\mathbb{L}$ is a product of fields. The key difference is the presence of zero-divisors in $\mathbb{L}$: a nonzero element of $\mathbb{L}$ is not necessarily invertible. Since the procedures in Section 8.5 use zero-tests and divisions, their direct application is not possible.

### 9.2.1 Directed evaluation.
The technique of directed evaluation, due to van der Hoeven and Lecerf [33], is an efficient version of the classical dynamic evaluation process [19].

In dynamic evaluation, prior to each zero-test or inversion, say by a quantity $q \in \mathbb{L}$, the computation of $h_1 = \gcd(q, h)$ gives the factorization $h = h_1 h_2$. Since $h$ is separable, $h_1$ and $h_2$ are coprime, and $\mathbb{L}$ can be decomposed as the product $\mathbb{L}_1 \times \mathbb{L}_2$, with $q = 0$ in $\mathbb{L}_1 = \mathbb{K}[\theta]/\langle h_1 \rangle$ and $q$ invertible in $\mathbb{L}_2 = \mathbb{K}[\theta]/\langle h_2 \rangle$. Under the dynamic evaluation paradigm, the calculation can then be continued in two branches, working modulo $h_1$ and $h_2$ separately.

In directed evaluation, the idea is rather to run the entire program in a unique branch, then to apply the process recursively in residual branches after reduction of input data modulo the corresponding polynomial. We do not detail the underlying techniques, for which we refer to Sections 3 and 4 of [33], and simply apply their *panoramic evaluation* procedure [33, Algo. 2]. It takes as input a computation tree $\mathcal{T}$ over $\mathbb{K}$ (see Section 2), a defining separable polynomial $h$ of degree $d$ for $\mathbb{L}$, and $\lambda = (\lambda_1, \dots, \lambda_s)$ in $\mathbb{K}[\theta]_<^s$ (representing an input to $\mathcal{T}$ in $\mathbb{L}^s$); it then returns a *panoramic value*, defined as follows.

**Definition 9.1 ([33, Def. 1 and Lem. 2]).** Given an input $(h, \lambda, \mathcal{T})$ as above, a *panoramic value* of $\mathcal{T}$ at $\lambda$ is a set of pairs $\{(h_i, \varepsilon_i), \dots, (h_t, \varepsilon_t)\}$, where
* $h_1, \dots, h_t$ are polynomials in $\mathbb{K}[\theta]$ that satisfy $h = h_1 \cdots h_t$ (thus $\mathbb{L} \simeq \mathbb{L}_1 \times \cdots \times \mathbb{L}_t$, with $\mathbb{L}_i = \mathbb{K}[\theta]/\langle h_i \rangle$);

* for all $i$, $\epsilon_i$ is a vector in $\mathbb{K}[\theta]^{\ell_i}_{<d_i}$ (representing an output in $\mathbb{L}_i^{\ell_i}$), with $d_i = \deg(h_i)$ and $\ell_i$ in $\mathbb{N}$;

* for all $1 \le i \le t$, let $h_{i,1}, \dots, h_{i,k_i}$ be the factorization of $h_i$ into irreducibles. For $1 \le j \le k_i$, let $\mathbb{L}_{i,j}$ be the the field $\mathbb{K}[\theta]/\langle h_{i,j} \rangle$, and denote by $\pi_{i,j} : \mathbb{K}[\theta] \to \mathbb{L}_{i,j}$ the canonical projection $a \mapsto a \bmod h_{i,j}$ (the notation carries over to vectors over $\mathbb{K}[\theta]$). Then $\mathcal{T}$ is supposed to be evaluable at $\pi_{i,j}(\lambda) \in \mathbb{L}_{i,j}^s$ for all $i, j$, and $\pi_{i,j}(\epsilon_i) \in \mathbb{L}_{i,j}^{\ell_i}$ is the result of evaluating $\mathcal{T}$ (seen as a computation tree over $\mathbb{L}_{i,j}$) at $\pi_{i,j}(\lambda)$, using the same branch of $\mathcal{T}$ for all $j$.

The application of this method requires that one uses computation trees as the underlying computational model, which is the case here (Section 2). Crucially, the cost overhead is then $\tilde{O}(d)$ [33, Thm. 1], i.e. similar (up to logarithmic factors) to the one incurred if $h$ were irreducible.

**9.2.2 Algorithm.** With Algorithm **COMPOSITIONMODULOINSEPARABLE-PRODUCTOFFIELDS** we apply panoramic evaluation (called **PANORAMIC** in our pseudo-code) to Algorithm **COMPOSITIONMODULOINSEPARABLE** for modular composition over $\mathbb{K}$. Note that in addition to field elements, the latter algorithm also takes two integers $e, \ell$ as input. Panoramic evaluation can still be used in this context, since each choice of the parameters $e, \ell$ corresponds to a computation tree, to which the techniques described above apply. This yields a factorization of $h$, and performs the compositions modulo the corresponding factors; the final result is then reconstructed using Chinese remaindering.

<caption>**Algorithm 9.1 COMPOSITIONMODULOINSEPARABLE-PRODUCTOFFIELDS(h, c, e, ℓ, A, G, r)**</caption>

**Input:** $h$ separable of degree $d$ in $\mathbb{K}[\theta]$, $c$ in $\mathbb{K}[\theta]_{<d}$, $e$ in $\mathbb{N}$ and $\ell$ in $\mathbb{N}_{>0}$ such that $f = (z^{p^e} - c)^\ell$ has degree $n = \ell p^e$, $A \in \mathbb{K}[\theta, z]_{<(d,n)}$, $G \in \mathbb{K}[\theta, y]_{<(d,n)}$, $r \in \mathbb{K}^{n+[n^\eta]}$
**Output:** $B = G(\theta, A)$ rem $\langle h, f \rangle$, or **FAIL**
1: $\mathbf{Splitting}$ $\mathbb{L} \simeq \mathbb{K}[\theta]/\langle h_1 \rangle \times \cdots \times \mathbb{K}[\theta]/\langle h_t \rangle$ and reductions of $B$, accordingly, using [33, Algo. 2] $\{(h_1, B_1), \dots, (h_t, B_t)\} \leftarrow \mathbf{PANORAMIC}(\mathbf{COMPOSITIONMODULOINSEPARABLE}, h, c, e, \ell, A, G, r)$
2: **if** any of the $B_i$'s equals **FAIL** **then return** **FAIL**
3: **return** **CHINESEREMAINDERING**($B_1, \dots, B_t$), $(h_1, \dots, h_t)$)

**PROPOSITION 9.2.** For a field $\mathbb{K}$ of characteristic $p$, given $h \in \mathbb{K}[\theta]$ separable of degree $d$, $c$ in $\mathbb{K}[\theta]_{<d}$, integers $e$ in $\mathbb{N}$ and $\ell$ in $\mathbb{N}_{>0}$, $A$ in $\mathbb{K}[\theta, z]_{<(d,n)}$, $G$ in $\mathbb{K}[\theta, y]_{<(d,n)}$, $r$ in $\mathbb{K}^{n+[n^\eta]}$ with $n = \ell p^e$ and $\eta$ from Eq. (3), Algorithm **COMPOSITIONMODULOINSEPARABLE-PRODUCTOFFIELDS** uses $\tilde{O}(d(\ell p^e)^\kappa) = \tilde{O}(dn^\kappa)$ operations in $\mathbb{K}$, with $\kappa < 1.43$ as in Eq. (1).

It returns either $G(\theta, A)$ rem $\langle h, f \rangle \in \mathbb{K}[\theta, z]_{<(d,n)}$ or **FAIL**, with $f = (z^{p^e} - c)^\ell$. If the entries of $r$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then the algorithm returns $G(\theta, A)$ rem $\langle h, f \rangle$ with probability at least $1 - 2dn^4/\text{card}(S)$.

The complexity bound $\tilde{O}(dn^\kappa)$ indicates that the overhead coming from operations modulo $h(\theta)$ is $\tilde{O}(d)$, as pointed out above.

**9.2.3 Proof of Proposition 9.2.** Combined with our Proposition 8.7, Theorem 1 in [33] gives the runtime estimate. In the pseudo-code, the output of the panoramic evaluation is written as $\{(h_1, B_1), \dots, (h_t, B_t)\}$, where $h_1 \cdots h_t$ is a factorization of $h$ (not necessarily into irreducibles), and for all $i$, either $B_i \in \mathbb{K}[\theta, z]_{<(d_i, n)}$ with $d_i = \deg(h_i)$, or $B_i = \text{FAIL}$. At the level of computation trees, a flag such as **FAIL** is obtained by setting a dedicated output value to 1 (and 0 otherwise); call
flag$_i$ this value, for $1 \le i \le t$. If flag$_i$ = 1 (failure), we set $B_i = 0$ by convention, so in the rest of this proof, $B_i$ is an element of $\mathbb{K}[\theta, z]$ for all $i$.

We use the following notation: for $1 \le i \le t$, the irreducible factors of $h_i$ are written $h_{i,1}, \dots, h_{i,k_i}$. For $1 \le j \le k_i$, we then define $\bar{c}_{i,j}, \bar{A}_{i,j}, \bar{G}_{i,j}$ by taking $c, A, G$ modulo $h_{i,j}$ and seeing them over the field $\mathbb{L}_{i,j} = \mathbb{K}[\theta]/\langle h_{i,j} \rangle$, so $\bar{c}_{i,j}$ is in $\mathbb{L}_{i,j}$, $\bar{A}_{i,j}$ in $\mathbb{L}_{i,j}[z]$ and $\bar{G}_{i,j}$ in $\mathbb{L}_{i,j}[y]$. The elements in the vector $r$ are already in $\mathbb{K}$, and thus in $\mathbb{L}_{i,j}$.

Finally, we let $\tilde{B}_{i,j}$ be the polynomial obtained by taking $B_i \in \mathbb{K}[\theta, z]$ and projecting it to $\mathbb{L}_{i,j}[z]$ through reduction modulo $h_{i,j}$, and we set flag$_{i,j}$ = flag$_i$ (recall that flag$_i \in \mathbb{K}$ is either 0 or 1).

Then, from Definition 9.1, the key property of the output of the first step is that for all indices $i, j$, flag$_{i,j}$ and $\tilde{B}_{i,j}$ are the result of calling Algorithm COMPOSITIONMODULOINSEPARABLE on input $\bar{c}_{i,j}, e, \ell, \bar{A}_{i,j}, \bar{G}_{i,j}, r$ over the field $\mathbb{L}_{i,j}$. This implies in particular that our algorithm returns FAIL if and only if the computation fails over one of the fields $\mathbb{L}_{i,j}$.

To quantify the probability of this event, we apply Proposition 8.7 over all fields $\mathbb{L}_{i,j}$. For any given $i, j$, Proposition 8.7 shows that flag$_{i,j}$ = 1 occurs with probability at most $2n^4/\text{card}(S)$. Since there are at most $d$ such indices $i, j$, the probability that this happens for at least one pair of indices is at most $2dn^4/\text{card}(S)$.

Assume none of the flag$_{i,j}$'s is 1, so that the algorithm does not return FAIL. Then, for all $i, j$, $\tilde{B}_{i,j} \in \mathbb{L}_{i,j}[z]_{<n}$ is equal to $\tilde{G}_{i,j}(\bar{A}_{i,j})$ rem $(z^{p^e} - \bar{c}_{i,j})^\ell$. In terms of bivariate polynomials, the Chinese Remainder Theorem then implies that for all $i, B_i$ itself is equal to $G(\theta, A)$ rem $\langle h_i, (z^{p^e} - c)^\ell \rangle \in \mathbb{K}[\theta]_{<(d_i, n)}$. In the last step of the algorithm, we further apply the Chinese Remainder Theorem coefficient-wise to the $B_i$'s with respect to $z$; this gives us $G(\theta, A)$ rem $\langle h, (z^{p^e} - c)^\ell \rangle$ as a polynomial in $\mathbb{K}[\theta, z]_{<(d, n)}$. The cost of this last step is in $\tilde{O}(d\ell p^e)$, so the proof is complete.

## 9.3 Untangling and tangling

In this subsection, we give the main tools (tangling, untangling and bivariate reduction) that are needed for reducing composition modulo powers of separable polynomials to the situation of the previous subsection. The central results are due to van der Hoeven and Lecerf [30] with $f = h(x)^\ell$ and $h$ separable (Sections 9.3.1 and 9.3.2). We slightly generalize them to the case $f = h(x^{p^e})^\ell$ with $e > 0$ (Sections 9.3.3 and 9.3.4).

### 9.3.1 Tangling and untangling. The starting point is the following observation.

**LEMMA 9.3 ([30, §4.2]).** For $h$ of degree $d$ in $\mathbb{K}[x]$ and for a positive integer $\ell$, there exists a $\mathbb{K}$-algebra homomorphism

$$ \psi_{h,\ell} : \mathbb{K}[x]/\langle h(x)^{\ell} \rangle \to \mathbb{K}[\theta, z]/\langle h(\theta), (z-\theta)^{\ell} \rangle \\ x \mapsto z. $$

If moreover $h$ is separable then $\psi_{h,\ell}$ is an isomorphism.

This homomorphism is a variant of the homomorphism $\pi_{h,\ell}$ considered by van der Hoeven and Lecerf, that maps $u \in \mathbb{K}[x]/\langle h(x)^\ell \rangle$ to $u(z + \theta) \in \mathbb{K}[\theta, z]/\langle h(\theta), z^\ell \rangle$. The morphism $\psi_{h,\ell}$ is obtained by composing $\pi_{h,\ell}$ with a translation $z \mapsto z - \theta$. It turns out that $\psi_{h,\ell}$ is more convenient than $\pi_{h,\ell}$ for our generalization in Section 9.3.3. van der Hoeven and Lecerf call UNTANGLING($h, \ell, u$) the algorithm which implements $\pi_{h,\ell}$; we use that terminology for the algorithm which implements $\psi_{h,\ell}$: given $u$ in $\mathbb{K}[x]_{<d\ell}$, it computes $U \in \mathbb{K}[\theta, z]_{<(d, \ell)}$ such that $U = u(z)$ rem $\langle h(\theta), (z - \theta)^\ell \rangle$. When $h$ is separable, the inverse operation is called TANGLING($h, \ell, U$). Again, we use their terminology for the inverse of $\psi_{h,\ell}$.

**LEMMA 9.4.** UNTANGLING and TANGLING (when defined) take $\tilde{O}(d\ell)$ operations in $\mathbb{K}$.
PROOF. This is mostly in [30]. First, it is easy to check that the algorithms 4.3 and 4.5 and the proofs of Prop. 4.6 and 4.10 of that reference do not make use of the separability of $h$. Next, translation can be performed in quasi-linear complexity over an arbitrary ring [24, Thm. 4.5], so that the complexity estimate is unchanged for our variant of these algorithms. □

### 9.3.2 Bivariate reduction.
The computation of the composition $g(a)$ rem $h(x)^\ell$ for a separable $h$ reduces to computing $\psi_{h,\ell}^{-1}(g(\psi_{h,\ell}(a \mod h(x)^\ell)))$, where the inner composition is performed as a univariate composition in $\mathbb{L}[z]$ modulo $(z - \bar{\theta})^\ell$, with $\mathbb{L} = \mathbb{K}[\theta]/\langle h \rangle$.

In order to make use of the algorithms of the previous sections to perform this composition, it is necessary to first reduce the degree of $g$. Denote by $\mathcal{A}$ the canonical lift of $\psi_{h,\ell}(a \mod h(x)^\ell)$, and by $\bar{\mathcal{A}}$ its projection in $\mathbb{L}[z]$. The idea is to reduce $g$ modulo the characteristic polynomial $(y - \bar{\mathcal{A}}(\bar{\theta}))^\ell \in \mathbb{L}[y]$ of $\bar{\mathcal{A}}(z)$ modulo $(z - \bar{\theta})^\ell$.

This is achieved in two steps. For $h$ of degree $d$, we let $\alpha \in \mathbb{K}[\theta]_{<d}$ be the canonical lift of $\bar{\mathcal{A}}(\bar{\theta}) \in \mathbb{L}$. First, one computes the canonical lift of $\psi_{\mu,\ell}(g \mod \mu^\ell)$, where $\mu$ is an annihilating polynomial of $\alpha \mod h$. This produces $\tilde{G}(z, y) \in \mathbb{K}[z, y]_{<\deg \mu, \ell}$ such that

$$ \tilde{G}(z, y) = \sum_{i=0}^{\ell-1} \tilde{G}_i(z)y^i = g(y) + \tilde{U}(z, y)\mu(z) + \tilde{V}(z, y)(y-z)^\ell $$

for some polynomials $\tilde{U}, \tilde{V}$ in $\mathbb{K}[z, y]$.

Next, in view of $\mu(\alpha) \equiv 0 \mod h$, a modular composition of each of the $\ell$ coefficients of this polynomial $\tilde{G}$ in $y$ with $\alpha(\theta)$ modulo $h(\theta)$ gives $G(\theta, y) \in \mathbb{K}[z, y]_{<\deg \mu, \ell}$ such that

$$ G(\theta, y) = g(y) + U(\theta, y)h(\theta) + V(\theta, y)(y - \alpha(\theta))^\ell, \quad (34) $$

for some polynomials $U, V$ in $\mathbb{K}[\theta, y]$. Eq. (34) may also be read as $\tilde{G}(\bar{\mathcal{A}}) = g(\bar{\mathcal{A}}) \text{ rem } (z - \bar{\theta})^\ell$ over $\mathbb{L}$.

These two steps are detailed in Algorithm BIVARIATEREDUCTION below and correspond to Steps (2)-(4) of [30, Algo. 4.2]. The runtime and probability analyses are new; they are based on the results of the previous sections.

<caption>**Algorithm 9.2 BIVARIATEREDUCTION(h, ℓ, α, g, r)**</caption>

Input: $h$ separable, monic, of degree $d$ in $\mathbb{K}[\theta]$, $\ell$ in $\mathbb{N}_{>0}$, $\alpha$ in $\mathbb{K}[\theta]_{<d}$, $g$ in $\mathbb{K}[y]$, $r$ in $\mathbb{K}^{d+\lfloor d\eta \rfloor}$

Output: $G(\theta, y) = g(y)$ rem $\langle h(\theta), (y - \alpha(\theta))^\ell \rangle \in \mathbb{K}[\theta, y]_{<\langle d, \ell \rangle}$, or FAIL

1: ▶ Either $\mu = \text{FAIL}$, or $\mu$ is nonzero in $\mathbb{K}[y]_{\le 4d}$ and $\mu(\alpha) \equiv 0 \mod h$
   $\mu \leftarrow \text{ANNIHILATINGPOLYNOMIAL}(h, \alpha, r)$
   ▶ Algorithm 8.2

if $\mu = \text{FAIL}$ then return FAIL

2: $\tilde{G} \leftarrow \text{UNTANGLING}(\mu, \ell, g \text{ rem } \mu^\ell)$
   ▶ $\tilde{G}(\gamma, z) \in \mathbb{K}[\gamma, y]_{<\langle \deg(\mu), \ell \rangle}$, Lemma 9.4

3: Write $\tilde{G} = \sum_{0 \le i < \ell} \tilde{G}_i(\gamma)y^i$
   ▶ $\tilde{G}_i \in \mathbb{K}[\gamma]_{<\deg(\mu)}$

4: **for** $i = 0, \dots, \ell - 1$ **do**
   $G_i \leftarrow \text{MODULARCOMPOSITIONBASECASE}(h, \alpha, \tilde{G}_i, r)$
   ▶ $G_i = \tilde{G}_i(\alpha)$ rem $h$ or FAIL, Algorithm 8.1

if $G_i = \text{FAIL}$ then return FAIL

5: $G \leftarrow \sum_{0 \le i < \ell} G_i y^i$
   ▶ $G$ is in $\mathbb{K}[\theta, y]_{<\langle d, \ell \rangle}$

6: **return** $G$

**LEMMA 9.5.** Given $h$ in $\mathbb{K}[\theta]$ monic, separable and of degree $d$, $\alpha$ in $\mathbb{K}[\theta]_{<d}$, $g$ in $\mathbb{K}[y]$, $r$ in $\mathbb{K}^{d+\lfloor d\eta \rfloor}$ with $\eta$ from Eq. (3), and $\ell$ in $\mathbb{N}_{>0}$, Algorithm BIVARIATEREDUCTION uses $\tilde{O}(\deg(g) + d^\kappa \ell)$ operations in $\mathbb{K}$ with $\kappa < 1.43$ as in Eq. (1), and returns either $g$ rem $\langle h(\theta), (y - \alpha(\theta))^\ell \rangle$ or FAIL. If the entries of $r$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then the algorithm returns $g$ rem $\langle h, (y - \alpha)^\ell \rangle$ with probability at least $1 - 6(\ell + 1)d^2/\text{card}(S)$.
PROOF. The reduction of $g \mod \mu^l$ is justified by the fact that $\mu(a)^l = 0 \mod h^l$. The correction of the rest of the algorithm when Step 6 is reached follows from the discussion above.

Since $h$ is separable, Proposition 8.3 applies; it shows that the first step computes an annihilating polynomial for $\alpha$ modulo $h$ with probability at least $1 - 6d^2/\text{card}(S)$. It also shows that each call to Algorithm MODULARCOMPOSITIONBASECASE succeeds with at least the same probability. Altogether, the probability of success of the whole algorithm is thus at least $1 - 6(\ell + 1)d^2/\text{card}(S)$.

By Corollary 8.2, the first step uses $\tilde{O}(d^\kappa)$ operations in $\mathbb{K}$. Since $\text{deg}(\mu)$ is in $O(d)$, computing $g$ mod $\mu^l$ takes $\tilde{O}(\text{deg}(g) + d\ell)$ operations in $\mathbb{K}$, and Lemma 9.4 shows that deducing $\tilde{G}$ takes a further $\tilde{O}(d\ell)$ cost. Finally, by Proposition 8.1, each pass in the loop at Step 4 takes $\tilde{O}(d^\kappa)$ operations, so that the overall runtime is $\tilde{O}(\text{deg}(g) + d^\kappa \ell)$. $\square$

9.3.3 General Tangling and Untangling. In fields of positive characteristic, the isomorphism of Lemma 9.3 and the complexity of its realization generalize as follows.

PROPOSITION 9.6. Let $f = h(x^{p^e})^l$ be of degree $n$, with $h$ of degree $d$ in $\mathbb{K}[x]$, and $\mathbb{K}$ of characteristic $p$ ($e = 0$ if $p = 0$). There exists a $\mathbb{K}$-algebra homomorphism

$$ \Psi_{h,\ell} : \mathbb{K}[x]/\langle f \rangle \to \mathbb{K}[\theta, z]/\langle h(\theta), (z^{p^e} - \theta)^l \rangle \\ x \mapsto z. $$

If moreover $h$ is separable then $\Psi_{h,\ell}$ is an isomorphism. Applying $\Psi_{h,\ell}$ or its inverse when the latter is defined takes quasi-linear time $\tilde{O}(n) = \tilde{O}(d\ell p^e)$ over $\mathbb{K}$.

PROOF. Write $\mathbb{A} = \mathbb{K}[x]/\langle f \rangle$ and $\mathbb{B} = \mathbb{K}[\theta, z]/\langle h(\theta), (z^{p^e} - \theta)^l \rangle$. When $h$ is separable, we prove that the minimal polynomial of $z$ in the $\mathbb{K}$-algebra $\mathbb{B}$ is $f$. This implies that $\mathbb{A}$ is $\mathbb{K}$-isomorphic (as a $\mathbb{K}$-algebra) to the sub-algebra of $\mathbb{B}$ generated by $z$. Since $\mathbb{B}$ has $\mathbb{K}$-dimension $n = \text{deg}(f)$, this sub-algebra is $\mathbb{B}$ itself, and the first claim will follow.

To determine the minimal polynomial of $z$, we can work in $\bar{\mathbb{K}} = \bar{\mathbb{K}}[\theta, z]/\langle h(\theta), (z^{p^e} - \theta)^l \rangle$, where $\bar{\mathbb{K}}$ is an algebraic closure of $\mathbb{K}$. If we let $\xi_1, \dots, \xi_d$ be the roots of $h$ in $\bar{\mathbb{K}}$ (which are pairwise distinct), then $\bar{\mathbb{B}}$ is isomorphic, as a $\bar{\mathbb{K}}$-algebra, to the product

$$ \bar{\mathbb{K}}[\theta, z]/\langle \theta - \xi_1, (z^{p^e} - \xi_1)^l \rangle \times \dots \times \bar{\mathbb{K}}[\theta, z]/\langle \theta - \xi_d, (z^{p^e} - \xi_d)^l \rangle. $$

The minimal polynomial of $z$ in the $i$th factor above is $\mu_i = (x^{p^e} - \xi_i)^l$ for $1 \le i \le d$. These polynomials are pairwise coprime: since $t \mapsto t^{p^e}$ is a bijection in $\bar{\mathbb{K}}$, $\mu_i$ has a unique root in $\bar{\mathbb{K}}$, which is the $p^e$-th root of $\xi_i$, and these roots are pairwise distinct, since the $\xi_i$'s are. As a result, the minimal polynomial of $z$ in $\bar{\mathbb{B}}$, or equivalently in $\mathbb{B}$, is the product $\mu_1 \cdots \mu_d = f$.

For the second claim, we take $a$ in $\mathbb{K}[x]$ of degree less than $n$, and write it as $a = \sum_{0 \le i < p^e} a_i(x^{p^e})x^i$, with all $a_i$'s of degree less than $n/p^e = d\ell$. Then,

$$ \begin{align} \Psi_{h,\ell}(a \mod f) &\equiv \sum_{0 \le i < p^e} a_i (z^{p^e})^i z^i \quad \mod \langle h(\theta), (z^{p^e} - \theta)^l \rangle, \notag \\ &\equiv \sum_{0 \le i < p^e} \tilde{A}_i(\theta, z^{p^e}) z^i \quad \mod \langle h(\theta), (z^{p^e} - \theta)^l \rangle, \tag{35} \end{align} $$

where $\tilde{A}_i(\theta, z) = a_i(z) \text{ rem } \langle h(\theta), (z - \theta)^l \rangle$ is in $\mathbb{K}[\theta, z]_{<(d,\ell)}$; these degree bounds show that the expression in Eq. (35) is indeed reduced modulo $\langle f(\theta), (z^{p^e} - \theta)^l \rangle$. Each $\tilde{A}_i = \psi_{h,\ell}(a_i)$ can be computed in time $\tilde{O}(d\ell)$ by Lemma 9.4, so that one application of $\Psi_{h,\ell}$ takes $\tilde{O}(d\ell p^e) = \tilde{O}(n)$ operations in $\mathbb{K}$, as claimed.
Conversely, any element $B$ in $\mathbb{K}[\theta, z]_{<(d, \ell p^e)}$ can be written as in Eq. (35), for some $\tilde{B}_i$'s in $\mathbb{K}[\theta, z]_{<(d, \ell p^e)}$. Applying $\psi_{h, \ell}^{-1}$ to each of them allows us to recover $b = \Psi_{h, \ell}^{-1}(B)$, by reversing the steps above. The cost analysis is similar to the one for $\Psi_{h, \ell}$. $\square$

We call UNTANGLING-GENERAL($h, e, \ell, a$) the algorithm outlined in this proof that applies $\Psi_{h, \ell}$ to (the class modulo $f$ of) $a \in \mathbb{K}[x]_{<n}$, and returns the canonical lift of $\Psi_{h, \ell}(a \text{ mod } f)$ to $\mathbb{K}[\theta, z]_{<(d, \ell p^e)}$; equivalently, $A(\theta, z) = a(z) \text{ rem } \langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$. For $B$ in $\mathbb{K}[\theta, z]_{<(d, \ell p^e)}$, the inverse operation is written TANGLING-GENERAL($h, e, \ell, B$).

### 9.3.4 Main reduction.
A more general form of bivariate reduction is needed in Section 9.4. With $h$ of degree $d$ as before, given $g$ in $\mathbb{K}[y]$ and now a bivariate $A$ in $\mathbb{K}[\theta, z]_{<(d, \ell p^e)}$, the aim is to reduce the degree of $g$ before performing the composition in $\mathbb{L}[z]$ modulo $(z^{p^e} - \bar{\theta})^\ell$ with $\mathbb{L} = \mathbb{K}[\theta]/\langle h \rangle$. Denoting by $\bar{A}$ the projection of $A$ in $\mathbb{L}[z]$, the idea is to compute $\tilde{G} = g \text{ rem } \chi_{\bar{A}}$ in $\mathbb{L}[z]$, where $\chi_{\bar{A}} \in \mathbb{L}[y]$ is the characteristic polynomial of $\bar{A} \in \mathbb{L}[z]$ in the extension $\mathbb{L} \to \mathbb{L}[z]/\langle (z^{p^e} - \bar{\theta})^\ell \rangle$. Thus, $\tilde{G} \in \mathbb{L}[y]$ has degree less than $\ell p^e$; its canonical lift $G \in \mathbb{K}[\theta, y]_{< (d, \ell p^e)}$ is the output.

The computation of $g$ rem $\chi_{\bar{A}}$ is made easy by an explicit formula for the characteristic polynomial $\chi_{\bar{A}}$. In the following lemma, we let $\sigma : \mathbb{L} \to \mathbb{L}$ be the $p^e$th-power operator; we write the image of $\Lambda \in \mathbb{L}$ as $\Lambda^\sigma$. This notation is extended to the coefficient-wise action on polynomial rings over $\mathbb{L}$.

**LEMMA 9.7.** *The characteristic polynomial of $\bar{A}$ relative to the extension $\mathbb{L} \to \mathbb{L}[z]/\langle (z^{p^e} - \bar{\theta})^\ell \rangle$ is* $\chi_{\bar{A}} = (y^{p^e} - \bar{\alpha})^\ell \in \mathbb{L}[y]$, where $\bar{\alpha} = \bar{A}^\sigma(\bar{\theta}) \in \mathbb{L}$.

**PROOF.** The characteristic polynomial $\chi_{\bar{A}}$ can be computed relative to the extension $\mathbb{L}^* \to \mathbb{L}^*[z]/\langle (z^{p^e} - \bar{\theta})^\ell \rangle$, where we set $\mathbb{L}^* = \mathbb{L}[w]/\langle w^{p^e} - \bar{\theta} \rangle$. In $\mathbb{L}^*[z]$, we have the factorization

$$ (z^{p^e} - \bar{\theta})^\ell = (z^{p^e} - w^{p^e})^\ell = (z - w)^{\ell p^e}, $$

so the characteristic polynomial of $\bar{A}$ in $\mathbb{L}^*[z]/\langle (z^{p^e} - \bar{\theta})^\ell \rangle$ is

$$ (y - \bar{A}(w))^{\ell p^e} = (y^{p^e} - \bar{A}(w))^{p^e} = (y^{p^e} - \bar{A}^{\sigma}(\bar{\theta}))^{\ell}. \quad \square $$

The reduction of $g$ by this characteristic polynomial is described in Algorithm MAINREDUCTION. First, the canonical lift $\alpha \in \mathbb{K}[\theta]_{<d}$ of $\bar{\alpha} \in \mathbb{L}$ from Lemma 9.7 is computed. Next, in Step 3, the polynomial $g$ is rewritten as a polynomial in $y$ of degree less than $p^e$, with coefficients $g_i(y^{p^e})$. Each of these polynomials $g_i(y)$ can then be reduced modulo $\langle h, (y-\alpha)^\ell \rangle$ by Algorithm BIVARIATEREDUCTION, producing a polynomial $G_i(\theta, y)$ (Step 4). Thus, $G_i(\theta, y) \equiv g_i(y) \text{ mod } \langle h, (y-\alpha)^\ell \rangle$, whence $G_i(\bar{\theta}, y^{p^e}) \equiv g_i(y^{p^e}) \text{ mod } \chi_{\bar{A}}$. Recombining these coefficients yields $G(\theta, y)$ such that $G(\bar{\theta}, y) \equiv g(y) \text{ mod } \chi_{\bar{A}}$. Finally, since $\chi_{\bar{A}}(\bar{A}) \equiv 0$ in $\mathbb{L}[z]/\langle (z^{p^e} - \bar{\theta})^\ell \rangle$, it follows that $G(\theta, A) \equiv g(A) \text{ mod } \langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$.

**PROPOSITION 9.8.** *Given $h$ separable, monic, of degree $d$ in $\mathbb{K}[x]$, $e$ in $\mathbb{N}$, $\ell$ in $\mathbb{N}_{>0}$, $A$ in $\mathbb{K}[\theta, z]_{< (d, \ell p^e)}$, $g$ in $\mathbb{K}[y]$, and $r$ in $\mathbb{K}^{d+\lceil d^n \rceil}$, Algorithm MAINREDUCTION uses $\tilde{O}(\deg(g) + n^\kappa)$ operations in $\mathbb{K}$, with $n = d\ell p^e$ and $\kappa < 1.43$ as in Eq. (1). It returns $G \in \mathbb{K}[\theta, y]_{< (d, \ell p^e)}$ such that $G(\theta, A) \equiv g(A) \text{ mod } \langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$, or FAIL.*

If the entries of $r$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then the algorithm returns $G$ with probability at least $1 - 6(\ell + 1)d^2p^e/\text{card}(S)$.

**PROOF.** The correction of the algorithm when it does not return FAIL follows from the discussion above.

Working coefficient-wise, since $e = O(\log(p^e))$ the computation of $\alpha$ at Step 2 takes $\tilde{O}(\ell p^e)$ operations on polynomials modulo $h$ of degree $d$, so $\tilde{O}(n)$ operations in $\mathbb{K}$; reducing it modulo $h$ has the same complexity bound. The cost is thus governed by the loop, which uses $\tilde{O}(\deg(g) + d^\kappa \ell p^e) = \tilde{O}(\deg(g) + (n/d)d^\kappa)$ operations by Lemma 9.5. The latter also allows us to quantify the probability
**Algorithm 9.3 MAINREDUCTION(h, e, ℓ, A, g, r)**

**Input:** $h$ separable, monic, of degree $d$ in $\mathbb{K}[x]$, $e$ in $\mathbb{N}$, $\ell$ in $\mathbb{N}_{>0}$, $A$ in $\mathbb{K}[\theta, z]_{<(d, \ell p^e)}$, $g$ in $\mathbb{K}[y]$, $r$ in $\mathbb{K}^{d+\lceil d^\eta \rceil}$

**Output:** $G \in \mathbb{K}[\theta, y]_{<(d, \ell p^e)}$ such that $G(\theta, A) \equiv g(A) \mod (h(\theta), (z^{p^e} - \theta)^\ell)$, or FAIL

1: Write $A = \sum_{0 \le i < \ell p^e} A_i z^i$ \qquad $\triangleright A_i \in \mathbb{K}[\theta]_{<d}$

2: $\triangleright$ Compute $\alpha$ s.t. the characteristic polynomial of $\tilde{A}$ is $(y^{p^e} - \alpha)^\ell$ (see Lemma 9.7)
$\alpha \leftarrow \sum_{0 \le i < \ell p^e} A_i p^e \theta^i$; $\alpha \leftarrow \alpha \text{ rem } h$ \qquad $\triangleright \alpha \in \mathbb{K}[\theta]_{<d}$

3: Write $g = \sum_{0 \le i < p^e} g_i (y^{p^e}) y^i$ \qquad $\triangleright \deg(g_i) \le \deg(g)/p^e$

4: **for** $i = 0, \dots, p^e - 1$ **do**
$G_i \leftarrow \text{BIVARIATE REDUCTION}(h, \ell, \alpha, g_i, r)$ \qquad $\triangleright G_i \in \mathbb{K}[\theta, y]_{<(d, \ell)}$

5: $G \leftarrow \sum_{0 \le i < p^e} G_i(\theta, y^{p^e})y^i$ \qquad $\triangleright G \in \mathbb{K}[\theta, y]_{<(d, \ell p^e)}$

6: **return** $G$

of success: each of the $p^e$ calls to Algorithm **BIVARIATE REDUCTION** succeeds with probability at least $1 - 6(\ell + 1)d^2/\text{card}(S)$. $\square$

## 9.4 Composition modulo powers

We now consider $f = h(x^{p^e})^\ell$, with $h$ separable of degree $d$ and integers $e, \ell$, with $\ell$ positive and not divisible by $p$ (and $e = 0$ if $p = 0$); the degree of $f$ is $n = d\ell p^e$. Algorithm **MODULARCOMPOSITIONMODULOPOWER** computes $g(a)$ rem $f$, extending to $e \ne 0$ the approach of van der Hoeven and Lecerf [30] outlined in Section 9.3.2.

We first compute $A(\theta, z) = a(z)$ rem $\langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$; this is done using the general untangling operation of Section 9.3.3. The reduction of the degree of $g$ is done by Algorithm **MAINREDUCTION**, giving $G$ in $\mathbb{K}[\theta, y]_{<(d, \ell p^e)}$, such that $G(\theta, A) \equiv g(A) \mod \langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$; the construction of $A$ then implies $G(\theta, A) \equiv g(a(z)) \mod \langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$. The quantity $B = G(\theta, A)$ rem $\langle h(\theta), (z^{p^e} - \theta)^\ell \rangle$ is obtained by Algorithm **COMPOSITIONMODULOINSEPARABLE-PRODUCTOFFIELDS**. We finally apply the general tangling procedure of Section 9.3.3 to $B$; since tangling is a $\mathbb{K}$-algebra isomorphism, the outcome is $b = g(a)$ rem $h(x^{p^e})^\ell$.

**Algorithm 9.4 MODULARCOMPOSITIONMODULOPOWER(h, e, ℓ, a, g, r)**

**Input:** $h$ separable, monic, of degree $d$ in $\mathbb{K}[x]$, $e$ in $\mathbb{N}$, $\ell$ in $\mathbb{N}_{>0}$, such that $f = h(x^{p^e})^\ell$ has degree $n = d\ell p^e$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[y]$, $r$ in $\mathbb{K}^{p+\lceil p^\eta \rceil}$ where $\rho = \max(d, n/d)$

**Output:** $b = g(a)$ rem $f$ or FAIL

1: $\triangleright$ Conversion of $a \in \mathbb{K}[x]$ to a bivariate polynomial (Proposition 9.6)
$A \leftarrow \text{UNTANGLING-GENERAL}(h, e, \ell, a)$ \qquad $\triangleright A \in \mathbb{K}[\theta, z]_{<(d, \ell p^e)}$

2: $\triangleright$ Reduction of $g$ modulo the characteristic polynomial of $\tilde{A}$ (Proposition 9.8)
$G \leftarrow \text{MAINREDUCTION}(h, e, \ell, A, g, (r_k)_{0 \le k < d + \lceil d^\eta \rceil})$ \qquad $\triangleright G \in \mathbb{K}[\theta, y]_{<(d, \ell p^e)}$

if $G = \text{FAIL}$ then return \text{FAIL}

3: $\triangleright$ Modular composition, $B = G(\theta, A)$ rem $\langle h(\theta), (z^{p^e} - \theta)^\ell \rangle \in \mathbb{K}[\theta, z]_{<(d, \ell p^e)}$ or FAIL
$B \leftarrow \text{COMPOSITIONMODULOINSEPARABLE-PRODUCTOFFIELDS}(h, \theta, e, \ell, A, G, (r_k)_{0 \le k < \frac{n}{d} + \lceil \frac{n}{d} \rceil \eta})$
if $B = \text{FAIL}$ then return \text{FAIL}

4: $\triangleright$ Recovery of $b$ over $\mathbb{K}$ (Proposition 9.6)
$b \leftarrow \text{TANGLING-GENERAL}(h, e, \ell, B)$

5: **return** $b$
**PROPOSITION 9.9.** For a field $\mathbb{K}$ of characteristic $p$, given $h$ separable, monic and of degree $d$ in $\mathbb{K}[x]$, integers $e$ in $\mathbb{N}$ and $\ell$ in $\mathbb{N}_{>0}$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[y]$, $r$ in $\mathbb{K}^{\rho+[\rho^n]}$, with $n = d\ell p^e$, $\rho = \max(d, \ell p^e)$ and $\eta$ from Eq. (3), Algorithm *MODULARCOMPOSITIONMODULOPOWER* uses $\tilde{O}(\deg(g) + n^\kappa)$ operations in $\mathbb{K}$, with $\kappa < 1.43$ as in Eq. (1), and returns $g(a)$ rem $h(x^{p^e})^\ell$ or *FAIL*.

If the entries of $r$ are chosen uniformly and independently from a finite subset $S$ of $\mathbb{K}$, then the algorithm returns $g(a)$ rem $h(x^{p^e})^\ell$ with probability at least $1 - (2n^4 + 12n^2)/\text{card}(S)$.

**PROOF.** That the output of the algorithm is $g(a)$ rem $h(x^{p^e})^\ell$ or *FAIL* follows from the previous discussion. By Proposition 9.6, with $n = d\ell p^e$, the first and last step both take $\tilde{O}(n)$ operations in $\mathbb{K}$. Proposition 9.8 shows that Step 2 takes $\tilde{O}(\deg(g) + n^\kappa)$ operations in $\mathbb{K}$. Finally, Proposition 9.2 shows that Step 3 takes $\tilde{O}(d(\ell p^e)^\kappa) = \tilde{O}(d(n/d)^\kappa)$ operations in $\mathbb{K}$, so the runtime estimate is proved.

The steps that may output *FAIL* are the computation of $G$ at Step 2 and that of $B$ at Step 3. By Proposition 9.8, the former happens with probability at most $6(\ell + 1)d^2p^e/\text{card}(S) \le 12n^2/\text{card}(S)$; by Proposition 9.2, the latter happens with probability at most $2d(\ell p^e)^4/\text{card}(S) \le 2n^4/\text{card}(S)$. $\square$

## 9.5 Main algorithm and its analysis

We can now give Algorithm *MODULARCOMPOSITION* performing modular composition with general polynomials, and prove Theorem 1.1.

The separable decomposition $f_1 \cdots f_s$ of $f$ allows us to reduce the problem to compositions modulo the $f_i$'s, which are powers of polynomials as in Section 9.4. The polynomials $a$ and $g$ are first reduced so that the compositions modulo the $f_i$'s are called with inputs of appropriate degrees, then the result $b = g(a)$ rem $f$ is recovered using Chinese remaindering. The number of random elements in $\mathbb{K}$ we use is an *a priori* bound that can be refined if the separable decomposition of $f$ is known.

**PROOF OF THEOREM 1.1.** First we prove correctness. Suppose that none of the subroutines returns *FAIL*; we show that the output is $g(a)$ mod $f$.

Using the same notation for $p^e$-th powering as in Lemma 9.7, at the $i$-th pass in the loop at Step 4, the polynomial $\mu_i$ satisfies $\mu_i(\alpha_i) \equiv 0 \bmod h_i$, with $\alpha_i = a_i^\sigma$ (that is, the coefficients of $\alpha_i$ are the $p^{e_i}$-th powers of those of $a_i$). Raising this equality to the power $\ell_i$ gives $\mu_i^{\ell_i}(\alpha_i) \equiv 0 \bmod h_i^{\ell_i}$. Evaluation at $y^{p^{e_i}}$ using the facts that $\alpha_i(y^{p^{e_i}}) = a_i^{p^{e_i}}$ and $\chi_i = \mu_i(y^{p^{e_i}})^{\ell_i}$ finally gives $\chi_i(a_i) \equiv 0 \bmod f_i$. The degree bound $\deg(\mu_i) \le 4d_i$ follows from the specifications of Algorithm *ANNIHILATINGPOLYNOMIAL*, and the degree bound for $\chi_i$ follows.

In the second for-loop at Step 6, $b_i$ satisfies $b_i \equiv g_i(a_i) \bmod h_i(x^{p^{e_i}})^{\ell_i} \equiv g_i(a_i) \bmod f_i$. Since $g_i = g$ rem $\chi_i$, and $\chi_i$ cancels $a_i$ modulo $f_i$, $b_i$ is also equal to $g(a_i)$ rem $f_i$, and thus to $g(a)$ rem $f_i$. It follows that the return value, obtained by Chinese remaindering, is indeed $g(a)$ rem $f$.

Next, we bound the overall cost. The call to *SEPARABLEDECOMPOSITION*($f$) takes $\tilde{O}(n)$ operations in $\mathbb{K}$ [53, Prop. 5]. Using repeated squaring, the polynomials $f_1, \dots, f_s$ can be computed in quadratic time as well, and the same holds for the remainders $a_1, \dots, a_s$.

Consider a fixed index $i$ in the loop at Step 4, and denote $d_i \ell_i p^{e_i}$ by $n_i$. Working coefficient-wise, computing $\alpha_i = a_i^\sigma$ takes $\tilde{O}(n_i)$ operations since $e_i = O(\log(n_i))$, and reducing it modulo $h_i$ has the same complexity bound. By Proposition 8.1, Algorithm *ANNIHILATINGPOLYNOMIAL* uses $\tilde{O}(d_i^\kappa)$ operations in $\mathbb{K}$. If it does not fail, $\chi_i$ is then deduced in $\tilde{O}(n_i)$ operations again, hence the cost of the loop is $\tilde{O}(n^\kappa)$.

When Step 5 is reached, since all $\chi_i$'s have respective degrees at most $4n_i$, fast multiple remaindering gives the polynomials $g_i$ in $\tilde{O}(n)$ operations, with $\deg(g_i) < 4n_i$. Then, by Proposition 9.9, each call to Algorithm *MODULARCOMPOSITIONMODULOPOWER* uses $\tilde{O}(n_i^\kappa)$ operations in $\mathbb{K}$, so their
<caption>**Algorithm 9.5 MODULARCOMPOSITION(f, a, g, r)**</caption>

Input: $f$ of degree $n$ in $\mathbb{K}[x]$, $a$ in $\mathbb{K}[x]_{<n}$, $g$ in $\mathbb{K}[y]_{<n}$, $r \in \mathbb{K}^{n+\lceil n^\eta \rceil}$

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

1: $\triangleright$ Decomposition of $f$ [53, Algo. 3]
   ($h_1, e_1, \ell_1), \dots, (h_s, e_s, \ell_s) \leftarrow \text{SEPARABLEDECOMPOSITION}(f)$ $\triangleright$ $h_i$ monic of degree $d_i$ in $\mathbb{K}[x]$

2: ($f_1, \dots, f_s) \leftarrow (h_1(x^{p^{e_1}})^{\ell_1}, \dots, h_s(x^{p^{e_s}})^{\ell_s})$ $\triangleright$ $f_i$ of degree $n_i = d_i \ell_i p^{e_i}$ in $\mathbb{K}[x]$

3: $\triangleright$ Degree reduction, $\deg(a_i) < n_i$
($a_1, \dots, a_s) \leftarrow (a \text{ rem } f_1, \dots, a \text{ rem } f_s)$

4: $\triangleright$ Annihilating polynomials of the $a_i$ modulo $f_i$

for $i = 1, \dots, s$ do

Write $a_i = \sum_{0 \le k < n_i} a_{i,k} x^k$
$\alpha_i \leftarrow \sum_{0 \le k < n_i} a_{i,k} p^{e_i} x^k$; $\alpha_i \leftarrow \alpha_i \text{ rem } h_i$
$\mu_i \leftarrow \text{ANNIHILATINGPOLYNOMIAL}(h_i, \alpha_i, (r_k)_{0 \le k < d_i + \lceil d_i^\eta \rceil})$ $\triangleright$ $\mu_i(\alpha_i) \equiv 0 \bmod h_i$, $\deg(\mu_i) \le 4d_i$

if $\mu_i = \text{FAIL}$ then return \text{FAIL}

$\chi_i \leftarrow \mu_i(y^{p^{e_i}})^{\ell_i}$ $\triangleright$ $\chi_i(a_i) \equiv 0 \bmod f_i$, $\deg(\chi_i) \le 4n_i$

5: $\triangleright$ Degree reduction, $\deg(g_i) < 4n_i$
($g_1, \dots, g_s) \leftarrow (g \text{ rem } \chi_1, \dots, g \text{ rem } \chi_s)$

6: $\triangleright$ Modular compositions, either $b_i \equiv g(a)$ mod $f_i$ or FAIL

for $i = 1, \dots, s$ do

$\rho_i \leftarrow \max(d_i, n_i/d_i)$
$b_i \leftarrow \text{MODULARCOMPOSITIONMODULOPOWER}(h_i, e_i, \ell_i, a_i, g_i, (r_k)_{0 \le k < \rho_i + \lceil \rho_i^\eta \rceil})$

if $b_i = \text{FAIL}$ then return \text{FAIL}

7: return CHINESEREMAINDERING((b₁, ..., b₅), (f₁, ..., f₅))

total cost is $\tilde{O}(n^\kappa)$ again. Finally, the cost of the last step (if reached) is $\tilde{O}(n)$. Altogether, the cost is $\tilde{O}(n^\kappa)$, as claimed.

It remains to discuss the probability of failure. By Proposition 8.3, the $i$th call to Algorithm **ANNIHILATINGPOLYNOMIAL** fails with probability at most $6d_i^2/\text{card}(S)$; hence, the probability that we successfully exit the first for-loop is at least $1-6n^2/\text{card}(S)$. Then, by Proposition 9.9, the $i$th call to Algorithm **MODULARCOMPOSITIONMODULOPOWER** fails with probability at most $(2n_i^4+12n_i^2)/\text{card}(S)$, so the probability that we successfully exit the second for-loop is at least $1-(2n^4+12n^2)/\text{card}(S)$. Altogether this gives a failure probability of at most $(2n^4+18n^2)/\text{card}(S)$. $\square$

# 10 APPLICATIONS

We now list several variants of the modular composition problem and related ones and sketch how the algorithms presented above can improve the best known complexity.

## 10.1 Annihilating polynomials

10.1.1 Annihilating polynomial. A by-product of Algorithm MODULARCOMPOSITION is a Las Vegas algorithm that takes $\tilde{O}(n^\kappa)$ ($\kappa$ from Theorem 1.1) arithmetic operations for computing an annihilating polynomial for $a$ of degree at most $4n$.

Indeed, with the notation of the algorithm, for all $1 \le i \le s$, since $\chi_i(a_i) \equiv 0 \bmod f_i$ we have $\chi_i(a) = r_i f_i$ for some $r_i \in \mathbb{K}[x]$. Hence $\prod_{i=1}^s \chi_i$ is an annihilating polynomial for $a$ modulo $f = \prod_{i=1}^s f_i$, whose degree is at most $4 \sum_{i=1}^s n_i = 4n$.
### 10.1.2 Minimal and characteristic polynomial.
In general, our knowledge of the minimal and characteristic polynomial depends on whether we have a certified basis of relations.

**PROPOSITION 10.1.** Let $R \in \mathbb{K}[y]_{\le 2d}^{m \times m}$ be the matrix produced by Algorithm CANDIDATEBASIS. If $R$ is a basis of $M_m^{(a,f)}$, then the first $m$ invariant factors of $yI_n - M_a$, hence in particular the minimal polynomial $\mu_a \in \mathbb{K}[y]$ of a modulo $f$, can be computed in $\tilde{O}(m^\omega d)$ operations in $\mathbb{K}$. If furthermore CERT is returned (implying that $R$ is a basis of $M_m^{(a,f)}$), then the product of these invariant factors gives the characteristic polynomial $\chi_a \in \mathbb{K}[y]$ of a modulo $f$.

**PROOF.** If $R$ is a basis of $M_m$, Proposition 4.1 shows that the Hermite normal form of $R$ is a triangular basis of $M_m$ whose diagonal entries are the first invariant factors $\sigma_1, \dots, \sigma_m$ of $yI_n - M_a$; in particular $\mu_a = \sigma_1$. If CERT is returned, then $R$ is a basis of $M_m$ and $v_m = n$ (Proposition 5.6). Hence $\deg \det(R) = n$ and all the invariant factors are known; the characteristic polynomial is their product. The Hermite normal form of $R$ can be computed in $\tilde{O}(m^\omega d)$ operations [49, Thm. 1.2]. $\square$

One case of certification of the minimal polynomial is when CERT is returned by Algorithm CANDIDATEBASIS, which occurs in particular for any $f$ in $\mathbb{K}[x]$ with $f(0) \ne 0$ and a generic $a$ in $\mathbb{K}[x]_{<n}$ (see Section 7.3.2). Using Proposition 5.6 and a shift as in Remark 3.8, this establishes the complexity bound $\tilde{O}(n^\kappa)$ for computing a basis of relations and the minimal polynomial in the case of a generic $a \in \mathbb{K}[x]_{<n}$.

Under the assumptions of Proposition 8.3 with the additional hypothesis $v_m^{(a,f)} = n$ for $m = [n^\eta]$, a call to Algorithm CANDIDATEBASIS instead of a call to Algorithm MATRIXOFRELATIONS in Algorithm MODULARCOMPOSITIONBASECASE, leads to a certified basis of relations of $M_m^{(\alpha, \mu_\gamma)}$ with good probability (use Proposition 5.6 instead of Proposition 5.8 in the proof of Proposition 8.3). From Proposition 10.1, this also allows one to compute and certify the minimal and characteristic polynomials in time $\tilde{O}(n^\kappa)$ when $f$ is separable and $v_{[n^\eta]}^{(a,f)} = n$.

The latter can be extended to the case $f$ irreducible and separable since then the minimal polynomial $\mu_a$ must be irreducible as well, and therefore $yI_n - M_a$ has $r$ nontrivial invariant factors all equal to $\mu_a$. If for $m = [n^\eta]$ the minimal polynomial satisfies $\delta = \deg(\mu_a) \ge n/m$, then $r \le m$ and $v_m^{(a,f)} = n$, hence the above certification when $f$ is separable leads to the minimal polynomial. The low degree case $\delta < n/m$ can be treated directly using Lemma 3.2, allowing to compute $\mu_a$ in time $\tilde{O}(n\delta^{(\omega_2/2)^{-1}})$, which is $\tilde{O}(n^\kappa)$ since $\delta < [n^{1-\eta}]$.

However, a matrix $R$ returned by Algorithm CANDIDATEBASIS might not be a basis of $M_m^{(\alpha, \mu_\gamma)}$ without an efficient certification of this property, Proposition 10.1 only gives a minimal polynomial algorithm of the Monte Carlo kind. Proceeding as done above, with a call to Algorithm CANDIDATEBASIS instead of a call to Algorithm MATRIXOFRELATIONS in Algorithm MODULARCOMPOSITIONBASECASE, a Monte Carlo minimal polynomial algorithm in $\tilde{O}(n^\kappa)$ can be derived under the assumptions of Proposition 8.3.

## 10.2 Power series reversion and power series equations

In this subsection, the characteristic of $\mathbb{K}$ is 0.

For a given $a \in \mathbb{K}[x]$ with $a(0) = 0$ and $a'(0) \ne 0$, power series reversion (or functional inversion) asks for a power series $g \in \mathbb{K}[[x]]$ such that

$$a(g) = g(a) \equiv x \mod x^n.$$

By Newton's iteration, a composition algorithm in $\tilde{O}(n^c)$ operations for some $c > 1$ induces a reversion algorithm in $\tilde{O}(n^c)$ operations as well [15]. Thus, we get a Las Vegas algorithm for power
series reversion in $\tilde{O}(n^\kappa)$ operations in $\mathbb{K}$. Note that the converse reduction, from reversion to composition, also holds in this situation [15].

The approach for reversion extends partially to the resolution of a class of power series equations. The aim is to solve an equation

$$g(x, y) = b \mod x^n \qquad (36)$$

for $y \in \mathbb{K}[[x]]_{<n}$, when $g \in \mathbb{K}[[x]][y]$ satisfies $g(0, 0) = b(0)$ and its partial derivative with respect to $y$ is not 0 at $(0, 0)$.

By Proposition 8.7, Algorithm COMPOSITIONMODULOINSEPARABLE computes a composition $g(x, a)$ in $\tilde{O}(n^\kappa)$ operations for $g$ in $\mathbb{K}[x, y]_{<(m, n)}$ with $m = O(n^\eta)$ and $\eta \approx 0.313$ from Eq. (3). Together with Newton's iteration, this gives a Las Vegas algorithm solving Eq. (36) in $\tilde{O}(n^\kappa)$ operations for $g \in \mathbb{K}[x, y]_{<(n^\eta, n)}$. Reversion is the special case with $b = x$ and $\deg_x(g) = 0$.

*Note.* It is known that the complexity of composition of power series (in terms of nonscalar operations) is essentially that of computing the coefficient of $x^{n-1}$ of $g(a)$ [64]. By contrast, computing the coefficient of $x^{n-1}$ in the reverse of $a$ costs only $\tilde{O}(n)$ arithmetic operations [15].

## 10.3 Bivariate composition

In this subsection, the characteristic of $\mathbb{K}$ is 0.

Brent and Kung gave an algorithm that computes

$$g(a, b) \text{ rem } x^n$$

for $g \in \mathbb{K}[x, y]_{< (n, n)}$ and truncated power series $a, b \in \mathbb{K}[[x]]$ in only $\tilde{O}(n^2)$ operations [14]. This is quasi-optimal, since the number of coefficients of $g$ is $\Theta(n^2)$ in general. In the simple situation where $a(0) = 0$ and $a'(0) = 1$, the algorithm is as follows:

(1) by power series reversion, compute $s(x)$ such that $a(s) = s(a) \equiv x \mod x^n$;

(2) by univariate composition, compute $c = b(s) \text{ rem } x^n$;

(3) by uni-bivariate composition, compute $d = g(x, c) \text{ rem } x^n$;

(4) by univariate composition, compute $d(a) \text{ rem } x^n$.

The complexity is dominated by the uni-bivariate composition in Step (3), which can be performed by Horner evaluation in $\tilde{O}(n^2)$ operations.

We obtain a Las Vegas algorithm with a complexity reduced to $\tilde{O}(n^\kappa)$ when $g \in \mathbb{K}[x, y]_{< (n^\eta, n)}$, where the uni-bivariate composition is done in $\tilde{O}(n^\kappa)$ as discussed in the case of power series equations, and all the other steps are univariate compositions that are also performed in $\tilde{O}(n^\kappa)$ by our algorithm.

This method extends to the computation of

$$g(a, b) \text{ rem } f$$

with $f$ of degree $n$ in $\mathbb{K}[x]$, and $a, b$ in $\mathbb{K}[x]_{<n}$. The algorithm becomes

(1) compute an annihilating polynomial $\chi$ of $a$ modulo $f$;

(2) by inverse modular composition, compute $c$ such that $c(a) \equiv b \mod f$;

(3) by uni-bivariate composition, compute $d = g(x, c) \text{ rem } \chi$;

(4) by univariate composition, compute $d(a) \text{ rem } f$.

At least for generic $a$, this is again a Las Vegas algorithm in $\tilde{O}(n^\kappa)$ operations when $g \in \mathbb{K}[x, y]_{< (n^\eta, n)}$.

## REFERENCES

[1] S. Abelard, A. Couvreur, and G. Lecerf. 2020. Sub-quadratic time for Riemann-Roch spaces: case of smooth divisors over nodal plane projective curves. In Proc. ISSAC. ACM Press, 14–21. https://doi.org/10.1145/3373207.3404053
[2] S. Abelard, A. Couvreur, and G. Lecerf. 2021. *Efficient computation of Riemann-Roch spaces for plane curves with ordinary singularities*. HAL Report hal-03110135. https://hal.archives-ouvertes.fr/hal-03110135

[3] J. Alman and V. Vassilevska Williams. 2021. A refined laser method and faster matrix multiplication. In Proc. SODA. SIAM, 522–539. https://doi.org/10.1137/1.9781611976465.32

[4] B. Beckermann and G. Labahn. 1994. A uniform approach for the fast computation of matrix-type Padé approximants. SIAM J. Matrix Analysis and Applications 15, 3 (1994), 804–823. https://doi.org/10.1137/S0895479892230031

[5] B. Beckermann, G. Labahn, and G. Villard. 1999. Shifted Normal Forms of Polynomial Matrices. In Proc. ISSAC. ACM Press, 189–196. https://doi.org/10.1145/309831.309929

[6] D. J. Bernstein. 1998. Composing power series over a finite ring in essentially linear time. J. Symb. Comput. 26 (1998), 339–341. https://doi.org/10.1006/jsco.1998.0216

[7] D. Bini and V. Y. Pan. 1994. *Polynomial and matrix computations*. Birkhäuser. https://doi.org/10.1007/978-1-4612-0265-3

[8] A. Bostan, F. Chyzak, M. Giusti, R. Lebreton, G. Lecerf, B. Salvy, and É. Schost. 2017. *Algorithmes efficaces en calcul formel*. In French. Edited by the authors. https://hal.archives-ouvertes.fr/AECF/

[9] A. Bostan, F. Chyzak, F. Ollivier, B. Salvy, É. Schost, and A. Sedoglavic. 2007. Fast computation of power series solutions of systems of differential equations. In Proc. SODA. SIAM, 1012–1021. https://dl.acm.org/doi/10.5555/1283383.1283492

[10] A. Bostan, P. Flajolet, B. Salvy, and É. Schost. 2006. Fast computation of special resultants. J. Symb. Comput. 41, 1 (2006), 1–29. https://doi.org/10.1016/j.jsc.2005.07.001

[11] A. Bostan, G. Lecerf, and É. Schost. 2003. Tellegen’s principle into practice. In Proc. ISSAC. ACM Press, 37–44. https://doi.org/10.1145/860854.860870

[12] A. Bostan, B. Salvy, and É. Schost. 2008. Power series composition and change of basis. In Proc. ISSAC. ACM Press, 269–276. https://doi.org/10.1145/1390768.1390806

[13] R. P. Brent, F. G. Gustavson, and D. Y. Y. Yun. 1980. Fast solution of Toeplitz systems of equations and computation of Padé approximants. J. Algorithms 1, 3 (1980), 259–295. https://doi.org/10.1016/0196-6774(80)90013-9

[14] R. P. Brent and H. T. Kung. 1977. Fast algorithms for composition and reversion of multivariate power series. In Proc. Conference on Theoretical Computer Science (Waterloo ON, August 15–17, 1977). University of Waterloo, 149–158.

[15] R. P. Brent and H. T. Kung. 1978. Fast algorithms for manipulating formal power series. J. ACM 25, 4 (1978), 581–595. https://doi.org/10.1145/322092.322099

[16] M. W. Buck, R. A. Coley, and D. P. Robbins. 1992. A generalized Vandermonde determinant. J. Algebraic Combin. 1, 2 (1992), 105–109. https://doi.org/10.1023/A:1022468019197

[17] M. Bürgisser, P. and Clausen and M. A. Shokrollahi. 1997. *Algebraic complexity theory*. Grundlehren der mathematischen Wissenschaften, Vol. 315. Springer. https://doi.org/10.1007/978-3-662-03338-8

[18] D. Coppersmith. 1994. Solving Homogeneous Linear Equations over GF(2) via Block Wiedemann Algorithm. Math. Comput. 62, 205 (1994), 333–350. https://doi.org/10.2307/2153413

[19] J. Della Dora, C. Discrescenzo, and D. Duval. 1985. About a new method for computing in algebraic number fields. In EUROCAL'85 (LNCS), Vol. 204. Springer, 289–290. https://doi.org/10.1007/3-540-15984-3_279

[20] W. Eberly, M. Giesbrecht, P. Giorgi, A. Storjohann, and G. Villard. 2007. Faster inversion and other black box matrix computation using efficient block projections. In Proc. ISSAC. ACM Press, 143–150. https://doi.org/10.1145/1277548.1277569

[21] M. Gasca and J. J. Martínez. 1987. On the computation of multivariate confluent Vandermonde determinants and its applications. In Proc. Mathematics of Surfaces II. Vol. 11. Oxford Univ. Press, 101–114.

[22] J. von zur Gathen and J. Gerhard. 1999. *Modern computer algebra*. Third edition 2013. Cambridge University Press. https://doi.org/10.1017/CBO9781139856065

[23] J. von zur Gathen and V. Shoup. 1992. Computing Frobenius maps and factoring polynomials. Comput. Complex. 2 (1992), 187–224. https://doi.org/10.1007/BF01272074

[24] J. Gerhard. 2004. *Modular algorithms in symbolic summation and symbolic integration*. Springer. https://doi.org/10.1007/b104035

[25] P. Gianni and Trager B. 1996. Square-free algorithms in positive characteristic. Appl. Algebra Eng. Commun. Comput. 7, 1 (1996), 1–14. https://doi.org/10.1007/BF01613611

[26] M. Giesbrecht, A. Jamshidpey, and É. Schost. 2021. Subquadratic-time algorithms for normal bases. Comput. Complex. 30, 5 (2021). https://doi.org/10.1007/s00037-020-00204-9

[27] P. Giorgi, C. Jeannerod, and G. Villard. 2003. On the complexity of polynomial matrix computations. In Proc. ISSAC. ACM Press, 135–142. https://doi.org/10.1145/860854.860889

[28] J. van der Hoeven. 2002. Relax, but Don't be Too Lazy. J. Symb. Comput. 34, 6 (2002), 479–542. https://doi.org/10.1006/jsco.2002.0562

[29] J. van der Hoeven and R. Larrieu. 2019. Fast Gröbner basis computation and polynomial reduction for generic bivariate ideals. Appl. Algebra Eng. Comm. 30, 6 (2019), 509–539. https://doi.org/10.1007/s00200-019-00389-9
[30] J. van der Hoeven and G. Lecerf. 2017. Composition modulo powers of polynomials. In Proc. *ISSAC*. ACM Press, 445-452. https://doi.org/10.1145/3087604.3087634

[31] J. van der Hoeven and G. Lecerf. 2018. Modular composition via factorization. *J. Complexity* 48, 36-68. https://doi.org/10.1016/j.jco.2018.05.002

[32] J. van der Hoeven and G. Lecerf. 2019. Accelerated tower arithmetic. *J. Complexity* 55 (2019). https://doi.org/10.1016/j.jco.2019.03.002

[33] J. van der Hoeven and G. Lecerf. 2020. Directed evaluation. *J. Complexity* 60 (2020). https://doi.org/10.1016/j.jco.2020.101498

[34] J. van der Hoeven and G. Lecerf. 2021. Amortized bivariate multi-point evaluation. In Proc. *ISSAC*. ACM Press, 179-185. https://doi.org/10.1145/3452143.3465531

[35] J. van der Hoeven and G. Lecerf. 2021. Fast amortized multi-point evaluation. *J. Complexity* (2021), 101574. https://doi.org/10.1016/j.jco.2021.101574

[36] J. van der Hoeven and G. Lecerf. 2021. Fast computation of generic bivariate resultants. *J. Complexity* 62 (2021). https://doi.org/10.1016/j.jco.2020.101499

[37] X. Huang and V. Y. Pan. 1998. Fast rectangular matrix multiplication and applications. *J. Complexity* 14 (1998), 257-299. https://doi.org/10.1006/jcom.1998.0476

[38] C.-P. Jeannerod, V. Neiger, and G. Villard. 2020. Fast computation of approximant bases in canonical form. *J. Symb. Comput.* 98 (2020), 192-224. https://doi.org/10.1016/j.jsc.2019.07.011

[39] T. Kailath. 1980. *Linear Systems*. Prentice-Hall.

[40] E. Kaltofen. 1992. On computing determinants without divisions. In Proc. *ISSAC*. ACM Press, 342-349. https://doi.org/10.1145/143242.143350

[41] E. Kaltofen. 2000. Challenges of symbolic computation: my favorite open problems. *J. Symb. Comput.* 29, 6 (2000), 891-919. https://doi.org/10.1006/jSCO.2000.0370

[42] E. Kaltofen and V. Y. Pan. 1991. Processor efficient parallel solution of linear systems over an abstract field. In Proc. SPAA. ACM, 180-191. https://doi.org/10.1145/113379.113396

[43] E. Kaltofen and D. Saunders. 1991. On Wiedemann's method of solving sparse linear systems. In *AAECC-9 (LNCS)*, Vol. 539. Springer Verlag, 29-38. https://doi.org/10.1007/3-540-54522-0_93

[44] E. Kaltofen and V. Shoup. 1997. Fast polynomial factorization over high algebraic extensions of finite fields. In Proc. *ISSAC*. ACM Press, 184-188. https://doi.org/10.1145/258726.258777

[45] E. Kaltofen and V. Shoup. 1998. Subquadratic-time factoring of polynomials over finite fields. *Math. Comp.* 67, 233 (1998), 1179-1197. https://doi.org/10.1090/S0025-5718-98-00944-2

[46] E. Kaltofen and G. Villard. 2005. On the complexity of computing determinants. *Comput. Complex.* 13, 3 (2005), 91-130. https://doi.org/10.1007/s00037-004-0185-3

[47] E. Kaltofen and G. Yuhasz. 2013. On the matrix Berlekamp-Massey algorithm. *ACM Trans. Algorithms* 9, 4 (2013), 33:1-33:24. https://doi.org/10.1145/2500122

[48] K. S. Kedlaya and C. Umans. 2011. Fast polynomial factorization and modular composition. *SIAM J. on Computing* 40, 6 (2011), 1767-1802. https://doi.org/10.1137/08073408X

[49] G. Labahn, V. Neiger, and W. Zhou. 2017. Fast, deterministic computation of the Hermite normal form and determinant of a polynomial matrix. *J. Complexity* 42 (2017), 44-71. https://doi.org/10.1016/j.jco.2017.03.003

[50] D. Lazard. 1985. Ideal bases and primary decomposition: case of two variables. *J. Symb. Comput.* 1, 3 (1985), 261-270. https://doi.org/10.1007/s00747-7171(85)80035-3

[51] F. Le Gall. 2014. Powers of tensors and fast matrix multiplication. In Proc. *ISSAC*. ACM Press, 296-303. https://doi.org/10.1145/2608628.2608664

[52] F. Le Gall and F. Urrutia. 2018. Improved rectangular matrix multiplication using powers of the Coppersmith-Winograd tensor. In Proc. *SODA*. SIAM, 1029-1046. https://doi.org/10.1137/1.9781611975031.67

[53] G. Lecerf. 2008. Fast separable factorization and applications. *Appl. Algebra Eng. Commun. Comput.* 19, 2 (2008), 135-160. https://doi.org/10.1007/s00200-008-0062-4

[54] J. D. Lipson. 1976. Newton's method: a great algebraic algorithm. In Proc. *SYMSAC*. ACM Press, 260-270. https://doi.org/10.1145/800205.806344

[55] T. Mulders and A. Storjohann. 2003. On lattice reduction for polynomial matrices. *J. Symb. Comput.* 35 (2003), 377-401. Issue 4. https://doi.org/10.1016/S0747-7171(02)00139-6

[56] V. Neiger. 2016. *Bases of relations in one or several variables: fast algorithms and applications*. Ph.D. Dissertation. École Normale Supérieure de Lyon. https://tel.archives-ouvertes.fr/tel-01431413/

[57] V. Neiger, J. Rosenkilde, and G. Solomatov. 2020. Generic bivariate multi-point evaluation, interpolation and modular composition with precomputation. In Proc. *ISSAC*. ACM Press, 388-395. https://doi.org/10.1145/3373207.3404032

[58] M. Newman. 1972. *Integral Matrices*. Academic Press. First edition.
[59] M. Nüsken and M. Ziegler. 2004. Fast multipoint evaluation of bivariate polynomials. In *Algorithms – ESA 2004*. Springer, Berlin, Heidelberg, 544–555. https://doi.org/10.1007/978-3-540-30140-0_49

[60] M. Paterson and L. J. Stockmeyer. 1973. On the number of nonscalar multiplications necessary to evaluate polynomials. *SIAM J. Comput.* 2, 1 (1973), 60–66. https://doi.org/10.1137/0202007

[61] A. Poteaux and É. Schost. 2013. Modular composition modulo triangular sets and applications. *Comput. Complex.* 22, 3 (2013), 463–516. https://doi.org/10.1007/s00037-013-0063-y

[62] A. Poteaux and É. Schost. 2013. On the complexity of computing with zero-dimensional triangular sets. *J. Symb. Comput.* 50 (2013), 110–138. https://doi.org/10.1016/j.jsc.2012.05.008

[63] D. Reischert. 1997. Asymptotically fast computation of subresultants. In *Proc. ISSAC*. ACM Press, 233–240. https://doi.org/10.1145/258726.258792

[64] P. Ritzmann. 1986. A fast numerical algorithm for the composition of power series with complex coefficients. *Theoret. Comput. Sci.* 44, 1 (1986), 1–16. https://doi.org/10.1016/0304-3975(86)90107-6

[65] V. Shoup. 1994. Fast construction of irreducible polynomials over finite fields. *J. Symb. Comput.* 17, 5 (1994), 371–391. https://doi.org/10.1006/jsco.1994.1025

[66] V. Shoup. 1995. A new polynomial factorization algorithm and its implementation. *J. Symb. Comput.* 20, 4 (1995), 363–397. https://doi.org/10.1006/jsco.1995.1055

[67] V. Shoup. 1999. Efficient computation of minimal polynomials in algebraic extensions of finite fields. In *Proc. ISSAC*. ACM Press, 53–58. https://doi.org/10.1145/309831.309859

[68] C. Umans. 2008. Fast polynomial factorization and modular composition in small characteristic. In *Proc. STOC*. ACM Press, 481–490. https://doi.org/10.1145/1374376.1374445

[69] M. Van Barel and A. Bultheel. 1992. A general module theoretic framework for vector M-Padé and matrix rational interpolation. *Numer. Algorithms* 3 (1992), 451–462. https://doi.org/10.1007/BF02141952

[70] G. Villard. 1997. *A study of Coppersmith's block Wiedemann algorithm using matrix polynomials*. RR 975 IM IMAG. http://perso.ens-lyon.fr/gilles.villard/BIBLIOGRAPHIE/PDF/rr0497.pdf

[71] G. Villard. 2018. On computing the resultant of generic bivariate polynomials. In *Proc. ISSAC*. ACM Press, 391–398. https://doi.org/10.1145/3208976.3209020

[72] D. Wiedemann. 1986. Solving sparse linear equations over finite fields. *IEEE Trans. Information Theory* 32, 1 (1986), 54–62. https://doi.org/10.1109/TIT.1986.1057137

[73] W. A. Wolovich. 1974. *Linear Multivariable Systems*. Applied Mathematical Sciences, Vol. 11. Springer-Verlag New-York. https://doi.org/10.1007/978-1-4612-6392-0

[74] W. Zhou and G. Labahn. 2012. Efficient algorithms for order basis computation. *J. Symb. Comput.* 47, 7 (2012), 793–819. https://doi.org/10.1016/j.jsc.2011.12.009

[75] W. Zhou, G. Labahn, and A. Storjohann. 2012. Computing minimal nullspace bases. In *Proc. ISSAC*. ACM Press, 366–373. https://doi.org/10.1145/2442829.2442881