Computing Smith Forms Modulo $p^2$ of Sparse Matrices Faster Than Matrix Multiplication

Let $p$ be a prime and $R=\mathbb{Z}/p^2\mathbb{Z}$ the ring of integers modulo $p^2$. Any $A\in R^{n\times n}$ is unimodularly equivalent to its Smith form $$S=diag\bigl(\underbrace{1,\ldots,1}_{r_0}, \underbrace{p,\ldots,p}_{r_1}, \underbrace{0,\ldots,0}_{r_2}\bigr) \in R^{n\times n},$$ i.e., there exist $U,V\in R^{n\times n}$ such that $UAV=S$, with $\det U,\det V\in R^*$ (where $R^*$ is the set of units in $R$, elements not equivalent to $0\bmod p$). Our goal in this paper is to determine $r_0,r_1,r_2$ quickly when $A$ is sparse or structured. By ``sparse'' we mean $A$ is given by a black box such that for any $v\in R^{n\times 1}$ we can compute $v\mapsto Av$ with $\tilde{O}(n)$ operations in $R$, which captures having few nonzero elements or a multiplicative structure (e.g., Hankel or Toeplitz matrices). We present a randomized algorithm which requires an expected number of $$\tilde{O}\left(n^{3-1/(\omega-1)}\right)$$ operations in $R$ to compute the Smith form, where $\omega$ is the exponent of dense matrix multiplication. Using standard cubic matrix multiplication ($\omega=3$) our algorithm thus requires $\tilde{O}(n^{2.5})$ operations in $R$, while using the current asymptotically fastest matrix multiplication, with $\omega<2.371339$, our algorithm requires $\tilde{O}(n^{2.270786})$ operations in $R$. Our algorithm is probabilistic of the Monte Carlo type, meaning it fails on any invocation with controllably small probability. We employ iterative block-Wiedemann-style matrix techniques and structured preconditioners. To our knowledge, this is the first algorithm to compute the modular Smith Normal Form modulo $p^2$ requiring fewer than $\tilde{O}(n^\omega)$ operations in $R$, i.e., faster than any dense algorithm.