Implementations of efficient univariate polynomial matrix algorithms and application to bivariate resultants

TL;DR

This paper explores the practical implementation of efficient univariate polynomial matrix algorithms, demonstrating improved performance in computing bivariate resultants, with implications for error-correcting codes and polynomial systems.

Contribution

It provides the first practical implementation of several fundamental polynomial matrix operations and applies these to enhance bivariate resultant computations.

Findings

Improved performance over existing methods for large parameters

Successful implementation of key polynomial matrix operations

Enhanced algorithms for bivariate resultants

Abstract

Complexity bounds for many problems on matrices with univariate polynomial entries have been improved in the last few years. Still, for most related algorithms, efficient implementations are not available, which leaves open the question of the practical impact of these algorithms, e.g. on applications such as decoding some error-correcting codes and solving polynomial systems or structured linear systems. In this paper, we discuss implementation aspects for most fundamental operations: multiplication, truncated inversion, approximants, interpolants, kernels, linear system solving, determinant, and basis reduction. We focus on prime fields with a word-size modulus, relying on Shoup's C++ library NTL. Combining these new tools to implement variants of Villard's algorithm for the resultant of generic bivariate polynomials (ISSAC 2018), we get better performance than the state of the art for large parameters.

Equations

$${\mathcal{I}}_{\mathbf{\alpha }}(\mathbf{E})=\{\mathbf{p}\in \mathbb{K}[x{]}^{1\times m}\mid \mathbf{p}({\alpha }_i){\mathbf{E}}_i=0\text{ for }1\le i\le d\}.$$

$$h_{i,j,k}:=\mathrm{coeff}(h_{i,j},x^k)=\mathrm{coeff}(-z^j{A}^{-k-1}\mathrm{mod}P,z^i).$$

Full text

Polynomial matrices, algorithms, implementation, resultant.

††copyright: rightsretained††journalyear: 2019††conference: International Symposium on Symbolic and AlgebraicComputation; July 15–18, 2019; Beijing, China††booktitle: International Symposium on Symbolic and Algebraic Computation(ISSAC ’19), July 15–18, 2019, Beijing, China††price: 15.00††doi: 10.1145/3326229.3326272††isbn: 978-1-4503-6084-5/19/07

1. Introduction

Hereafter, $\mathbb{K}$ is a field and $\mathbb{K}[x]$ is the algebra of univariate polynomials over $\mathbb{K}$. Recent years have witnessed a host of activity on fast algorithms for polynomial matrices and their applications:

• •

Minimal approximant bases (Giorgi et al., 2003; Zhou and Labahn, 2012) were used to compute kernel bases (Zhou et al., 2012), giving the first efficient deterministic algorithm for linear system solving over $\mathbb{K}[x]$.

• •

Basis reduction (Giorgi et al., 2003; Gupta et al., 2012) played a key role in accelerating the decoding of one-point Hermitian codes (Nielsen and Beelen, 2015) and in designing deterministic determinant and Hermite form algorithms (Labahn et al., 2017).

• •

Progress on minimal interpolant bases (Jeannerod et al., 2017, 2016) led to the best known complexity bound for list-decoding Reed-Solomon codes and folded Reed-Solomon codes (Jeannerod et al., 2017, Sec. 2.4 to 2.7).

• •

Coppersmith’s block Wiedemann algorithm and its extensions (Coppersmith, 1994; Kaltofen, 1995; Villard, 1997) were used in a variety of contexts, from integer factorization (Team, 2017) to polynomial system solving (Villard, 2018; Hyun et al., 2017b).

At the core of these improvements, one also finds techniques such as high-order lifting (Storjohann, 2003) and partial linearization (Storjohann, 2006),(Gupta et al., 2012, Sec. 6).

For many of these operations, no implementation of the latest algorithms is available and no experimental evidence has been given regarding their practical behavior. Our goal is to partly remedy this issue, by providing and discussing implementations for a core of fundamental algorithms such as multiplication, approximant and interpolant bases, etc., upon which one may implement higher level algorithms. As an illustration, we describe the performance of slightly modified versions of Villard’s recent breakthroughs on bivariate resultant and characteristic polynomial computation (Villard, 2018).

Our implementation is based on Shoup’s Number Theory Library (NTL) (Shoup, 2018), and is dedicated to polynomial matrix arithmetic over $\mathbb{K}=\mathbb{F}_{p}$ for a word-size prime $p$. Particular attention was paid to performance issues, so that our library compares favorably with previous work for those operations where comparisons were possible. Our code is available at https://github.com/vneiger/pml.

Overview.

Polynomial matrix algorithms rely on efficient arithmetic in $\mathbb{K}[x]$ and for matrices over $\mathbb{K}$; in Section 2, we review some related algorithms and their NTL implementations. Then, we describe our implementation of a key building block: multiplication.

Section 3 presents the next major part of our work, concerning algorithms for approximant bases (Beckermann and Labahn, 1994; Giorgi et al., 2003; Zhou and Labahn, 2012; Jeannerod et al., 2018) and interpolant bases (Van Barel and Bultheel, 1992; Beckermann and Labahn, 2000; Jeannerod et al., 2016, 2017). We focus on a version of interpolants which is less general than in these references but allows for a more efficient algorithm. In particular, we show that with this version, both interpolant and approximant bases can be used interchangeably in several contexts, with interpolants sometimes achieving better performance than approximants. In Section 4, we discuss algorithms for minimal kernel bases, linear system solving, determinant, and basis reduction. Finally, using these tools, we study the practical behavior of the bivariate resultant algorithm of (Villard, 2018) (Section 5).

Below, cost bounds are given in an algebraic complexity model, counting all operations in the base field at unit cost. While standard, this point of view fails to describe parts of the implementation (CRT-based algorithms, such as the 3-primes FFT, cannot be described in such a manner), but we believe that this is a minor issue.

Implementation choices.

NTL is a C++ library for polynomial and matrix arithmetic over rings such as $\mathbb{Z}$, $\mathbb{Z}/n\mathbb{Z}$, etc., and is often seen as a reference point for fast implementations in such contexts. Other libraries for these operations include for example FLINT (Hart et al., 2015) as well as FFLAS-FFPACK and LinBox (The FFLAS-FFPACK Group, 2019; The LinBox Group, 2018). Currently, our implementation relies solely on NTL; this choice was based on comparisons of performance for the functionalities we need.

In our implementation, the base field is a prime finite field $\mathbb{F}_{p}$; we rely on NTL’s lzz_p class. At the time of writing, on standard x86_64 platforms, NTL v11.3.1 uses unsigned long’s as its primary data type for lzz_p, supporting moduli up to 60 bits long.

For such fields, one can directly compare running times and cost bounds, since in the literature most polynomial matrix algorithms are analyzed in the algebraic complexity model. Besides, computations over $\mathbb{F}_{p}$ are at the core of a general approach consisting in solving problems over $\mathbb{Z}$ or $\mathbb{Q}$ by means of reduction modulo sufficiently many primes, which are chosen so as to satisfy several, partly conflicting, objectives. We may want them to support Fourier transforms of high orders. Linear algebra modulo each prime should be fast, so we may wish them to be small enough to support vectorized matrix arithmetic with SIMD instructions. On the other hand, using larger primes allows one to use fewer of them, and reduces the likelihood of unlucky choices in randomized algorithms.

As a result, while all NTL lzz_p moduli are supported, our implementation puts an emphasis on three families: small FFT primes that support AVX-based matrix multiplication (such primes have at most 23 bits); arbitrary size FFT primes (at most 60 bits); arbitrary moduli (at most 60 bits). Very small fields such as $\mathbb{F}_{2}$ or $\mathbb{F}_{3}$ are supported, but we did not make any specific optimization for them.

Experiments.

All runtimes below are in seconds and were measured on an Intel Core i7-4790 CPU with 32GB RAM, using the version 11.3.1 of NTL. Unless specified otherwise, timings are obtained modulo a random 60 bit prime. Runtimes were measured on a single thread; currently, most parts of our code do not explicitly exploit multi-threading. Tables below only show a few selected timings, with the best time(s) in bold; for more timings, see https://github.com/vneiger/pml/tree/master/benchmarks.

2. Basic polynomial and matrix arithmetic

We review basic algorithms for polynomials and matrices, and related complexity results that hold over an abstract field $\mathbb{K}$, and we describe how we implemented these operations. Hereafter, for $d\geq 0$, $\mathbb{K}[x]_{d}$ is the set of elements of $\mathbb{K}[x]$ of degree less than $d$.

2.1.   Polynomial multiplication.

Multiplication in $\mathbb{K}[x]$ and Fast Fourier Transform (FFT) are cornerstones of most algorithms in this paper. Let $\mathsf{M}:\mathbb{N}\to\mathbb{N}$ be a function such that polynomials of degree at most $d$ in $\mathbb{K}[x]$ can be multiplied in $\mathsf{M}(d)$ operations in $\mathbb{K}$. If $\mathbb{K}$ supports FFT, we can take $\mathsf{M}(d)\in O(d\log(d))$, and otherwise, $\mathsf{M}(d)\in O(d\log(d)\log\log(d))$ (Gathen and Gerhard, 2013, Chapter 8); as in this reference, we assume that $d\mapsto\mathsf{M}(d)/d$ is increasing. A useful variant of multiplication is the middle product (Hanrot et al., 2004; Bostan et al., 2003): for integers $c$ and $d$, and $F$ in $\mathbb{K}[x]_{c}$ and $G$ in $\mathbb{K}[x]_{c+d}$, ${\textsc{MiddleProduct}}(F,G,c,d)$ returns the slice of the product $FG$ with coefficients of degrees $c,\dots,c+d-1$; a common case is with $c=d$. The direct approach computes the whole product and extracts the slice. Yet, the transposition principle (Kaminski et al., 1988) yields a more efficient approach, saving a constant factor (roughly a factor $2$ when $c=d$, if FFT multiplication is used).

Polynomial matrix algorithms frequently use fast evaluation and interpolation at multiple points. In general, subproduct tree techniques (Gathen and Gerhard, 2013, Chapter 10) allow one to do evaluation and interpolation of polynomials in $\mathbb{K}[x]_{d}$ at $d$ points in $O(\mathsf{M}(d)\log(d))$ operations. For special sets of points, one can do better: if we know $\alpha$ in $\mathbb{K}$ of order at least $d$, then evaluation and interpolation at the geometric progression $(1,\alpha,\ldots,\alpha^{d-1})$ can both be done in time $O(\mathsf{M}(d))$ (Bostan and Schost, 2005).

In NTL, multiplication in $\mathbb{F}_{p}[x]$ uses either naive, Karatsuba, or FFT techniques, depending on $p$ and on the degree (NTL provides FFT primes with roots of unity of order $2^{25}$, and supports arbitrary user-chosen FFT primes). FFT multiplication uses the TFT algorithm of (Harvey, 2009) and Harvey’s improvements on arithmetic mod $p$ (Harvey, 2014). For primes $p$ that do not support Fourier transforms, multiplication is done by means of either 3-primes FFT techniques (Gathen and Gerhard, 2013, Chapter 8) or Schönhage and Strassen’s algorithm. We implemented middle products for naive, Karatsuba and FFT multiplication, closely following (Hanrot et al., 2004; Bostan et al., 2003), as well as evaluation/interpolation algorithms for general sets of points and for geometric progressions.

2.2.   Matrix multiplication.

Let $\omega$ be such that $n\times n$ matrices over any ring can be multiplied by a bilinear algorithm doing $O(n^{\omega})$ ring operations. The naive algorithm does exactly $n^{3}$ multiplications. First improvements due to Winograd and Waksman (Winograd, 1968; Waksman, 1970) reduced the number of operations to $n^{3}/2+O(n^{2})$ if $2$ is a unit. Strassen’s and Winograd’s recursive algorithms (Strassen, 1969; Winograd, 1971) have $\omega=\log_{2}(7)$; the best known bound is $\omega\leq 2.373$ (Coppersmith and Winograd, 1990; Le Gall, 2014). Note that, using blocking, rectangular matrices of sizes $(m\times n)$ and $(n\times p)$ can be multiplied in $O(mnp\min(m,n,p)^{\omega-3})$ ring operations. NTL implements its own arithmetic for matrices over $\mathbb{F}_{p}$ and chooses one of several implementations depending on the bitsize of $p$, the matrix dimensions, the available processor instructions, etc.

2.3.   Polynomial matrix multiplication.

In what follows, we write $\mathsf{MM}(n,d)$ for a function such that two $n\times n$ matrices of degree at most $d$ can be multiplied in $\mathsf{MM}(n,d)$ operations in $\mathbb{K}$; we make the assumption that $d\mapsto\mathsf{MM}(n,d)/d$ is increasing for all $n$.

From the definitions above we obtain $\mathsf{MM}(n,d)\in O(n^{\omega}\mathsf{M}(d))$, which is in $\mathchoice{\tilde{O}\left(n^{\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega}d)}$. Using evaluation/interpolation at $1,\alpha,\ldots,\alpha^{2d}$ or at roots of unity, one obtains the following bounds on $\mathsf{MM}(n,d)$:

• •

$O(n^{\omega}d+n^{2}\mathsf{M}(d))$ if an element $\alpha$ in $\mathbb{K}$ of order more than $2d$ is known (Bostan and Schost, 2005, Thm. 2.4).

• •

$O(n^{\omega}d+n^{2}d\log(d))$ if $\mathbb{K}$ supports FFT in degree $2d$.

We also mention a polynomial analogue of an integer matrix multiplication algorithm from (Doliskani et al., 2018) which uses evaluation/interpolation, done plainly via multiplication by (inverse) Vandermonde matrices. Then, the corresponding part of the cost (e.g. $O(n^{2}\mathsf{M}(d))$ for geometric progressions) is replaced by the cost of multiplying matrices over $\mathbb{K}$ in sizes roughly $(d\times d)$ by $(d\times n^{2})$; this is in $O(n^{2}d^{\omega-1})$ if $d\leq n^{2}$. For moderate values of $d$, where $\mathsf{M}(d)$ is not in the FFT regime, this allows us to leverage fast matrix multiplication over $\mathbb{K}$.

We implemented and compared various algorithms for matrix multiplication over $\mathbb{F}_{p}[x]$. For matrices of degree less than 5, we use dedicated routines based on Karatsuba’s and Montgomery’s formulas (Montgomery, 2005); for matrices of small size (up to 10, depending on $p$), we use Waksman’s algorithm. For other inputs, most of our efforts were spent on variants of the evaluation/interpolation scheme.

For FFT primes, we use evaluation/interpolation at roots of unity. For general primes, we use either evaluation/interpolation at geometric progressions (if such points exist in $\mathbb{F}_{p}$), or our adaptation of the algorithm of (Doliskani et al., 2018), or 3-primes multiplication (as for polynomials, we lift the product from $\mathbb{F}_{p}[x]$ to $\mathbb{Z}[x]$, where it is done modulo up to 3 FFT primes). No single variant outperformed or underperformed all others for all sizes and degrees, so thresholds were experimentally determined to switch between these options, with different values for small (less than 23 bits) and for large primes.

Middle product versions of these algorithms were implemented, and are used in approximant basis algorithms (Section 3.1) and Newton iteration (Section 4.3). Multiplier classes are available: they do precomputations on a matrix $\mathbf{{A}}$ to accelerate repeated multiplications by $\mathbf{{A}}$; they are used in Dixon’s algorithm (Section 4.3).

The table below shows timings for our multiplication and LinBox’ one, for random $m\times m$ matrices of degree $d$ and two choices of prime $p$. The global comparison showed running times that are either similar or in favor of our implementation.

20 bit FFT prime 60 bit prime

$m$ $d$ ours Linbox ratio ours Linbox ratio

8 131072 1.1198 1.5930 0.70 3.577 13.59 0.26

32 4096 0.4283 0.5092 0.84 2.000 5.330 0.38

128 1024 1.7292 2.1126 0.82 15.73 23.13 0.68

512 128 4.3533 4.3837 0.99 41.57 50.62 0.82

3. Approximant bases and interpolant bases

These bases are matrix generalizations of Padé approximation and play an important role in many higher-level algorithms. For $\mathbf{{F}}$ in $\mathbb{K}[x]^{m\times n}$ and $M$ non-constant in $\mathbb{K}[x]$, they are bases of the $\mathbb{K}[x]$-module $\mathcal{A}_{M}(\mathbf{{F}})$ of all $\boldsymbol{p}$ in $\mathbb{K}[x]^{1\times m}$ such that $\boldsymbol{p}\mathbf{{F}}=0\bmod M$. Specifically, approximant bases are for $M=x^{d}$ and interpolant bases for $M=\prod_{i}(x-\alpha_{i})$ for $d$ distinct points $\alpha_{1},\ldots,\alpha_{d}$ in $\mathbb{K}$. (Here, we do not consider more general cases from the literature, for example with several moduli $M_{1},\ldots,M_{n}$, one for each column of $\boldsymbol{p}\mathbf{{F}}$.)

Since $\mathcal{A}_{M}(\mathbf{{F}})$ is free of rank $m$, such a basis is represented row-wise by a nonsingular $\mathbf{{P}}$ in $\mathbb{K}[x]^{m\times m}$. The algorithms below return $\mathbf{{P}}$ in $\boldsymbol{s}$-ordered weak Popov form (also known as $\boldsymbol{s}$-quasi Popov form (Beckermann et al., 1999)), for a given shift $\boldsymbol{s}=(s_{1},\dots,s_{m})$ in $\mathbb{Z}^{m}$. Shifts allow us to set degree constraints on the sought basis $\mathbf{{P}}$, and they inherently occur in a general approach for finding bases of solutions to equations (approximants, interpolants, kernels, etc.). Approximant basis algorithms often require $\mathbf{{P}}$ to be in $\boldsymbol{s}$-reduced form (Van Barel and Bultheel, 1992); although the $\boldsymbol{s}$-ordered weak Popov form is stronger, obtaining it involves minor changes in these algorithms, without impact on performance according to our experiments. Recent literature shows that this stronger form reveals valuable information for further computations with $\mathbf{{P}}$ (Jeannerod et al., 2016; Jeannerod et al., 2018), in particular for finding $\boldsymbol{s}$-Popov bases (Beckermann et al., 1999).

From the shift $\boldsymbol{s}$, the $\boldsymbol{s}$-degree of $\boldsymbol{p}=[p_{i}]_{i}\in\mathbb{K}[x]^{1\times m}$ is defined as $\mathrm{rdeg}_{\boldsymbol{s}}(\boldsymbol{p})=\max_{1\leq i\leq m}(\deg(p_{i})+s_{i})$, which extends to matrices: $\mathrm{rdeg}_{\boldsymbol{s}}(\mathbf{{P}})$ is the list of $\boldsymbol{s}$-degrees of the rows of $\mathbf{{P}}$. Then, the $\boldsymbol{s}$-pivot of $\boldsymbol{p}$ is its rightmost entry $p_{i}$ such that $\mathrm{rdeg}_{\boldsymbol{s}}(\boldsymbol{p})=\deg(p_{i})+s_{i}$, and a nonsingular matrix $\mathbf{{P}}$ is in $\boldsymbol{s}$-ordered weak Popov form if the $\boldsymbol{s}$-pivots of its rows are located on the diagonal.

To simplify cost bounds below, we make use of the function $\mathsf{MM}^{\prime}(m,d)={\sum}_{i=0}^{\log_{2}(d)}2^{i}\mathsf{MM}(m,d/2^{i})\in O(\mathsf{MM}(m,d)\log(d))$.

3.1.   Approximant bases.

For $\mathbf{{F}}$ in $\mathbb{K}[x]^{m\times n}$ and $d$ in $\mathbb{Z}_{>0}$, an approximant basis for $(\mathbf{{F}},d)$ is a nonsingular $m\times m$ matrix whose rows form a basis of $\mathcal{A}_{x^{d}}(\mathbf{{F}})$. We implemented minor variants of the algorithms M-Basis (iterative, via matrix multiplication) and PM-Basis (divide and conquer, via polynomial matrix multiplication) from (Giorgi et al., 2003). The lowest-level function (M-Basis-1 with the signature in Algorithm 1), handles order $d=1$ in time $O(\mathrm{rank}(\mathbf{{F}})^{\omega-2}mn)$; here, working modulo $X$, the matrix $\mathbf{{F}}$ is over $\mathbb{K}$. Our implementation follows (Jeannerod et al., 2018, Algo. 1), which returns an $\boldsymbol{s}$-Popov basis, using only an additional row permutation compared to the algorithm in (Giorgi et al., 2003).

This form of the output of M-Basis-1 suffices to ensure that M-Basis and PM-Basis return bases in $\boldsymbol{s}$-ordered weak Popov form. Our implementation of M-Basis follows the original design (Giorgi et al., 2003) with $d$ iterations, each computing the residual $\mathbf{{R}}$ and updating $\mathbf{{P}}$ via multiplication by a basis $\mathbf{{Q}}$ obtained by M-Basis-1 on $\mathbf{{R}}$. We also follow (Giorgi et al., 2003) for PM-Basis, using a threshold $T$ such that M-Basis is called if $d\leq T$. Building PM-Basis directly upon M-Basis-1, i.e. choosing $T=1$, has the same cost bound but is slower in practice.

These algorithms use $O((m^{\omega}+m^{\omega-1}n)d^{2})$ and $O((1+\frac{n}{m})\mathsf{MM}^{\prime}(m,d))$, respectively (Giorgi et al., 2003). Some implementation details are discussed in Section 3.2. The next table compares timings for LinBox’ and our implementations of PM-Basis for a 20 bit FFT prime (LinBox’ implementation is not optimized for large primes and general primes).

[TABLE]

We also implemented (Jeannerod et al., 2018, Algo. 3) which returns $\boldsymbol{s}$-Popov bases and is about twice slower than PM-Basis; making this overhead negligible for some usual cases is future work. For completeness, we handle general approximants (with one modulus per column of $\mathbf{{F}}$) by an iterative approach from (Van Barel and Bultheel, 1992; Beckermann and Labahn, 1994); faster algorithms are more complex (Jeannerod et al., 2016, 2017; Jeannerod et al., 2018) and use partial linearization techniques.

These techniques from (Storjohann, 2006; Zhou and Labahn, 2012) yield cost bounds in $\mathchoice{\tilde{O}\left(m^{\omega-1}nd\right)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega-1}nd)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega-1}nd)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega-1}nd)}$, which is a $\Theta(\frac{m}{n})$ speedup compared to PM-Basis. Implementing them is work in progress. experimental code, which focuses for simplicity on $n=1$ and “generic” inputs for which the degrees in $\mathbf{{P}}$ can be predicted, revealed significant speedups:

[TABLE]

Approximant bases are often applied to solve block-Hankel systems (Labahn et al., 1990). In two specific settings, we have compared this approach to the one which uses structured matrix algorithms; we are not aware of previous comparisons of this kind. We obtain the following running times, using the NTL-based solver from (Hyun et al., 2017a).

[TABLE]

Setting 1: we call PM-Basis on $[\mathbf{{F}}^{\mathsf{T}}\;\;-\mathbf{{I}}_{m}]^{\mathsf{T}}$ at order $2d$ with shift $(0,\ldots,0)$, where $\mathbf{{F}}$ is an $m\times m$ matrix of degree $2d-1$, and we solve a system with $m\times m$ Hankel blocks of size $d\times d$ (the structured solver returns a random solution to the system). Our experiments show a clear advantage for approximant algorithms. The asymptotic costs being similar, the effects at play here are constant factor differences: approximant basis algorithms seem to be somewhat simpler and to better leverage the main building blocks (matrix arithmetic over $\mathbb{K}$ and univariate polynomial arithmetic).

Setting 2 is the vector rational reconstruction problem. We call PM-Basis on $[\mathbf{{F}}^{\mathsf{T}}\;\;-\mathbf{{I}}_{m}]^{\mathsf{T}}$ at order $(m+1)d$ with shift $(0,\ldots,0)$, where $\mathbf{{F}}$ is a $1\times m$ vector of degree $(m+1)d-1$, and we solve a block system with $1\times m$ Hankel blocks of size $md\times d$. The cost bounds are $\mathchoice{\tilde{O}\left(m^{\omega+1}d\right)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega+1}d)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega+1}d)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega+1}d)}$ and $\mathchoice{\tilde{O}\left(m^{\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}d)}$, respectively. Approximants are faster up to dimension about 15, which is explained by the arguments in the previous paragraph. For larger dimensions, as predicted by the cost estimates, the block-Hankel solver is more efficient.

3.2.   Interpolant bases.

For matrices $\mathbf{{E}}=(\mathbf{{E}}_{1},\dots,\mathbf{{E}}_{d})$ in $\mathbb{K}^{m\times n}$ and pairwise distinct points $\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{d})$ in $\mathbb{K}$, consider

[TABLE]

An interpolant basis for $(\mathbf{{E}},\boldsymbol{\alpha})$ is a matrix whose rows form a basis of the $\mathbb{K}[x]$-module $\mathcal{I}_{\boldsymbol{\alpha}}(\mathbf{{E}})$. Note that $\mathcal{I}_{\boldsymbol{\alpha}}(\mathbf{{F}}(\alpha_{1}),\dots,\mathbf{{F}}({\alpha_{d}}))$ coincides with $\mathcal{A}_{M}(\mathbf{{F}})$, for $\mathbf{{F}}$ in $\mathbb{K}[x]^{m\times n}$ and $M=\Pi_{i=1}^{d}(x-\alpha_{i})$.

This definition is a specialization of those in (Beckermann and Labahn, 2000; Jeannerod et al., 2017), which consider $n$ sets of points, one for each of the $n$ columns of $\mathbf{{E}}_{1},\dots,\mathbf{{E}}_{d}$: here, these sets are all equal. This more restrictive problem allows us to give faster algorithms than those in these references, by direct adaptations of the approximant basis algorithms presented above. Besides, Sections 4.1 and 4.2 will show that interpolant bases can often play the same role as approximant bases in applications.

These adaptations are described in Algorithms 4 and 5, where $\boldsymbol{\alpha}_{i\ldots j}$ stands for the sublist $(\alpha_{i},\alpha_{i+1},\ldots,\alpha_{j})$. In the next proposition, we assume that $\mathsf{MM}(n,d)$ is in $\Omega(n^{2}\mathsf{M}(d))$ (instead, one may add an extra term $O(n^{2}\mathsf{M}(d)\log(d))$ in the cost).

Proposition 3.0.

Algorithm 5* is correct. For input evaluation points in geometric progression, it costs $O(\mathsf{MM}^{\prime}(m,d))$ if $n\leq m$ and $O(\mathsf{MM}^{\prime}(m,d)+m^{\omega-1}nd\log(d))$ otherwise. For general evaluation points, an extra cost $O(m^{2}\mathsf{M}(d)\log^{2}(d))$ is incurred.*

(Correctness follows from Items (i) and (iii) of (Jeannerod et al., 2018, Lem. 2.4); the cost analysis is standard for such divide and conquer algorithms.)

Our current code uses the threshold $T=32$ in the divide and conquer PM-Basis and PM-IntBasis: beyond this point, they are faster than the iterative M-Basis and M-IntBasis. Unlike in most other functions, where elements of $\mathbb{K}[x]^{m\times n}$ are represented as matrices of polynomials (Mat<Vec<zz_p>> in NTL), in M-Basis and M-IntBasis we see them as polynomials with matrix coefficients (Vec<Mat<zz_p>>). Indeed, since these algorithms involve only matrix arithmetic over $\mathbb{K}$ (recall that $\deg(\mathbf{{Q}})\leq 1$), this turns out to be more cache-friendly and faster.

We implemented two variants for approximant bases: either the residual $\mathbf{{R}}$ is computed from $\mathbf{{P}}$ and $\mathbf{{F}}$ at each iteration, or we initialize a list of residuals with a copy of $\mathbf{{F}}$ and we update the whole list at each iteration using $\mathbf{{Q}}$. The second variant improves over the first when $n>m/2$, with significant savings when $n$ is close to $m$. For interpolant bases, this does not lead to any gain.

Timings are showed in the next table, for Algorithms M-Basis (M), M-IntBasis (M-I), PM-Basis (PM), PM-IntBasis for general points (PM-I) and for geometric points (PM-Ig). For approximants, we take a random input in $\mathbb{K}[x]^{m\times n}$ of degree $d-1$; for interpolants, we take $d$ random matrices in $\mathbb{K}^{m\times n}$. We focus on the common case $m\simeq 2n$, which arises for example in kernel algorithms (Section 4.2) and in fraction reconstruction, itself used in basis reduction (Section 4.5) and in the resultant algorithm of (Villard, 2018) (Section 5).

[TABLE]

Concerning iterative algorithms, we observe that interpolants are slightly faster than approximants, which is explained by the cost of computing the residual $\mathbf{{R}}$: it uses one Horner evaluation of $\mathbf{{P}}$ and one matrix product for interpolants, whereas for approximants it uses about $\min(k,\deg(\mathbf{{P}}))$ matrix products at iteration $k$.

As for the divide and conquer algorithms, interpolant bases with general points are slower, in some cases significantly, than the other two algorithms: although the complexity analysis predicted a disadvantage, we believe that our implementation of multipoint evaluation at general points could be improved to reduce this gap. For the other two algorithms, the comparison is less clear. There could be many factors at play here, but the main differences lie in the base case (Step 1) which calls the iterative algorithm, and in the computation of residuals (Step 3) which uses either middle products or geometric evaluation. It seems that FFT-based polynomial multiplication performs slightly better than geometric evaluation for small matrices and slightly worse for large matrices.

4. Higher-level algorithms

In this section we consider kernel bases, system solving, determinants, and basis reduction; we discuss algorithms which rely on multiplication, through approximant/interpolant bases and lifting techniques. For many of these algorithms, we are not aware of previous implementations or experimental comparisons.

4.1.   A note on matrix fraction reconstruction.

Given $\mathbf{{H}}$ in $\mathbb{K}(x)^{n\times n}$, a left fraction description of $\mathbf{{H}}$ is a pair of polynomial matrices $(\mathbf{{Q}},\mathbf{{R}})$ in $\mathbb{K}[x]^{n\times n}$ such that $\mathbf{{H}}=\mathbf{{Q}}^{-1}\mathbf{{R}}$. It is minimal if $\mathbf{{Q}}$ and $\mathbf{{R}}$ have unimodular left matrix GCD and $\mathbf{{Q}}$ is in reduced form (right fraction descriptions are defined similarly). Besides, $\mathbf{{H}}$ is said to be strictly proper if the numerator of each of its entries has degree less than the corresponding denominator.

Such a description of $\mathbf{{H}}$ is often computed from the power series expansion of $\mathbf{{H}}$ at sufficient precision, using an approximant basis. Yet, for resultant computations in Section 5.2, we would like to use an interpolant basis to obtain this description from sufficiently many values of $\mathbf{{H}}$. We now state the validity of this approach; this is a matrix version of rational function reconstruction (Gathen and Gerhard, 2013, Chap. 5.7).

Proposition 4.0.

*Let $\mathbf{{H}}$ be in $\mathbb{K}(x)^{n\times n}$ be strictly proper and suppose $\mathbf{{H}}$ admits left and right fraction descriptions of degrees at most $D$, for some $D\in\mathbb{Z}_{>0}$. For $M$ in $\mathbb{K}[x]$ of degree at least $2D$ and such that all denominators in $\mathbf{{H}}$ are invertible modulo $M$, define the matrix $\mathbf{{F}}=[\mathbf{{H}}\bmod M\;\;-\mathbf{{I}}_{n}]^{\mathsf{T}}\in\mathbb{K}[x]^{2n\times n}.$ Then, if $\mathbf{{P}}\in\mathbb{K}[x]^{2n\times 2n}$ is a $\boldsymbol{0}$-ordered weak Popov basis of $\mathcal{A}_{M}(\mathbf{{F}})$, the first $n$ rows of $\mathbf{{P}}$ form a matrix $[\mathbf{{Q}}\;\;\mathbf{{R}}]$ such that $(\mathbf{{Q}},\mathbf{{R}})$ is a minimal left fraction description of $\mathbf{{H}}$, with $\mathbf{{Q}}$ in $\boldsymbol{0}$-ordered weak Popov form. *

The proof given in (Giorgi et al., 2003, Lem 3.7) for the specific $M=x^{2D+1}$ extends to any modulus $M$; using an ordered weak Popov form (rather than a reduced form) allows us both to know a priori that the first $n$ rows are those of degree at most $D$, and to use degree $2D$ instead of $2D+1$ (since $\deg(\mathbf{{R}})<\deg(\mathbf{{Q}})$ is ensured by this form).

In particular, if $M=\prod_{i=1}^{2D}(x-\alpha_{i})$ for pairwise distinct points $(\alpha_{1},\ldots,\alpha_{2D})$, the interpolant basis algorithms in Section 3.2 give a minimal left fraction description of $\mathbf{{H}}$ from $\mathbf{{H}}(\alpha_{1}),\dots,\mathbf{{H}}(\alpha_{2D})$.

4.2.   Kernel basis.

We implemented two kernel basis algorithms: the first one, based on Lemma 4.0, finds the kernel basis from a single approximant basis at sufficiently large order; the second one, from (Zhou et al., 2012), uses several approximant bases at small order and combines recursively obtained kernel bases via multiplication. With a minor modification and no performance impact, the latter returns an $\boldsymbol{s}$-ordered weak Popov basis. In both cases, we designed variants which rely on interpolant bases instead of approximant bases.

Lemma 4.0.

Let $\mathbf{{F}}$ be in $\mathbb{K}[x]^{m\times n}$ of degree $d\geq 0$, let $\boldsymbol{s}$ be in $\mathbb{N}^{m}$, and let $\delta$ in $\mathbb{Z}_{>0}$ be an upper bound on the $\boldsymbol{s}$-degree of any $\boldsymbol{s}$-reduced left kernel basis of $\mathbf{{F}}$; for example, $\delta=nd+\max(\boldsymbol{s})-\min(\boldsymbol{s})+1$. Let $M$ be in $\mathbb{K}[x]$ of degree at least $\delta+d$, and $\mathbf{{P}}$ in $\mathbb{K}[x]^{m\times m}$ be an $\boldsymbol{s}$-reduced basis of $\mathcal{A}_{M}(\mathbf{{F}})$. Then, the submatrix of $\mathbf{{P}}$ formed by its rows of $\boldsymbol{s}$-degree less than $\delta$ is an $\boldsymbol{s}$-reduced left kernel basis for $\mathbf{{F}}$.

For a proof, see https://hal.archives-ouvertes.fr/hal-01995873v1/document.

In particular, one may find $\mathbf{{P}}$ via PM-IntBasis at $d+\delta$ points or via PM-Basis at order $d+\delta$; for $n\leq m$, this costs $O(\mathsf{MM}^{\prime}(m,d+\delta))$. The approximant-based direct approach is folklore (Zhou et al., 2012, Sec. 2.3), yet explicit statements in the literature focus on shifts linked to the degrees in $\mathbf{{F}}$, with better bounds $\delta$ (see (Zhou et al., 2012, Lem. 3.3),(Neiger et al., 2018, Lem. 4.3)).

The algorithm of (Zhou et al., 2012) is more efficient, at least when the entries of $\boldsymbol{s}$ are close to the corresponding row degrees of $\mathbf{{F}}$; for a uniform shift, it costs $\mathchoice{\tilde{O}\left(m^{\omega}\lceil nd/m\rceil\right)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}\lceil nd/m\rceil)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}\lceil nd/m\rceil)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}\lceil nd/m\rceil)}$ operations. We obtained significant practical improvements over the plain implementation of (Zhou et al., 2012, Algo. 1) thanks to the following observation: if $n\leq m/2$, for a vast majority of input $\mathbf{{F}}$, the approximant basis at Step 2 of (Zhou et al., 2012, Algo. 1), computed at order more than $2s$, contains the sought kernel basis. Furthermore, this can be easily tested by checking well-chosen degrees, and then the algorithm can exit early, avoiding the further recursive calls. We took advantage of this via the following modifications: we use order $2s+1$ rather than $3s$ (see (Zhou et al., 2012, Rmk. 3.5) for a discussion on this point), and when $n>m/2$ we directly reduce the number of columns via the divide and conquer scheme in (Zhou et al., 2012, Thm. 3.15).

The use of approximants here follows the idea in Lemma 4.0: row vectors of small degree which are in $\mathcal{A}_{M}(\mathbf{{F}})$ for a large degree $M$ must be in the kernel of $\mathbf{{F}}$. Thus, one can directly replace approximant bases with interpolant bases in (Zhou et al., 2012, Algo. 1), up to modifying Step 8 accordingly (dividing by the appropriate polynomial $M$).

Timings for both approaches are showed in the next table, for a random $\mathbf{{F}}$ of degree $d$. Except for the last row, the shift is uniform and, as expected, (Zhou et al., 2012, Algo. 1) is faster than the direct approach; the differences between interpolant and approximant variants follow those observed in Section 3. The last row corresponds to a shift yielding the kernel basis in Hermite form and shows, as expected, that the direct approach is faster for shifts that are far from uniform.

We note that (Zhou et al., 2012, Algo. 1) may use partial linearization if it computes matrix products with unbalanced degrees or approximant bases with $n\ll m$. We have not yet implemented this part of the algorithm, which may lead to slowdowns for some rare inputs.

[TABLE]

4.3.   Linear system solving.

For systems $\mathbf{{A}}\boldsymbol{\upsilon}=\boldsymbol{b}$, with $\mathbf{{A}}$ in $\mathbb{K}[x]^{m\times n}$, $\boldsymbol{b}$ in $\mathbb{K}[x]^{m\times 1}$ and $\boldsymbol{\upsilon}$ in $\mathbb{K}(x)^{n\times 1}$, we implemented two families of algorithms. The first one uses lifting techniques, assuming $\mathbf{{A}}$ is square, nonsingular, with $\mathbf{{A}}(0)$ invertible; the algorithm returns a pair $(\boldsymbol{u},f)$ in $\mathbb{K}[x]^{n\times 1}\times\mathbb{K}[x]$ such that $\mathbf{{A}}\boldsymbol{u}=f\boldsymbol{b}$ and $f$ has minimal degree. The second one uses a kernel basis and works for any input $\mathbf{{A}}$; under the assumptions above, it has a similar output.

Lifting techniques.

Under the above assumptions, our lifting algorithm is standard: if $\mathbf{{A}}$ and $\boldsymbol{b}$ have degree at most $d$, we first compute the truncated inverse $\mathbf{{S}}=\mathbf{{A}}^{-1}\bmod x^{d+1}$ by matrix Newton iteration (Schulz, 1933). Then, we use Dixon’s algorithm (Dixon, 1982) to compute $\boldsymbol{\upsilon}\bmod x^{2nd}=\mathbf{{A}}^{-1}\boldsymbol{b}\bmod x^{2nd}$; it consists of roughly $2n$ steps, each involving a matrix-vector product using either $\mathbf{{A}}$ or $\mathbf{{S}}$. Then, vector rational reconstruction is applied to recover $(\boldsymbol{u},f)$ from $\boldsymbol{\upsilon}$. The cost of this algorithm is $O(\mathsf{MM}(n,d))$ for the truncated inverse of $\mathbf{{A}}$ and $O(n^{3}\mathsf{M}(d))$ for Dixon’s algorithm; overall this is in $\mathchoice{\tilde{O}\left(n^{3}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{3}d)}{O\tilde{\leavevmode\nobreak\ }(n^{3}d)}{O\tilde{\leavevmode\nobreak\ }(n^{3}d)}$.

To reduce the exponent in $n$, Storjohann introduced the high-order lifting algorithm (Storjohann, 2003). The core of this algorithm is the computation of $\Theta(\log(n))$ slices $\mathbf{{S}}_{0},\mathbf{{S}}_{1},\dots$ of the power series expansion of $\mathbf{{A}}^{-1}$, where the coefficients of $\mathbf{{S}}_{i}$ are the coefficients of degree $(2^{i}-1)d-2^{i}+1,\dots,(2^{i}+1)d-2^{i}-1$ in $\mathbf{{A}}^{-1}$. These matrices are computed recursively, each step involving 4 matrix products; the other steps of the algorithm, that use these $\mathbf{{S}}_{i}$ to compute $\boldsymbol{\upsilon}\bmod x^{2nd}$, are cheaper, so the runtime is $O(\mathsf{MM}(n,d)\log(n))\subset\mathchoice{\tilde{O}\left(n^{\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega}d)}$.

Using kernel bases.

For this second approach, let $\mathbf{{A}}$ be any matrix in $\mathbb{K}[x]^{m\times n}$ and $\boldsymbol{b}$ be in $\mathbb{K}[x]^{m\times 1}$. The algorithm simply computes $\mathbf{{K}}\in\mathbb{K}[x]^{(n+1)\times k}$, a right kernel basis of the augmented matrix $[\mathbf{{A}}\;|\;\boldsymbol{b}]\in\mathbb{K}[x]^{m\times(n+1)}$. The matrix $\mathbf{{K}}$ generates, via $\mathbb{K}(x)$-linear combinations of its columns, all solutions $\boldsymbol{\upsilon}\in\mathbb{K}(x)^{n\times 1}$ to $\mathbf{{A}}\boldsymbol{\upsilon}=\boldsymbol{b}$.

In particular, if $\mathbf{{K}}$ is empty (i.e. $k=0$, which requires $m\geq n$), or if the last row of $\mathbf{{K}}$ is zero, then the system has no solution. Furthermore, if $\mathbf{{A}}$ is square and nonsingular, $\mathbf{{K}}$ has a single column $[\boldsymbol{u}^{\mathsf{T}}\;|\;f]^{\mathsf{T}}$, where $\boldsymbol{u}\in\mathbb{K}[x]^{n\times 1}$ and $f\in\mathbb{K}[x]$, with $f$ of minimal degree (otherwise, $\mathbf{{K}}$ would not be a basis).

In this context, the fastest known kernel algorithm is (Zhou et al., 2012, Algo. 1). To exploit it best, we choose the input shift $\boldsymbol{s}=(\boldsymbol{d},d)$, where $d=\deg(\boldsymbol{b})$ and $\boldsymbol{d}\in\mathbb{N}^{n}$ is the tuple of column degrees of $\mathbf{{A}}$ (zero columns of $\mathbf{{A}}$ are discarded while computing $\boldsymbol{d}$).

Implementation.

We implemented the approaches described above: lifting with Dixon’s algorithm, high-order lifting, and via kernel. The table below shows timings for randomly chosen $m\times m$ matrix $\mathbf{{A}}$ and $m\times 1$ vector $\boldsymbol{b}$, both of degree $d$. In this case the lifting algorithms apply (with high probability). On such inputs, Dixon’s algorithm usually does best. High-order lifting, although theoretically faster, is outperformed, mainly because it performs $\Theta(\log(n))$ matrix products (we will however see that this algorithm still plays an important role for basis reduction). The kernel based approach is moderately slower than Dixon’s algorithm, but has the advantage of working without any assumption on $\mathbf{{A}}$.

[TABLE]

4.4.   Determinant.

We implemented four algorithms, taking as input a square $m\times m$ matrix $\mathbf{{A}}$.

The most basic one uses expansion by minors, which turns out to be the fastest option up to dimension about 6.

The second one assumes that we have an element $\alpha$ in $\mathbb{K}$ of order at least $2md+1$, and uses evaluation/interpolation at the geometric progression $1,\alpha^{2},\ldots,\alpha^{2md}$; this costs $O(m^{3}\mathsf{M}(d)+m^{\omega+1}d)$ operations in $\mathbb{K}$. For dimensions between 7 and about 20, it is often the fastest variant, sometimes competing with the third.

The third one consists in solving a linear system with random right-hand side of degree $d$; this yields a solution $(\boldsymbol{u},f)$ with $f$ the sought determinant up to a constant (Pan, 1988). For dimensions exceeding 20, this is the fastest method (assuming that $\mathbf{{A}}(0)$ is invertible).

The last one, from (Labahn et al., 2017), is based on triangularization and runs in $\mathchoice{\tilde{O}\left(m^{\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}d)}$. Our implementation currently only supports the generic case where the Hermite form of $\mathbf{{A}}$ has diagonal $(1,\ldots,1,\det(\mathbf{{A}}))$, or in other words, all so-called row bases computed in that algorithm are the identity. This allows us to circumvent the temporary lack of a row basis implementation, while still being able to observe meaningful timings. Indeed, we believe that timings for a complete implementation called on such “generic” input will be similar to the timings presented here in many cases of interest, where one can easily detect whether the algorithm has made a wrong prediction about the row basis being identity (for example if the input matrix is reduced, which is the case if it is the denominator of a minimal fraction description). This recursive determinant algorithm calls expansion by minors as a base case for small dimensions; for larger dimensions, it is generally slightly slower than the third method.

The timings below are for a random $m\times m$ matrix $\mathbf{{A}}$ of degree $d$.

[TABLE]

4.5.   Basis reduction.

Our implementation of the algorithm of (Giorgi et al., 2003) takes as input a nonsingular matrix $\mathbf{{A}}\in\mathbb{K}[x]^{m\times m}$ of degree $d$ such that $\mathbf{{A}}(0)$ is invertible, and returns a reduced form of $\mathbf{{A}}$. The algorithm first computes a slice $\mathbf{{S}}$ of $2d$ consecutive coefficients of degree about $md$ in the power series expansion of $\mathbf{{A}}^{-1}$, then uses PM-Basis to reconstruct a fraction description $\mathbf{{S}}^{-1}=\mathbf{{R}}^{-1}\mathbf{{Q}}$, and then returns $\mathbf{{R}}$. A Las Vegas randomized version is mentioned in (Giorgi et al., 2003), to remove the assumption on $\mathbf{{A}}(0)$: we will implement it for large enough $\mathbb{K}$, but for smaller fields this requires to work in an extension of $\mathbb{K}$, which is currently beyond the scope of our work.

In our experiments, to create the input $\mathbf{{A}}$, we started from a random $m\times m$ matrix of degree $d/3$ (which is reduced with high probability), and we left-multiplied it by a lower unit triangular matrix and then by an upper one, both chosen at random of degree $d/3$. The following table shows timings for both steps, with the first step either based on Newton iteration or on high-order lifting; the displayed total time is when using the faster of the two. We conclude that for reduction, as opposed to the above observations for system solving, it is crucial to rely on high-order lifting. Indeed, it improves over Newton iteration already for dimension 8, and the gap becomes quite significant when the dimension grows.

[TABLE]

5. Applications to bivariate resultants

We conclude this paper with algorithms originating from Villard’s recent breakthrough on computing the determinant of structured polynomial matrices (Villard, 2018). Fix a field $\mathbb{K}$ and consider the two following questions: computing the resultant of two polynomials $F,G$ in $\mathbb{K}[x,z]$ with respect to $z$, and computing the characteristic polynomial of an element $A$ in $\mathbb{K}[z]/(P)$, for some $P$ in $\mathbb{K}[z]$.

The second problem is a particular case of the former, since the characteristic polynomial of $A$ modulo $P$ is the resultant of $x-A(z)$ and $P(z)$ with respect to $z$, up to a nonzero constant. Let $n$ be an upper bound on the degree in $z$ of the polynomials we consider, and $d$ be a bound on their degree in $x$ (so in the second problem, $d=1$). Villard proved that for generic inputs, both problems can be solved in $\mathchoice{\tilde{O}\left(n^{2-1/\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{2-1/\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{2-1/\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{2-1/\omega}d)}\subset\mathchoice{\tilde{O}\left(n^{1.58}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{1.58}d)}{O\tilde{\leavevmode\nobreak\ }(n^{1.58}d)}{O\tilde{\leavevmode\nobreak\ }(n^{1.58}d)}$ operations in $\mathbb{K}$. For the first problem, the best previous bound is $\mathchoice{\tilde{O}\left(n^{2}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{2}d)}{O\tilde{\leavevmode\nobreak\ }(n^{2}d)}{O\tilde{\leavevmode\nobreak\ }(n^{2}d)}$, obtained either by evaluation/interpolation techniques or Reischert’s algorithm (Reischert, 1997). For the second problem, the previous record was $\mathchoice{\tilde{O}\left(n^{\omega_{2}/2}\right)}{O\tilde{\leavevmode\nobreak\ }(n^{\omega_{2}/2})}{O\tilde{\leavevmode\nobreak\ }(n^{\omega_{2}/2})}{O\tilde{\leavevmode\nobreak\ }(n^{\omega_{2}/2})}$, where $\omega_{2}$ is the exponent of matrix multiplication in size $(s,s)\times(s,s^{2})$, with $\omega_{2}/2\leq 1.63$ (Le Gall and Urrutia, 2018). Note that these bounds apply to all inputs.

We show how the work we presented above allows us to put Villard’s ideas to practice, and outperform the previous state of the art for large input sizes. This is however not straightforward: in both cases, this required modifications of Villard’s original designs (for the second case, using an algorithm from (Neiger et al., 2019)).

5.1.   Overview of the approach.

In (Villard, 2018), Villard designed the following algorithm to find the determinant of a matrix $\mathbf{{P}}$ over $\mathbb{K}[x]$.

The parameter $m$ is chosen so as to minimize the theoretical cost. The correctness of the algorithm follows from the next properties, which do not hold for an arbitrary nonsingular $\mathbf{{P}}$: the matrix $\mathbf{{H}}$ is strictly proper and admits a left fraction description $\mathbf{{H}}=\mathbf{{Q}}^{-1}\mathbf{{R}}$ such that $\det(\mathbf{{P}})=\det(\mathbf{{Q}})$, for $\mathbf{{Q}}$ and $\mathbf{{R}}$ in $\mathbb{K}[x]^{m\times m}$ of degree at most $\lceil\nu/m\rceil d$ (see Section 4.1 for definitions). In (Villard, 2018), they are proved to hold for generic instances of the problems discussed here.

Once sufficiently many terms of the expansion of $\mathbf{{H}}$ have been obtained in Step 1, the denominator $\mathbf{{Q}}$ is recovered by an approximant basis algorithm and its determinant is computed by a general algorithm in Steps 2 and 3, which cost $\mathchoice{\tilde{O}\left(m^{\omega}(\nu d/m)\right)}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}(\nu d/m))}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}(\nu d/m))}{O\tilde{\leavevmode\nobreak\ }(m^{\omega}(\nu d/m))}$.

While the algorithm applies to any nonsingular matrix $\mathbf{{P}}$ satisfying the properties above, in general it does not improve over previously known methods (see Section 4.4). Indeed, the fastest known algorithm for obtaining $\bar{\mathbf{{H}}}$ costs $\mathchoice{\tilde{O}\left(\nu^{\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(\nu^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(\nu^{\omega}d)}{O\tilde{\leavevmode\nobreak\ }(\nu^{\omega}d)}$ operations via high-order lifting (see for example (Gupta et al., 2012, Thm. 1)).

However, sometimes $\mathbf{{P}}$ has some structure which helps to speed up the first step. Villard pointed out that when $\mathbf{{P}}$ is the Sylvester matrix of two bivariate polynomials, then $\mathbf{{P}}^{-1}$ is a Toeplitz-like matrix which can be described succinctly as $\mathbf{{P}}^{-1}=\mathbf{{L}}_{1}\mathbf{{U}}_{1}+\mathbf{{L}}_{2}\mathbf{{U}}_{2}$; here, $\mathbf{{L}}_{1}$, $\mathbf{{L}}_{2}$ (resp. $\mathbf{{U}}_{1}$, $\mathbf{{U}}_{2}$) are lower (resp. upper) triangular Toeplitz matrices with entries in $\mathbb{K}(x)$. Hence, we start by computing the first columns $\boldsymbol{c}_{1}$ of $\mathbf{{L}}_{1}$ and $\boldsymbol{c}_{2}$ of $\mathbf{{L}}_{2}$ as well as the first rows $\boldsymbol{r}_{1}$ of $\mathbf{{U}}_{1}$ and $\boldsymbol{r}_{2}$ of $\mathbf{{U}}_{2}$, all of them modulo $x^{2\lceil\nu/m\rceil d+1}$; then, $\bar{\mathbf{{H}}}$ is directly obtained via the above formula for $\mathbf{{P}}^{-1}$, using $\mathchoice{\tilde{O}\left(m\nu d\right)}{O\tilde{\leavevmode\nobreak\ }(m\nu d)}{O\tilde{\leavevmode\nobreak\ }(m\nu d)}{O\tilde{\leavevmode\nobreak\ }(m\nu d)}$ operations. Computing these rows and columns is done by solving systems with matrices $\mathbf{{P}}$ and $\mathbf{{P}}^{\mathsf{T}}$ and very simple right-hand sides (Villard, 2018, Prop. 5.1), with power series coefficients, in time $\mathchoice{\tilde{O}\left(\nu^{2}d/m\right)}{O\tilde{\leavevmode\nobreak\ }(\nu^{2}d/m)}{O\tilde{\leavevmode\nobreak\ }(\nu^{2}d/m)}{O\tilde{\leavevmode\nobreak\ }(\nu^{2}d/m)}$.

Altogether, taking $m=\nu^{1/\omega}$ minimizes the cost, yielding the runtime $\mathchoice{\tilde{O}\left(\nu^{2-1/\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(\nu^{2-1/\omega}d)}{O\tilde{\leavevmode\nobreak\ }(\nu^{2-1/\omega}d)}{O\tilde{\leavevmode\nobreak\ }(\nu^{2-1/\omega}d)}$. In the case of bivariate resultants described above, the Sylvester matrix of $F$ and $G$ has size $\nu=2n$, hence the cost bound $\mathchoice{\tilde{O}\left(n^{2-1/\omega}d\right)}{O\tilde{\leavevmode\nobreak\ }(n^{2-1/\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{2-1/\omega}d)}{O\tilde{\leavevmode\nobreak\ }(n^{2-1/\omega}d)}$.

5.2.   Resultant of generic bivariate polynomials.

We implemented the algorithm described in the previous section to compute the resultant of generic $F,G$ in $\mathbb{K}[x,z]$; first experiments showed that obtaining $\boldsymbol{c}_{1}$, $\boldsymbol{c}_{2}$, $\boldsymbol{r}_{1}$, $\boldsymbol{r}_{2}$ was a bottleneck. These vectors have power series entries and are solutions of linear systems whose matrix is the Sylvester matrix of $F$ and $G$ or its transpose: they were obtained via Hensel lifting techniques, following (Gathen and Gerhard, 2013, Ch. 15.4).

To get better performance, we designed a minor variant of Villard’s algorithm: instead of computing the power series expansion of $\mathbf{{H}}$ modulo $x^{\delta}$, where $\delta=2\lceil\nu/m\rceil d+1$, we compute values of $\mathbf{{H}}$ at $\delta$ points. We choose these points in geometric progression and use the interpolant basis algorithm of Section 3.2 to recover $\mathbf{{Q}}$ and $\mathbf{{N}}$, as detailed in Section 4.1. The value of $\mathbf{{H}}$ at $x=\alpha$ is computed following the same approach as above, but over $\mathbb{K}$ instead of $\mathbb{K}[[x]]$. In particular, our implementation directly relies on NTL’s extended GCD algorithm over $\mathbb{K}=\mathbb{F}_{p}$ to compute the vectors $\boldsymbol{c}_{1}$, $\boldsymbol{c}_{2}$, $\boldsymbol{r}_{1}$, $\boldsymbol{r}_{2}$.

The next table compares our implementation to the direct approach via evaluation/interpolation; note that the latter approach, while straightforward conceptually, is the state of the art.

[TABLE]

[TABLE]

For these running times, input polynomials were chosen at random with partial degree $n$ both in $x$ and in $z$; such polynomials have total degree $2n$, and their resultant has degree $2n^{2}$. The largest examples have quite significant sizes, but such degrees are not unheard-of in applications, as for instance in the genus-2 point counting algorithms of (Gaudry and Harley, 2000; Gaudry and Schost, 2004, 2012; Abélard et al., 2018). Overall, with $n=d$, we observe a crossover point around $n=400$. Besides, in a close match with the analysis above, the parameter $m$ was set to $\lceil n^{0.4}\rceil$ since this gave us the best runtimes. As an example, for $d=300$, the cost of each individual steps were 65s for computing structured inversions, and 40s for obtaining $\mathbf{{Q}}$ and its determinant, which is a good balance.

5.3.   Characteristic polynomial.

We consider the computation of the characteristic polynomial of an element $A$ in $\mathbb{K}[z]/(P)$, for some monic $P$ in $\mathbb{K}[z]$ of degree $n$. The algorithm we implemented, and which we sketch below, is from (Neiger et al., 2019) and assumes that $A$ and $P$ are generic.

As explained previously, this problem is a particular case of a bivariate resultant, but we rely on another point of view that allows for a better asymptotic cost. Indeed, the characteristic polynomial of $A$ modulo $P$ is by definition the characteristic polynomial of the matrix $\mathbf{{M}}$ of multiplication by $A$ modulo $P$. In other words, it is the determinant of the degree-1 matrix $\mathbf{{P}}=x\mathbf{{I}}-\mathbf{{M}}\in\mathbb{K}[x]^{n\times n}$.

The genericity assumption ensures that $\mathbf{{M}}$ is invertible, hence the power series expansion of $\mathbf{{P}}^{-1}$ is $\sum_{k\geq 0}-\mathbf{{M}}^{-k-1}x^{k}$. Here, we use the top-left $m\times m$ quadrant $\mathbf{{H}}$ of $\mathbf{{P}}^{-1}$; it has entries $h_{i,j}\in\mathbb{K}[[x]]$, where

[TABLE]

for $0\leq i,j<m$ and for all $k\geq 0$.

A direct implementation of this idea does not improve on the runtime given in Section 5.1, since it computes $A^{-k-1}\bmod P$ for all $0\leq k<\delta=2\lceil n/m\rceil$ and therefore costs $\Omega(n^{2}/m)$. It turns out that baby-steps giant-steps techniques allow one to compute $h_{i,j,k}$ for $0\leq i,j<m$ and $0\leq k<\delta$ in $\mathchoice{\tilde{O}\left(\delta^{(\omega-1)/2}n+mn\right)}{O\tilde{\leavevmode\nobreak\ }(\delta^{(\omega-1)/2}n+mn)}{O\tilde{\leavevmode\nobreak\ }(\delta^{(\omega-1)/2}n+mn)}{O\tilde{\leavevmode\nobreak\ }(\delta^{(\omega-1)/2}n+mn)}$ operations in $\mathbb{K}$. Taking $m=\lceil n^{1/3}\rceil$ minimizes the overall cost, resulting in the runtime $\mathchoice{\tilde{O}\left(n^{(\omega+2)/3}\right)}{O\tilde{\leavevmode\nobreak\ }(n^{(\omega+2)/3})}{O\tilde{\leavevmode\nobreak\ }(n^{(\omega+2)/3})}{O\tilde{\leavevmode\nobreak\ }(n^{(\omega+2)/3})}\subset O(n^{1.46})$.

The following table compares our implementation to NTL’s built-in characteristic polynomial algorithm, with random inputs $A$ and $P$. For such inputs, NTL uses Shoup’s algorithm for power projection (Shoup, 1994), which runs in time $\mathchoice{\tilde{O}\left(n^{(\omega+1)/2}\right)}{O\tilde{\leavevmode\nobreak\ }(n^{(\omega+1)/2})}{O\tilde{\leavevmode\nobreak\ }(n^{(\omega+1)/2})}{O\tilde{\leavevmode\nobreak\ }(n^{(\omega+1)/2})}$.

[TABLE]

Acknowledgements.

Hyun was supported by a Mitacs Globalink award; Neiger was supported by CNRS/INS2I’s program for young researchers; Schost was supported by an NSERC Discovery Grant.