alkahest-cas 3.12.0

//! [`MatGroup`] — a subgroup of `GL(d, q)` given by generators.

use std::collections::{HashMap, HashSet};
use std::sync::OnceLock;

use rug::ops::Pow;
use rug::Integer;

use super::chain::{build_chain, MatrixSiftResult, MatrixStabilizerChain, DEFAULT_BUDGET};
use super::element::{
    commutator, commutes, conjugate, entry_key, fe_sub, field_elements, field_order_u64,
    is_identity, projective_normalise, EntryKey,
};
use super::error::MatGroupError;
use super::{
    DEFAULT_MATGROUP_ELEMENT_CAP, MAX_MATGROUP_COMMUTANT_ELEMENTS, MAX_MATGROUP_DEGREE,
    MAX_MATGROUP_ORBIT,
};
use crate::ffield::{FieldElement, FiniteField, GfMatrix};
use crate::group::{Permutation, PermutationGroup};

/// A subgroup of `GL(d, q)`, given by a list of invertible `d × d` matrices.
///
/// The action is on **row** vectors, `v ↦ v·M`, so that a product of matrices
/// composes left to right and agrees with [`crate::group::Permutation::compose`]
/// — see the module docs.
///
/// The stabilizer chain is computed on first use and cached.
#[derive(Clone, Debug)]
pub struct MatGroup {
    field: FiniteField,
    degree: usize,
    generators: Vec<GfMatrix>,
    budget: u64,
    chain: OnceLock<Result<MatrixStabilizerChain, MatGroupError>>,
}

impl MatGroup {
    /// The subgroup of `GL(degree, q)` generated by `generators`.
    ///
    /// An empty generator list is the trivial group, which is a group — hence
    /// the explicit `degree`.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-005` for a degree above [`MAX_MATGROUP_DEGREE`];
    /// `E-MATGRP-001` / `-002` / `-003` for a generator that is not a `d × d`
    /// matrix over `field`; `E-MATGRP-004` for a singular generator.
    pub fn new(
        field: &FiniteField,
        degree: usize,
        generators: Vec<GfMatrix>,
    ) -> Result<Self, MatGroupError> {
        if degree == 0 || degree > MAX_MATGROUP_DEGREE {
            return Err(MatGroupError::DegreeTooLarge {
                degree,
                max: MAX_MATGROUP_DEGREE,
            });
        }
        for (index, g) in generators.iter().enumerate() {
            let (rows, cols) = g.shape();
            if rows != cols {
                return Err(MatGroupError::NotSquare { rows, cols });
            }
            if rows != degree {
                return Err(MatGroupError::DegreeMismatch {
                    expected: degree,
                    got: rows,
                });
            }
            if g.field() != field {
                return Err(MatGroupError::FieldMismatch {
                    expected: format!("{field}"),
                    got: format!("{}", g.field()),
                });
            }
            if g.determinant()?.is_zero() {
                return Err(MatGroupError::SingularGenerator { index });
            }
        }
        Ok(MatGroup {
            field: field.clone(),
            degree,
            generators,
            budget: DEFAULT_BUDGET,
            chain: OnceLock::new(),
        })
    }

    /// Like [`MatGroup::new`], taking the degree from the first generator.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-002` for an empty list: the degree of the trivial group is not
    /// recoverable from no generators, and guessing would silently build a
    /// group acting on the wrong space.
    pub fn from_generators(
        field: &FiniteField,
        generators: Vec<GfMatrix>,
    ) -> Result<Self, MatGroupError> {
        let degree =
            generators
                .first()
                .map(|g| g.nrows())
                .ok_or(MatGroupError::DegreeMismatch {
                    expected: 1,
                    got: 0,
                })?;
        Self::new(field, degree, generators)
    }

    /// The trivial subgroup of `GL(degree, q)`.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-005` for a degree outside `1..=`[`MAX_MATGROUP_DEGREE`].
    pub fn trivial(field: &FiniteField, degree: usize) -> Result<Self, MatGroupError> {
        Self::new(field, degree, Vec::new())
    }

    /// The same group with a different Schreier–Sims work budget.
    ///
    /// Clears the cached chain, so the next [`MatGroup::order`] rebuilds under
    /// the new budget.
    pub fn with_budget(&self, budget: u64) -> Self {
        MatGroup {
            field: self.field.clone(),
            degree: self.degree,
            generators: self.generators.clone(),
            budget: budget.max(1),
            chain: OnceLock::new(),
        }
    }

    /// The field.
    pub fn field(&self) -> &FiniteField {
        &self.field
    }

    /// The size of the matrices.
    pub fn degree(&self) -> usize {
        self.degree
    }

    /// The generators, as supplied.
    pub fn generators(&self) -> &[GfMatrix] {
        &self.generators
    }

    /// The Schreier–Sims work budget in force.
    pub fn budget(&self) -> u64 {
        self.budget
    }

    /// The identity matrix.
    ///
    /// # Errors
    ///
    /// `E-GFQ-012` for a shape FLINT cannot allocate.
    pub fn identity(&self) -> Result<GfMatrix, MatGroupError> {
        Ok(GfMatrix::identity(&self.field, self.degree)?)
    }

    /// Is every generator the identity?
    ///
    /// Read off the generators, with no stabilizer chain.
    pub fn is_trivial(&self) -> bool {
        self.generators.iter().all(|g| is_identity(&self.field, g))
    }

    /// The base and strong generating set, computed once and cached.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-006` when a basic orbit outgrows [`MAX_MATGROUP_ORBIT`];
    /// `E-MATGRP-007` when Schreier–Sims exhausts its budget.
    pub fn stabilizer_chain(&self) -> Result<&MatrixStabilizerChain, MatGroupError> {
        self.chain
            .get_or_init(|| build_chain(&self.field, self.degree, &self.generators, self.budget))
            .as_ref()
            .map_err(|e| e.clone())
    }

    /// `|G|`, exactly, from the stabilizer chain.
    ///
    /// This is computed from the **generators**, not from a closed-form
    /// formula: for the classical families the two are cross-checked against
    /// each other in this module's tests.
    ///
    /// ```
    /// use alkahest_cas::experimental::{FiniteField, MatGroup};
    ///
    /// let gf2 = FiniteField::prime(2).unwrap();
    /// // |GL(3,2)| = (8-1)(8-2)(8-4) = 168, found by Schreier-Sims.
    /// assert_eq!(MatGroup::general_linear(&gf2, 3).unwrap().order().unwrap(), 168);
    /// ```
    ///
    /// # Errors
    ///
    /// As [`MatGroup::stabilizer_chain`].
    pub fn order(&self) -> Result<Integer, MatGroupError> {
        Ok(self.stabilizer_chain()?.order())
    }

    /// Sift `m` through the chain.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::stabilizer_chain`], plus `E-MATGRP-001` / `-002` / `-003`
    /// for a matrix of the wrong shape or field.
    pub fn sift(&self, m: &GfMatrix) -> Result<MatrixSiftResult, MatGroupError> {
        self.stabilizer_chain()?.sift(m)
    }

    /// Is `m` an element of this group?
    ///
    /// A singular matrix is not in any subgroup of `GL(d, q)`, so it is
    /// `Ok(false)` rather than an error — but it must still be a `d × d` matrix
    /// over the right field, because a matrix of the wrong shape is a mistake
    /// rather than a non-member.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::sift`].
    pub fn contains(&self, m: &GfMatrix) -> Result<bool, MatGroupError> {
        let (rows, cols) = m.shape();
        if rows != cols {
            return Err(MatGroupError::NotSquare { rows, cols });
        }
        if rows != self.degree {
            return Err(MatGroupError::DegreeMismatch {
                expected: self.degree,
                got: rows,
            });
        }
        if m.field() != &self.field {
            return Err(MatGroupError::FieldMismatch {
                expected: format!("{}", self.field),
                got: format!("{}", m.field()),
            });
        }
        if m.determinant()?.is_zero() {
            return Ok(false);
        }
        Ok(self.sift(m)?.is_member())
    }

    /// The base, as standard-basis indices.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::stabilizer_chain`].
    pub fn base(&self) -> Result<Vec<usize>, MatGroupError> {
        Ok(self.stabilizer_chain()?.base())
    }

    /// Every strong generator, deduplicated.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::stabilizer_chain`].
    pub fn strong_generators(&self) -> Result<Vec<GfMatrix>, MatGroupError> {
        Ok(self.stabilizer_chain()?.strong_generators())
    }

    /// Every element of the group, refusing above
    /// [`DEFAULT_MATGROUP_ELEMENT_CAP`].
    ///
    /// # Errors
    ///
    /// `E-MATGRP-008` above the cap. The order is still exact in that case; it
    /// is the *list* that is refused.
    pub fn elements(&self) -> Result<Vec<GfMatrix>, MatGroupError> {
        self.elements_with_cap(DEFAULT_MATGROUP_ELEMENT_CAP)
    }

    /// [`MatGroup::elements`] with an explicit cap.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-008` above `cap`.
    pub fn elements_with_cap(&self, cap: u64) -> Result<Vec<GfMatrix>, MatGroupError> {
        self.stabilizer_chain()?.elements_with_cap(cap)
    }

    // -----------------------------------------------------------------------
    // Orbits
    // -----------------------------------------------------------------------

    fn check_vector(&self, v: &GfMatrix) -> Result<(), MatGroupError> {
        if v.nrows() != 1 {
            // Vectors here are 1 x d *rows*, because the action is v -> v*M.
            return Err(MatGroupError::DegreeMismatch {
                expected: 1,
                got: v.nrows(),
            });
        }
        if v.ncols() != self.degree {
            return Err(MatGroupError::DegreeMismatch {
                expected: self.degree,
                got: v.ncols(),
            });
        }
        if v.field() != &self.field {
            return Err(MatGroupError::FieldMismatch {
                expected: format!("{}", self.field),
                got: format!("{}", v.field()),
            });
        }
        Ok(())
    }

    /// The orbit of the row vector `v` under `v ↦ v·M`.
    ///
    /// Computed by breadth-first search over the generators, with no stabilizer
    /// chain, so it works whenever the orbit itself fits. The zero vector is
    /// fixed by everything and gives a one-point orbit.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-002` / `-003` for a vector of the wrong length or field;
    /// `E-MATGRP-006` when the orbit outgrows [`MAX_MATGROUP_ORBIT`].
    pub fn vector_orbit(&self, v: &GfMatrix) -> Result<Vec<GfMatrix>, MatGroupError> {
        self.check_vector(v)?;
        self.orbit_bfs(v.clone(), false)
    }

    /// The orbit of the projective point `[v]` — the 1-dimensional subspace
    /// spanned by `v` — under `[v] ↦ [v·M]`.
    ///
    /// Returned as canonical representatives: each vector scaled so its first
    /// non-zero coordinate is `1`. This is the action that keeps degrees
    /// manageable: `PG(d−1, q)` has `(q^d − 1)/(q − 1)` points rather than
    /// `q^d − 1`.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-012` for the zero vector; otherwise as
    /// [`MatGroup::vector_orbit`].
    pub fn projective_orbit(&self, v: &GfMatrix) -> Result<Vec<GfMatrix>, MatGroupError> {
        self.check_vector(v)?;
        let start = projective_normalise(&self.field, v)?;
        self.orbit_bfs(start, true)
    }

    fn orbit_bfs(&self, start: GfMatrix, projective: bool) -> Result<Vec<GfMatrix>, MatGroupError> {
        let mut points = vec![start.clone()];
        let mut seen: HashSet<EntryKey> = HashSet::new();
        seen.insert(entry_key(&self.field, &start));
        let mut i = 0;
        while i < points.len() {
            for g in &self.generators {
                let mut image = points[i].mul(g)?;
                if projective {
                    image = projective_normalise(&self.field, &image)?;
                }
                if seen.insert(entry_key(&self.field, &image)) {
                    if points.len() >= MAX_MATGROUP_ORBIT {
                        return Err(MatGroupError::OrbitTooLarge {
                            points: points.len() + 1,
                            max: MAX_MATGROUP_ORBIT,
                        });
                    }
                    points.push(image);
                }
            }
            i += 1;
        }
        Ok(points)
    }

    /// Every non-zero row vector of `GF(q)^d`, in the enumeration order this
    /// module uses for permutation actions.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-009` when `q` is above [`MAX_MATGROUP_FIELD_ORDER`](super::MAX_MATGROUP_FIELD_ORDER);
    /// `E-MATGRP-006` when `q^d − 1` is above [`MAX_MATGROUP_ORBIT`].
    pub fn nonzero_vectors(&self) -> Result<Vec<GfMatrix>, MatGroupError> {
        let alphabet = field_elements(&self.field)?;
        let q = alphabet.len();
        let total = q
            .checked_pow(u32::try_from(self.degree).unwrap_or(u32::MAX))
            .ok_or(MatGroupError::OrbitTooLarge {
                points: usize::MAX,
                max: MAX_MATGROUP_ORBIT,
            })?;
        if total - 1 > MAX_MATGROUP_ORBIT {
            return Err(MatGroupError::OrbitTooLarge {
                points: total - 1,
                max: MAX_MATGROUP_ORBIT,
            });
        }
        let mut out = Vec::with_capacity(total - 1);
        for code in 1..total {
            let mut rest = code;
            let mut entries = Vec::with_capacity(self.degree);
            for _ in 0..self.degree {
                entries.push(alphabet[rest % q].clone());
                rest /= q;
            }
            out.push(GfMatrix::from_elements(
                &self.field,
                1,
                self.degree,
                &entries,
            )?);
        }
        Ok(out)
    }

    /// The points of `PG(d−1, q)`: canonical representatives of the
    /// 1-dimensional subspaces, `(q^d − 1)/(q − 1)` of them.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::nonzero_vectors`].
    pub fn projective_points(&self) -> Result<Vec<GfMatrix>, MatGroupError> {
        let mut seen = HashSet::new();
        let mut out = Vec::new();
        for v in self.nonzero_vectors()? {
            let rep = projective_normalise(&self.field, &v)?;
            if seen.insert(entry_key(&self.field, &rep)) {
                out.push(rep);
            }
        }
        Ok(out)
    }

    /// The orbits of the group on the non-zero vectors of `GF(q)^d`.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::nonzero_vectors`].
    pub fn vector_orbits(&self) -> Result<Vec<Vec<GfMatrix>>, MatGroupError> {
        self.partition_into_orbits(self.nonzero_vectors()?, false)
    }

    /// The orbits of the group on the projective points.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::projective_points`].
    pub fn projective_orbits(&self) -> Result<Vec<Vec<GfMatrix>>, MatGroupError> {
        self.partition_into_orbits(self.projective_points()?, true)
    }

    fn partition_into_orbits(
        &self,
        points: Vec<GfMatrix>,
        projective: bool,
    ) -> Result<Vec<Vec<GfMatrix>>, MatGroupError> {
        let mut assigned: HashSet<EntryKey> = HashSet::new();
        let mut out = Vec::new();
        for p in points {
            let key = entry_key(&self.field, &p);
            if assigned.contains(&key) {
                continue;
            }
            let orbit = self.orbit_bfs(p, projective)?;
            for point in &orbit {
                assigned.insert(entry_key(&self.field, point));
            }
            out.push(orbit);
        }
        Ok(out)
    }

    // -----------------------------------------------------------------------
    // Induced permutation actions — the reuse of `crate::group`
    // -----------------------------------------------------------------------

    /// The induced action on the `q^d − 1` non-zero vectors, as a
    /// [`PermutationGroup`].
    ///
    /// This action is **faithful**: if `v·M = v` for every `v` then `M = I`. So
    /// the returned group has exactly [`MatGroup::order`] elements, and
    /// everything `crate::group` offers — orbits, its own Schreier–Sims,
    /// uniform random elements, element enumeration — applies unchanged.
    ///
    /// The point set is [`MatGroup::nonzero_vectors`], numbered from `0`.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::nonzero_vectors`]. Note that `q^d − 1` may be above
    /// [`crate::group::MAX_BSGS_DEGREE`], in which case the group is built but
    /// its `order()` refuses with `E-GRP-005`; use [`MatGroup::order`], which
    /// has no such limit.
    pub fn permutation_action_on_vectors(&self) -> Result<PermutationGroup, MatGroupError> {
        self.induced_action(self.nonzero_vectors()?, false)
    }

    /// The induced action on the `(q^d − 1)/(q − 1)` projective points, as a
    /// [`PermutationGroup`].
    ///
    /// This is `G` modulo the scalars it contains, so the permutation group's
    /// order is `|G| / |G ∩ Z(GL(d, q))|` and **not** `|G|` — that quotient is
    /// the point: `PSL(2, 7)` is the projective action of `SL(2, 7)` on the 8
    /// points of `PG(1, 7)`, of order `168`, while `|SL(2, 7)| = 336`.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::projective_points`].
    pub fn permutation_action_on_projective_points(
        &self,
    ) -> Result<PermutationGroup, MatGroupError> {
        self.induced_action(self.projective_points()?, true)
    }

    fn induced_action(
        &self,
        points: Vec<GfMatrix>,
        projective: bool,
    ) -> Result<PermutationGroup, MatGroupError> {
        let mut index: HashMap<EntryKey, usize> = HashMap::with_capacity(points.len());
        for (i, p) in points.iter().enumerate() {
            index.insert(entry_key(&self.field, p), i);
        }
        let mut permutations = Vec::with_capacity(self.generators.len());
        for g in &self.generators {
            let mut images = Vec::with_capacity(points.len());
            for p in &points {
                let mut image = p.mul(g)?;
                if projective {
                    image = projective_normalise(&self.field, &image)?;
                }
                let target =
                    *index
                        .get(&entry_key(&self.field, &image))
                        .ok_or(MatGroupError::Internal {
                            detail: "a generator moved a point outside the point set",
                        })?;
                images.push(target);
            }
            permutations.push(Permutation::from_images(images)?);
        }
        Ok(PermutationGroup::new(points.len(), permutations)?)
    }

    // -----------------------------------------------------------------------
    // Random elements — product replacement
    // -----------------------------------------------------------------------

    /// A pseudo-random element, by **product replacement** (the "rattle"
    /// variant of Celler–Leedham-Green–Murray–Niemeyer–O'Brien).
    ///
    /// Every value returned is a product of generators and their inverses, so it
    /// is an element of the group **by construction** — the randomness affects
    /// only the distribution, never membership. The distribution is not proved
    /// uniform; the algorithm is the standard heuristic and is used here for
    /// the same purpose it is used elsewhere, namely to get test inputs and
    /// candidate elements cheaply.
    ///
    /// `seed` is fed to a SplitMix64 counter, so the result is **reproducible**:
    /// the same seed on the same generating set always gives the same matrix.
    ///
    /// # Errors
    ///
    /// `E-GFQ-*` from the matrix arithmetic.
    pub fn random_element(&self, seed: u64) -> Result<GfMatrix, MatGroupError> {
        Ok(self
            .random_elements(seed, 1)?
            .into_iter()
            .next()
            .unwrap_or(self.identity()?))
    }

    /// `count` pseudo-random elements from a single product-replacement run.
    ///
    /// Cheaper than `count` calls to [`MatGroup::random_element`] with
    /// different seeds, because the accumulator state is shared.
    ///
    /// # Errors
    ///
    /// `E-GFQ-*` from the matrix arithmetic.
    pub fn random_elements(&self, seed: u64, count: usize) -> Result<Vec<GfMatrix>, MatGroupError> {
        let identity = self.identity()?;
        if self.generators.is_empty() {
            return Ok(vec![identity; count]);
        }
        // The standard recommendation is at least 10 accumulators and at least
        // twice the number of generators, with a burn-in of 50 steps per slot
        // before the first value is taken.
        let slots = self.generators.len().saturating_mul(2).max(10);
        let mut state: Vec<GfMatrix> = (0..slots)
            .map(|i| self.generators[i % self.generators.len()].clone())
            .collect();
        let mut accumulator = identity;
        let mut rng = SplitMix64::new(seed);

        for _ in 0..(50 * slots) {
            product_replacement_step(&mut state, &mut accumulator, &mut rng)?;
        }
        let mut out = Vec::with_capacity(count);
        for _ in 0..count {
            for _ in 0..slots {
                product_replacement_step(&mut state, &mut accumulator, &mut rng)?;
            }
            out.push(accumulator.clone());
        }
        Ok(out)
    }

    // -----------------------------------------------------------------------
    // Derived subgroup, normal closure, centre
    // -----------------------------------------------------------------------

    /// The normal closure of `subset` in this group: the smallest normal
    /// subgroup of `G` containing `⟨subset⟩`.
    ///
    /// Computed by closing `⟨subset⟩` under conjugation by the generators of
    /// `G` until nothing new appears. Conjugating by the generators alone is
    /// enough at the fixed point: if `xᵍ ∈ N` for every `x ∈ N` and every
    /// generator `g`, then conjugation by `g` is an injective endomorphism of
    /// the *finite* group `N`, hence onto, so `N` is also closed under
    /// conjugation by `g⁻¹`.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::stabilizer_chain`] — membership is tested with a chain,
    /// and it is rebuilt each time the closure grows.
    pub fn normal_closure(&self, subset: &[GfMatrix]) -> Result<MatGroup, MatGroupError> {
        let mut generators: Vec<GfMatrix> = Vec::new();
        let mut seen: HashSet<EntryKey> = HashSet::new();
        for m in subset {
            if is_identity(&self.field, m) {
                continue;
            }
            if seen.insert(entry_key(&self.field, m)) {
                generators.push(m.clone());
            }
        }
        let mut closure =
            MatGroup::new(&self.field, self.degree, generators.clone())?.with_budget(self.budget);
        loop {
            let mut grew = false;
            let current = closure.generators().to_vec();
            for g in &self.generators {
                for h in &current {
                    let c = conjugate(g, h)?;
                    if is_identity(&self.field, &c) {
                        continue;
                    }
                    if !closure.contains(&c)? {
                        generators.push(c);
                        closure = MatGroup::new(&self.field, self.degree, generators.clone())?
                            .with_budget(self.budget);
                        grew = true;
                    }
                }
            }
            if !grew {
                return Ok(closure);
            }
        }
    }

    /// The derived subgroup `G' = [G, G]`.
    ///
    /// The normal closure of the commutators of the generators, which is the
    /// standard identity `[G, G] = ⟨[gᵢ, gⱼ]⟩^G` — the commutators of a
    /// generating set generate `G'` only after normal closure, not on their own.
    ///
    /// For most `(d, q)`, `GL(d, q)' = SL(d, q)`; `GL(2, 2)` is the smallest
    /// exception, where the derived subgroup has order 3 and `SL(2, 2)` has
    /// order 6.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::normal_closure`].
    pub fn derived_subgroup(&self) -> Result<MatGroup, MatGroupError> {
        let mut commutators = Vec::new();
        for (i, g) in self.generators.iter().enumerate() {
            for h in self.generators.iter().skip(i + 1) {
                let c = commutator(g, h)?;
                if !is_identity(&self.field, &c) {
                    commutators.push(c);
                }
            }
        }
        self.normal_closure(&commutators)
    }

    /// Is this group perfect, i.e. `G' = G`?
    ///
    /// # Errors
    ///
    /// As [`MatGroup::derived_subgroup`].
    pub fn is_perfect(&self) -> Result<bool, MatGroupError> {
        Ok(self.derived_subgroup()?.order()? == self.order()?)
    }

    /// A basis for the algebra of matrices commuting with every generator —
    /// the endomorphism algebra of the natural module,
    /// `End_{GF(q)G}(GF(q)^d)`.
    ///
    /// This is a **linear** computation: `MX = XM` for each generator `M` is a
    /// homogeneous linear system in the `d²` entries of `X`, and the answer is
    /// its nullspace. It is what makes the centre computable without a
    /// backtrack search — `Z(G) = G ∩ End_{GF(q)G}(GF(q)^d)`, because commuting
    /// with the generators is the same as commuting with all of `G`.
    ///
    /// For an absolutely irreducible module the algebra is the scalars and this
    /// returns one matrix.
    ///
    /// # Errors
    ///
    /// `E-GFQ-*` from the nullspace computation.
    pub fn commutant_basis(&self) -> Result<Vec<GfMatrix>, MatGroupError> {
        let d = self.degree;
        let cells = d * d;
        if self.generators.is_empty() {
            // Everything commutes with nothing: the algebra is all of M_d(q),
            // spanned by the matrix units.
            let mut out = Vec::with_capacity(cells);
            for a in 0..d {
                for b in 0..d {
                    let mut entries = vec![self.field.zero(); cells];
                    entries[a * d + b] = self.field.one();
                    out.push(GfMatrix::from_elements(&self.field, d, d, &entries)?);
                }
            }
            return Ok(out);
        }

        // Row (g, i, j) says (GX - XG)[i][j] = 0; column (a, b) is the unknown
        // X[a][b]. The coefficient of X[a][b] is  delta(b,j)*G[i][a] - delta(a,i)*G[b][j].
        let mut rows: Vec<FieldElement> = Vec::with_capacity(self.generators.len() * cells * cells);
        for g in &self.generators {
            let ge = g.to_elements();
            for i in 0..d {
                for j in 0..d {
                    for a in 0..d {
                        for b in 0..d {
                            let left = if b == j {
                                ge[i * d + a].clone()
                            } else {
                                self.field.zero()
                            };
                            let right = if a == i {
                                ge[b * d + j].clone()
                            } else {
                                self.field.zero()
                            };
                            rows.push(fe_sub(&self.field, &left, &right)?);
                        }
                    }
                }
            }
        }
        let system =
            GfMatrix::from_elements(&self.field, self.generators.len() * cells, cells, &rows)?;
        let kernel = system.nullspace()?;
        let columns = kernel.ncols();
        let ke = kernel.to_elements();
        let mut out = Vec::with_capacity(columns);
        for c in 0..columns {
            let mut entries = Vec::with_capacity(cells);
            for r in 0..cells {
                entries.push(ke[r * columns + c].clone());
            }
            out.push(GfMatrix::from_elements(&self.field, d, d, &entries)?);
        }
        Ok(out)
    }

    /// Every element of the centre `Z(G)`.
    ///
    /// Two routes, whichever is affordable:
    ///
    /// * if `|G|` is at or below [`DEFAULT_MATGROUP_ELEMENT_CAP`], enumerate
    ///   `G` and keep the elements that commute with every generator;
    /// * otherwise intersect `G` with [`MatGroup::commutant_basis`]'s algebra,
    ///   enumerating its `q^dim` elements and keeping the invertible ones that
    ///   are in `G`.
    ///
    /// The second route is what makes `Z(GL(4, 3))` cheap: the algebra is
    /// 1-dimensional, so it is three candidates rather than `|GL(4, 3)|`.
    ///
    /// # Errors
    ///
    /// `E-MATGRP-011` when neither route fits — the group is too big to list
    /// and its centralizing algebra too big to enumerate.
    pub fn centre_elements(&self) -> Result<Vec<GfMatrix>, MatGroupError> {
        let order = self.order()?;
        if order <= DEFAULT_MATGROUP_ELEMENT_CAP {
            let mut out = Vec::new();
            for m in self.elements()? {
                let mut central = true;
                for g in &self.generators {
                    if !commutes(&m, g)? {
                        central = false;
                        break;
                    }
                }
                if central {
                    out.push(m);
                }
            }
            return Ok(out);
        }

        let basis = self.commutant_basis()?;
        let dimension = basis.len();
        let q = field_order_u64(&self.field)?;
        let total = Integer::from(q).pow(u32::try_from(dimension).unwrap_or(u32::MAX));
        let count = total.to_u64().unwrap_or(u64::MAX);
        if count > MAX_MATGROUP_COMMUTANT_ELEMENTS {
            return Err(MatGroupError::CommutantTooLarge {
                dimension,
                elements: total.to_string(),
                cap: MAX_MATGROUP_COMMUTANT_ELEMENTS,
            });
        }
        let alphabet = field_elements(&self.field)?;
        let mut out = Vec::new();
        for code in 0..count {
            let mut rest = code;
            let mut candidate = GfMatrix::zeros(&self.field, self.degree, self.degree)?;
            for b in &basis {
                let coefficient = &alphabet[(rest % q) as usize];
                rest /= q;
                if coefficient.is_zero() {
                    continue;
                }
                candidate = candidate.add(&b.scalar_mul(coefficient)?)?;
            }
            if candidate.determinant()?.is_zero() {
                continue;
            }
            if self.contains(&candidate)? {
                out.push(candidate);
            }
        }
        Ok(out)
    }

    /// The centre `Z(G)`, as a group generated by its own elements.
    ///
    /// # Errors
    ///
    /// As [`MatGroup::centre_elements`].
    pub fn centre(&self) -> Result<MatGroup, MatGroupError> {
        let elements = self.centre_elements()?;
        MatGroup::new(&self.field, self.degree, elements).map(|g| g.with_budget(self.budget))
    }
}

fn product_replacement_step(
    state: &mut [GfMatrix],
    accumulator: &mut GfMatrix,
    rng: &mut SplitMix64,
) -> Result<(), MatGroupError> {
    let slots = state.len();
    let i = (rng.next() % slots as u64) as usize;
    let mut j = (rng.next() % slots as u64) as usize;
    if j == i {
        j = (j + 1) % slots;
    }
    let invert = rng.next() & 1 == 1;
    let from_left = rng.next() & 1 == 1;
    let factor = if invert {
        state[j].inverse()?
    } else {
        state[j].clone()
    };
    let updated = if from_left {
        factor.mul(&state[i])?
    } else {
        state[i].mul(&factor)?
    };
    state[i] = updated;
    *accumulator = accumulator.mul(&state[i])?;
    Ok(())
}

/// SplitMix64, so that `random_element(seed)` is reproducible across platforms
/// without pulling in an RNG crate.
struct SplitMix64 {
    state: u64,
}

impl SplitMix64 {
    fn new(seed: u64) -> Self {
        SplitMix64 {
            state: seed.wrapping_add(0x9E37_79B9_7F4A_7C15),
        }
    }

    fn next(&mut self) -> u64 {
        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.state;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }
}

impl std::fmt::Display for MatGroup {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(
            f,
            "<matrix group of degree {} over {} with {} generator{}>",
            self.degree,
            self.field,
            self.generators.len(),
            if self.generators.len() == 1 { "" } else { "s" }
        )
    }
}