Linear algebra over generic rings and its specialisations

Dense linear algebra in FLINT comes in two layers, as polynomial arithmetic does (Polynomial arithmetic over generic rings and its specialisations). The generic module gr_mat implements matrix multiplication, Gaussian elimination and its fraction-free variant, triangular solving, determinants, ranks, echelon forms, nullspaces, inverses, characteristic and minimal polynomials and eigenvalues once, for entries in any ring given by a gr_ctx_t. The type-specific modules (nmod_mat, fmpz_mat, fmpq_mat, fmpz_mod_mat, fq_mat and its variants, arb_mat, acb_mat, ca_mat) provide kernels whose speed depends on the entry representation, above all matrix multiplication, and algorithms that exploit the structure of their ring, such as multimodular and p-adic methods over \mathbb{Z} and preconditioned elimination for balls.

Ring methods and the role of multiplication#

A few matrix operations are ring methods that a coefficient ring can override (Contexts, method tables and specialisation): multiplication (GR_METHOD_MAT_MUL), the LU decomposition, the determinant, the characteristic polynomial, pivot search, row reduction against an echelon basis (MAT_REDUCE_ROW, used by the minimal polynomial), triangular solving, the exponential and diagonalisation. The rest of gr_mat is written in terms of these and of ring arithmetic.

Matrix multiplication matters most. The recursive LU decomposition, triangular solving with many right-hand sides, inversion and the algorithms built on them do their work in matrix products, so a ring that overrides MAT_MUL with a fast algorithm makes them run in O(n^\omega) operations with the exponent of its multiplication. The default is the classical algorithm, with each entry one dot product (_gr_vec_dot); it thus benefits from fast dot products (delayed modular reduction, exact accumulation, one rounding per entry for floating-point and ball arithmetic). One-level Strassen-Winograd, Waksman and Rosowski multiplications are available for rings to install above a cutoff [Str1969; Win1971; Waksman1970; Rosowski2023]. The overrides by ring, from nmod_mat_mul over \mathbb{Z}/n\mathbb{Z} to the block algorithm of arb_mat_mul, are tabulated in Multiplying matrices over generic rings.

Gaussian elimination and LU decomposition#

Over a field, gr_mat_lu computes a decomposition PA = LU with a row permutation P, unit lower triangular L and upper triangular U (in reduced row echelon shape for singular or rectangular matrices). The default gr_mat_lu_generic uses classical elimination when the matrix has fewer than 5 rows or columns and otherwise the recursive algorithm, which splits the columns in two halves, decomposes the left half, updates the right half with one triangular solve and one matrix product, and recurses; with fast multiplication its cost is O(n^\omega) (Gaussian elimination and LU decomposition).

Pivots are chosen by the ring method MAT_FIND_NONZERO_PIVOT. The default takes the first entry that is provably nonzero, and if there is none returns GR_UNABLE when some entry cannot be decided and GR_DOMAIN otherwise; the floating-point rings (arf, acf, nfloat) choose the entry of largest magnitude for numerical stability, and ca works harder to decide zero tests. The rings fmpz_mod, mpn_mod, nfloat and nfloat_complex override the LU decomposition itself with their own cutoffs and kernels.

Over an integral domain that is not a field, the fraction-free LU decomposition (gr_mat_fflu) divides only exactly, using Sylvester's identity as in Bareiss' algorithm [Bar1968]: all intermediate entries are minors of the input, so their size is bounded and no fractions appear (Fraction-free elimination).

The generic functions built on these decompositions choose between them by the properties of the ring:

Determinants#

The determinant is a ring method. The generic choices are cofactor expansion up to 4 \times 4, and then:

• over any commutative ring, Berkowitz's division-free algorithm [Ber1984], with O(n^4) operations (gr_mat_det_generic, the default);
• over fields, the LU decomposition, with O(n^\omega) operations (gr_mat_det_generic_field, used by the floating-point rings);
• over integral domains, the fraction-free LU decomposition, with O(n^3) exact divisions (gr_mat_det_generic_integral_domain).

The integers and rationals use the multimodular algorithm of fmpz_mat_det (Multimodular and p-adic linear algebra over the integers), nmod calls nmod_mat_det, nmod8, nmod32, fmpz_mod and mpn_mod use LU with Berkowitz's algorithm as the fallback for composite moduli, the fq rings use LU, and balls use Gaussian elimination in ball arithmetic with approximate preconditioning for large matrices. Because the default is division-free, the determinant over a field without an override costs O(n^4) operations; Determinants and permanents gives the details. The permanent, by contrast, has no polynomial-time algorithm [Val1979]; gr_mat_permanent uses Ryser's and Glynn's formulas [Rys1963; Gly2010].

Characteristic and minimal polynomials#

The default characteristic polynomial is Berkowitz's algorithm, which is division-free and works over any commutative ring. gr_mat also implements Danilevsky's reduction to companion form [Dan1937], with O(n^3) operations over a field, the Faddeev-LeVerrier algorithm [FadSom1972] and its baby-step giant-step variant, which need division by 1, \ldots, n, and reduction to Hessenberg form by Gaussian or Householder transformations followed by the Hessenberg recurrence [Wil1965]. The rings override the method with their own choice: nmod uses Danilevsky's algorithm above dimension 8 and fmpz_mod from dimension 16, both with Berkowitz's algorithm as the fallback when a pivot is not invertible (composite moduli); the fq rings use Danilevsky's algorithm above dimension 12; and fmpz and fmpq use fmpz_mat_charpoly, which is Berkowitz's algorithm for small matrices and a multimodular algorithm otherwise.

The minimal polynomial over a field (gr_mat_minpoly_field) is computed from Krylov sequences v, Av, A^2v, \ldots for the unit vectors v, reducing each new vector against an echelon basis of the previous ones (the ring method MAT_REDUCE_ROW) to find the first linear relation, and taking the least common multiple of the relations.

Eigenvalues are the roots of the characteristic polynomial, computed with gr_poly_roots in the ring itself or in another ring (gr_mat_eigenvalues_other, for example the complex algebraic numbers for a rational matrix). For complex balls, acb overrides the diagonalisation method: an approximate eigendecomposition by the QR algorithm is certified with acb_mat_eig_simple, which succeeds when the eigenvalues are simple and well separated.

AlgorithmCharacteristic and minimal polynomialsBerkowitz's division-free algorithm and Danilevsky's reduction to companion form for the characteristic polynomial over \mathbb{Z}/n\mathbb{Z}, and the Krylov method with least common multiples that FLINT uses for the minimal polynomial.

Rings that are not fields, and inexact rings#

The generic functions never return a wrong answer when the ring lacks a property they need; they return a status instead. Over a ring with zero divisors such as \mathbb{Z}/6\mathbb{Z}, the elimination-based functions fail with GR_DOMAIN (or GR_UNABLE) when they meet a non-invertible pivot, while the division-free determinant, characteristic polynomial and adjugate still succeed. A ring \mathbb{Z}/n\mathbb{Z} that pretends to be a field runs the field algorithms and reports a factor of n if it finds one.

Over real and complex balls, a pivot whose ball contains zero cannot be used and cannot be excluded, so rank, solving and inversion return GR_UNABLE when the outcome depends on it, while an exactly zero pivot is recognised. Division-free algorithms and ball determinants always return an enclosure. Over floating-point rings (arf, acf, nfloat), values are treated as exact and pivoting by magnitude gives the usual numerical algorithms, without error bounds. The guide Linear algebra over general and inexact rings works through these cases with examples.

Specialised modules#

The specialised algorithms are described in Multimodular and p-adic linear algebra over the integers (integers and rationals), Multiplying matrices over the integers, Matrix multiplication modulo small integers, Lattice reduction (LLL), Hermite and Smith normal forms, Matrices of polynomials and Linear algebra over exact numbers (exact real and complex numbers). Over \mathbb{Z} and \mathbb{Q}, the ring methods make gr_mat_mul, gr_mat_det, gr_mat_charpoly and gr_mat_permanent call fmpz_mat_mul, fmpz_mat_det, _fmpz_mat_charpoly and fmpz_mat_permanent (or their fmpq_mat counterparts).

Some specialised algorithms are not reached through gr_mat: over arb and acb, gr_mat_nonsingular_solve and gr_mat_inv use the generic ball LU decomposition (as arb_mat_solve_lu does) rather than the preconditioned solver of arb_mat_solve, which gives much tighter enclosures for large well-conditioned systems; and over fmpz, gr_mat_nonsingular_solve_den uses fraction-free elimination rather than the p-adic solver of fmpz_mat_solve. For a fixed entry type and large inputs, the specialised module is then the better choice; the generic layer is needed for entry rings without a module of their own (polynomial rings, number fields, qqbar, ca), for code that works over several rings, and where the status codes are needed.

@misc{flint-doc-algorithms-linear-algebra-overview,
  author       = {{The FLINT team}},
  title        = {{Linear algebra over generic rings and its specialisations}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/linear-algebra/overview.html},
  note         = {Generated 2026-10-20}
}