Multiplying matrices over the integers

Multiplying an $m \times k$ by a $k \times n$ integer matrix with entries of $a$ and $b$ bits produces entries of up to $c = a + b + \log_2 k$ bits. Two parameters decide the best method: the dimension $d = \min(m, k, n)$ and the entry size. For small entries the cost is dominated by memory traffic and by the overhead of the fmpz representation; for large matrices of moderate entries, by the $mkn$ multiply-adds; for large entries, by the cost of each integer multiplication. fmpz_mat_mul implements the choice; the exact tests are in How the method is chosen.

The generic gr_mat_mul over the integers (gr_ctx_init_fmpz) calls fmpz_mat_mul. Conversely, several other modules reduce their matrix products to it: fmpq_mat_mul after clearing denominators of rows and columns, fmpz_mod_mat_mul before reducing, fmpz_poly_mat_mul_KS and nmod_poly_mat_mul_KS after Kronecker substitution, arb_mat_mul_block and nfloat's matrix multiplication on blocks of exponent-aligned mantissas, and padic_mat_mul.

The methods#

Machine words#

When all entries fit in a signed word minus two bits (small fmpz), the product can be computed without touching the multiprecision machinery: the kernels read the words, multiply them with one hardware multiplication, and accumulate in one, two or three words depending on the bound $c$. _fmpz_mat_mul_small_internal transposes $B$ into an array of words and processes blocks of 16 rows of $A$ against 128 entries of the inner dimension, which keeps the working set in cache, and splits the rows over threads. Entries of up to two words are handled likewise by _fmpz_mat_mul_double_word_internal, which multiplies two-word integers with four word products and accumulates in four or five words. For dimensions above about 1000 these lose to Strassen–Winograd (when $c$ fits in a word) or to multimodular multiplication.

Fewer multiplications#

For entries of thousands of bits, each multiplication of entries costs $M(b)$, much more than an addition, so algorithms that trade multiplications for additions win for small dimensions. Waksman's algorithm [Waksman1970] computes, for each pair of inner indices, the products $(a_{i,2l} + b_{2l+1,j})(a_{i,2l+1} + b_{2l,j})$, whose sum over $l$ is the dot product plus terms depending only on row $i$ or column $j$; it needs $mkn/2 + O(mk + kn)$ multiplications. Strassen–Winograd [Str1969; Win1971], in the operation order of Bodrato [Bodrato2010], needs 7 instead of 8 half-size products per level at the cost of 15 additions; FLINT applies one level and recurses through fmpz_mat_mul. See Strassen–Winograd multiplication.

Products modulo primes#

The entries of $C$ are less than $2^c$ in absolute value, so $C$ is determined by $C$ modulo moduli $p_1, \ldots, p_r$ with product at least $2^{c + 1}$ (or $2^c$ if all entries are nonnegative). The multimodular method reduces $A$ and $B$ modulo each $p_i$, multiplies with nmod_mat_mul, and reconstructs $C$ by the Chinese remainder theorem (Chinese remaindering). FLINT uses the modulus $2^{63}$, for which reduction is free, and primes above $2^{59}$, so $r \approx 1 + (c - 63)/59$; if $c \le 63$, the single modulus $2^c$.

Compared with the $d^3 M(b)$ of the classical algorithm, the matrix products cost $d^\omega$ word operations per prime, that is, about $d^\omega b/\beta$ in total, which for $b$ up to a few thousand bits is much less than $d^3 M(b)$, and the conversions are only quadratic in $d$. Since nmod_mat_mul runs at nearly the floating-point peak of the processor (Multiplying matrices over Z/nZ), this is the method of choice for large matrices.

Floating point#

fmpz_mat_mul_blas is a multimodular method with primes small enough for a double matrix product to be exact: residues in the symmetric range below $p/2$ with $k (p/2)^2 < 2^{53}$, so that every dot product is computed exactly (Exactness). That gives primes of about $26 - \frac12 \log_2 k$ bits (at most about 21 bits when $c > 200$, to keep the reconstruction cheap), and about $c/(26 - \frac12 \log_2 k)$ calls to flint_dgemm. If $a + b + \log_2 k \le 53$, a single dgemm on the entries suffices. With a tuned BLAS library or FLINT's own kernels (Floating-point matrix multiplication (gemm)) this is the fastest method for large matrices with entries up to a few hundred bits: fmpz_mat_mul uses it for $d > 50$ when $c \le 53$, and from about $d = 250$ to $1300$ depending on $c$.

Transformed entries#

For entries of $b \gg 1000$ bits, each product of entries costs a multiplication of large integers, which with an FFT means two forward transforms, a pointwise product and an inverse transform. In a matrix product each entry of $A$ takes part in $n$ products and each entry of $B$ in $m$, so transforming every entry once and adding the $k$ pointwise products of each entry of $C$ in the transformed domain saves almost all transforms.

fmpz_mat_mul_fft_small implements this with the transformed integer ring of gr, built on the small-prime FFT of fft_small.h (Sums of products and fixed operands). The transforms of the inputs must be stored, $(mk + kn) L$ numbers, which the function limits with flint_fft_small_max_transformed_ring_size, falling back to row and column blocks and, if necessary, to blocks of the inner dimension. fmpz_mat_mul uses it when the smaller entry size exceeds about $7000 + 0.15 d^2$ bits (10000 for $d \le 5$). The older fmpz_mat_mul_fft does the same with the Schönhage–Strassen transform of the fft.h module, for entries of at least 2000 bits.

Squaring and powers#

For $A^2$, the transformed method keeps one pool of transforms, and Bodrato's variant of Strassen–Winograd [Bodrato2010] uses squarings of half-size blocks where possible. fmpz_mat_sqr applies Bodrato's formulas for dimensions up to 3 (and 4 with large entries) and otherwise fmpz_mat_mul. fmpz_mat_pow uses binary powering.

Rational and polynomial matrices#

• fmpq_mat_mul multiplies $\operatorname{diag}(r)^{-1} A'$ and $B' \operatorname{diag}(s)^{-1}$, where $A'$ and $B'$ are $A$ and $B$ with rows, respectively columns, multiplied by their least common denominators $r_i$ and $s_j$, so that one integer product suffices.
• Polynomial matrices over $\mathbb{Z}$ and $\mathbb{Z}/n\mathbb{Z}$ are multiplied by Kronecker substitution into integer matrices, or by evaluation and interpolation; see Matrices of polynomials.

@misc{flint-doc-algorithms-linear-algebra-integer-matrix-multiplication,
  author       = {{The FLINT team}},
  title        = {{Multiplying matrices over the integers}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/linear-algebra/integer-matrix-multiplication.html},
  note         = {Generated 2026-11-02}
}