puremp 0.2.4

//! Finite (Galois) field extensions `GF(pᵏ)`.
//!
//! [`ModInt`](crate::mod_int::ModInt) already models the prime field `GF(p) =
//! ℤ/pℤ`. This module builds the *extension* fields `GF(pᵏ)` for `k ≥ 1`,
//! represented as `GF(p)[x] / (f)` for a monic irreducible modulus polynomial
//! `f` of degree `k` over `GF(p)`.
//!
//! A [`GaloisField`] is the shared field context: the prime `p`, the degree
//! `k`, and the modulus `f` (as coefficients low-to-high, each reduced into
//! `[0, p)`, leading coefficient `1`). The data is wrapped in an
//! [`Rc`] so [`GfElement`]s produced from one field share it
//! cheaply, exactly as [`ModInt`](crate::mod_int::ModInt) shares its modulus.
//!
//! A [`GfElement`] is the residue class of a polynomial, stored as its canonical
//! representative of degree `< k` (a length-`k` coefficient vector, each entry in
//! `[0, p)`). Elements of *different* fields must not be mixed: every binary
//! operation panics on a field mismatch, mirroring [`ModInt`](crate::mod_int::ModInt).
//!
//! Irreducibility is decided by **Rabin's test** (Lidl & Niederreiter, *Finite
//! Fields*, Thm. 3.x; Menezes et al., *Handbook of Applied Cryptography* §4.5.1):
//! a monic degree-`k` polynomial `f` is irreducible over `GF(p)` iff
//! `x^(pᵏ) ≡ x (mod f)` and `gcd(x^(p^{k/q}) − x, f) = 1` for every prime
//! divisor `q` of `k`. The powers `x^(pᵐ)` are built by iterating the Frobenius
//! (`p`-th power modulo `f`) `m` times.
//!
//! This is a clean-room implementation drawn from the open literature (Lidl &
//! Niederreiter; the HAC §2.6/§4.5; Cohen, *A Course in Computational Algebraic
//! Number Theory* §3.4–3.6; von zur Gathen & Gerhard, *Modern Computer
//! Algebra*); it copies no external source.

use alloc::rc::Rc;
use alloc::vec::Vec;
use core::fmt;

use crate::int::Int;

// ===========================================================================
// Internal `GF(p)[x]` arithmetic on `Vec<Int>` (coefficients low-to-high, each
// reduced into `[0, p)`, trailing zeros trimmed so the leading entry is nonzero;
// the empty vector is the zero polynomial).
// ===========================================================================

/// Highest index of a nonzero coefficient, or `None` for the zero polynomial.
fn deg_of(a: &[Int]) -> Option<usize> {
    a.iter().rposition(|c| !c.is_zero())
}

/// Drops trailing (high-degree) zero coefficients.
fn poly_trim(mut a: Vec<Int>) -> Vec<Int> {
    while a.last().is_some_and(Int::is_zero) {
        a.pop();
    }
    a
}

/// Pads (or truncates a trimmed poly) to exactly `len` coefficients with zeros.
fn pad_to(mut a: Vec<Int>, len: usize) -> Vec<Int> {
    a.truncate(len);
    while a.len() < len {
        a.push(Int::ZERO);
    }
    a
}

/// Coefficientwise `self − rhs` modulo `p`.
fn poly_sub(a: &[Int], b: &[Int], p: &Int) -> Vec<Int> {
    let n = a.len().max(b.len());
    let mut r = Vec::with_capacity(n);
    for i in 0..n {
        let x = a.get(i).cloned().unwrap_or(Int::ZERO);
        let y = b.get(i).cloned().unwrap_or(Int::ZERO);
        r.push(x.sub(&y).rem_euclid(p));
    }
    poly_trim(r)
}

/// Schoolbook `self · rhs`, reducing every coefficient modulo `p`.
fn poly_mul(a: &[Int], b: &[Int], p: &Int) -> Vec<Int> {
    if a.is_empty() || b.is_empty() {
        return Vec::new();
    }
    let mut r = alloc::vec![Int::ZERO; a.len() + b.len() - 1];
    for (i, ai) in a.iter().enumerate() {
        if ai.is_zero() {
            continue;
        }
        for (j, bj) in b.iter().enumerate() {
            r[i + j] = r[i + j].add(&ai.mul(bj));
        }
    }
    poly_trim(r.into_iter().map(|c| c.rem_euclid(p)).collect())
}

/// `self · scalar` modulo `p`.
fn poly_scale(a: &[Int], s: &Int, p: &Int) -> Vec<Int> {
    poly_trim(a.iter().map(|c| c.mul(s).rem_euclid(p)).collect())
}

/// Divides `a` by nonzero `b` over `GF(p)`, returning `(quotient, remainder)`
/// with `deg(remainder) < deg(b)`. The leading coefficient of `b` is inverted
/// modulo the prime `p`.
fn poly_divmod(a: &[Int], b: &[Int], p: &Int) -> (Vec<Int>, Vec<Int>) {
    let db = deg_of(b).expect("poly_divmod: division by zero polynomial");
    let inv_lc = b[db]
        .modinv(p)
        .expect("poly_divmod: leading coefficient must be invertible mod p");
    let mut r = poly_trim(a.to_vec());
    let rd = match deg_of(&r) {
        Some(d) if d >= db => d,
        _ => return (Vec::new(), r),
    };
    let mut q = alloc::vec![Int::ZERO; rd - db + 1];
    while let Some(rd) = deg_of(&r) {
        if rd < db {
            break;
        }
        let d = rd - db;
        let coef = r[rd].mul(&inv_lc).rem_euclid(p);
        for (j, bj) in b.iter().enumerate() {
            r[d + j] = r[d + j].sub(&coef.mul(bj)).rem_euclid(p);
        }
        q[d] = coef;
        r = poly_trim(r);
    }
    (poly_trim(q), r)
}

/// Remainder of `a` divided by nonzero `b` over `GF(p)`.
fn poly_rem(a: &[Int], b: &[Int], p: &Int) -> Vec<Int> {
    poly_divmod(a, b, p).1
}

/// Monic greatest common divisor of `a` and `b` over `GF(p)` (Euclid).
fn poly_gcd(a: &[Int], b: &[Int], p: &Int) -> Vec<Int> {
    let mut a = poly_trim(a.to_vec());
    let mut b = poly_trim(b.to_vec());
    while !b.is_empty() {
        let r = poly_rem(&a, &b, p);
        a = b;
        b = r;
    }
    // Normalize to monic.
    match deg_of(&a) {
        None => a,
        Some(d) => {
            if a[d].is_one() {
                a
            } else {
                let inv = a[d]
                    .modinv(p)
                    .expect("poly_gcd: leading coefficient invertible mod p");
                poly_scale(&a, &inv, p)
            }
        }
    }
}

/// `base^e mod modulus` in `GF(p)[x]` (square-and-multiply, exponent an [`Int`]).
fn poly_powmod(base: &[Int], e: &Int, modulus: &[Int], p: &Int) -> Vec<Int> {
    let mut result = poly_rem(&[Int::ONE], modulus, p);
    let mut b = poly_rem(base, modulus, p);
    for i in 0..e.bit_len() {
        if e.bit(i) {
            result = poly_rem(&poly_mul(&result, &b, p), modulus, p);
        }
        b = poly_rem(&poly_mul(&b, &b, p), modulus, p);
    }
    result
}

/// `x^(pᵐ) mod f`, built by iterating the Frobenius (`p`-th power mod `f`) `m`
/// times starting from `x`.
fn frobenius_pow(f: &[Int], m: usize, p: &Int) -> Vec<Int> {
    let x = alloc::vec![Int::ZERO, Int::ONE];
    let mut h = poly_rem(&x, f, p);
    for _ in 0..m {
        h = poly_powmod(&h, p, f, p);
    }
    h
}

/// Distinct prime divisors of `k` (ascending).
fn distinct_primes(k: usize) -> Vec<usize> {
    let mut out: Vec<usize> = Vec::new();
    for f in Int::from_u64(k as u64).factorize() {
        let v = f.to_u64().expect("small factor fits in u64") as usize;
        if out.last() != Some(&v) {
            out.push(v);
        }
    }
    out
}

/// Rabin's irreducibility test for a **monic** polynomial `f` of degree `≥ 1`
/// over `GF(p)`.
fn is_irreducible(f: &[Int], p: &Int) -> bool {
    let k = match deg_of(f) {
        Some(d) if d >= 1 => d,
        _ => return false,
    };
    if k == 1 {
        return true; // every monic linear polynomial is irreducible
    }
    let x = alloc::vec![Int::ZERO, Int::ONE];
    // (2) gcd(x^(p^{k/q}) − x, f) = 1 for every prime q | k.
    for q in distinct_primes(k) {
        let h = frobenius_pow(f, k / q, p);
        let diff = poly_sub(&h, &x, p);
        let g = poly_gcd(&diff, f, p);
        if deg_of(&g).is_some_and(|d| d >= 1) {
            return false; // a nontrivial common factor: reducible
        }
    }
    // (1) x^(pᵏ) ≡ x (mod f).
    let h = frobenius_pow(f, k, p);
    poly_sub(&h, &x, p).is_empty()
}

/// Extended Euclid in `GF(p)[x]`: the inverse `u` of `a` modulo `modulus`
/// (`u·a ≡ 1`), or `None` when `a ≡ 0` or `a` is not invertible.
fn poly_inv(a: &[Int], modulus: &[Int], p: &Int) -> Option<Vec<Int>> {
    let a = poly_trim(a.to_vec());
    if a.is_empty() {
        return None;
    }
    // Track only the coefficient of `a`: sᵢ·a ≡ rᵢ (mod modulus).
    let (mut r0, mut r1) = (modulus.to_vec(), a);
    let (mut s0, mut s1): (Vec<Int>, Vec<Int>) = (Vec::new(), alloc::vec![Int::ONE]);
    while !r1.is_empty() {
        let (q, r) = poly_divmod(&r0, &r1, p);
        let s = poly_sub(&s0, &poly_mul(&q, &s1, p), p);
        r0 = r1;
        r1 = r;
        s0 = s1;
        s1 = s;
    }
    // gcd = r0 must be a nonzero constant (guaranteed when `modulus` is
    // irreducible and `a ≠ 0`); otherwise `a` shares a factor and has no inverse.
    if deg_of(&r0) != Some(0) {
        return None;
    }
    let inv_c = r0[0].modinv(p)?;
    Some(poly_rem(&poly_scale(&s0, &inv_c, p), modulus, p))
}

/// Increments `digits` as a little-endian base-`p` counter; returns `false` on
/// wrap-around (every digit was `p−1`).
fn incr_base_p(digits: &mut [Int], p: &Int) -> bool {
    for d in digits.iter_mut() {
        *d = d.add(&Int::ONE);
        if *d < *p {
            return true;
        }
        *d = Int::ZERO;
    }
    false
}

// ===========================================================================
// Public types.
// ===========================================================================

/// The shared field context of `GF(pᵏ)`: the prime `p`, the degree `k`, and the
/// monic irreducible modulus polynomial.
struct FieldData {
    /// The characteristic (a prime).
    p: Int,
    /// The extension degree.
    k: usize,
    /// The modulus `f`, low-to-high, length `k+1`, monic, each entry in `[0, p)`.
    modulus: Vec<Int>,
    /// Precomputed reduction table `red[j] = x^{k+j} mod f` for `j = 0 … k−2`
    /// (i.e. `x^k … x^{2k−2}`), each a length-`k` coefficient vector in `[0, p)`.
    /// Empty for `k == 1` (products never reach degree `k`). Turns the post-
    /// multiply reduction into a fixed set of multiply-accumulates.
    red: Vec<Vec<Int>>,
}

impl FieldData {
    /// Assembles the shared field context, precomputing the `x^k … x^{2k-2} mod f`
    /// reduction table. `modulus` is monic of degree `k`, low-to-high, each entry
    /// already in `[0, p)`.
    fn build(p: Int, k: usize, modulus: Vec<Int>) -> FieldData {
        // red[j] = x^{k+j} mod f. Start at x^k mod f, then multiply by x (a shift)
        // and reduce the single degree-k overflow term with red[0], iterating.
        let mut red: Vec<Vec<Int>> = Vec::new();
        if k >= 2 {
            let mut xk = alloc::vec![Int::ZERO; k + 1];
            xk[k] = Int::ONE;
            let r0 = pad_to(poly_rem(&xk, &modulus, &p), k);
            red.push(r0.clone());
            let mut cur = r0;
            for _ in 1..(k - 1) {
                // Multiply `cur` (degree < k) by x: shift up by one. The new
                // coefficient at index k is `cur[k-1]`; fold it back via red[0].
                let top = cur[k - 1].clone();
                let mut next = alloc::vec![Int::ZERO; k];
                for i in (1..k).rev() {
                    next[i] = cur[i - 1].clone();
                }
                if !top.is_zero() {
                    for i in 0..k {
                        next[i].addmul(&top, &red[0][i]);
                    }
                }
                for c in next.iter_mut() {
                    *c = c.rem_euclid(&p);
                }
                red.push(next.clone());
                cur = next;
            }
        }
        FieldData { p, k, modulus, red }
    }

    /// Product of two canonical coefficient vectors (each length `k`, entries in
    /// `[0, p)`), reduced modulo the field modulus `f`. Returns the canonical
    /// residue as a length-`k` vector in `[0, p)` — bit-identical to
    /// `poly_rem(poly_mul(a, b, p), modulus, p)`.
    fn mul_coeffs(&self, a: &[Int], b: &[Int]) -> Vec<Int> {
        let k = self.k;
        let p = &self.p;
        if k == 1 {
            return alloc::vec![a[0].mul(&b[0]).rem_euclid(p)];
        }
        // Raw schoolbook product (no per-coefficient reduction): degree ≤ 2k−2.
        // Every term is non-negative, so the accumulators stay non-negative and
        // ride Int's inline i128 fast path for word-size primes.
        let mut prod = alloc::vec![Int::ZERO; 2 * k - 1];
        for (i, ai) in a.iter().enumerate() {
            if ai.is_zero() {
                continue;
            }
            for (j, bj) in b.iter().enumerate() {
                prod[i + j].addmul(ai, bj);
            }
        }
        self.reduce_raw(&prod)
    }

    /// Reduces a raw (unreduced) product polynomial `prod` — coefficients
    /// low-to-high, length `2k−1`, each a *non-negative* integer that may exceed
    /// `p` — modulo the field modulus `f` and modulo `p`, returning the canonical
    /// length-`k` residue in `[0, p)`.
    ///
    /// Bit-identical to reducing `prod` coefficient-wise mod `p` and then
    /// `poly_rem`-ing by `f`: it folds the high terms `x^k … x^{2k−2}` down
    /// through the precomputed reduction table, accumulating into the low part,
    /// then reduces every coefficient mod `p`. Because reduction mod `(p, f)` is a
    /// ring homomorphism, summing raw products and reducing once yields the same
    /// canonical residue as reducing each product and summing.
    fn reduce_raw(&self, prod: &[Int]) -> Vec<Int> {
        let k = self.k;
        let p = &self.p;
        if k == 1 {
            return alloc::vec![prod[0].rem_euclid(p)];
        }
        // Fold the high coefficients x^k … x^{2k−2} down through the reduction
        // table, accumulating into the low part.
        let mut acc: Vec<Int> = prod[..k].to_vec();
        for (j, red_j) in self.red.iter().enumerate() {
            let hi = prod[k + j].rem_euclid(p);
            if hi.is_zero() {
                continue;
            }
            for (ai, rj) in acc.iter_mut().zip(red_j.iter()) {
                ai.addmul(&hi, rj);
            }
        }
        for c in acc.iter_mut() {
            *c = c.rem_euclid(p);
        }
        acc
    }
}

/// A finite field `GF(pᵏ) = GF(p)[x] / (f)` for a monic irreducible `f`.
///
/// The field data (prime, degree, modulus) is reference-counted, so
/// [`GfElement`]s built from it share the context cheaply.
#[derive(Clone)]
pub struct GaloisField {
    data: Rc<FieldData>,
}

/// An element of a [`GaloisField`]: the residue class of a polynomial, stored as
/// its canonical representative of degree `< k`.
#[derive(Clone)]
pub struct GfElement {
    field: Rc<FieldData>,
    /// Coefficients low-to-high, length exactly `k`, each in `[0, p)`.
    coeffs: Vec<Int>,
}

impl GaloisField {
    /// Builds `GF(pᵏ)` from an explicit modulus.
    ///
    /// `p` must be prime and `modulus` a **monic** polynomial (low-to-high
    /// coefficients, leading coefficient `1` after reduction mod `p`) of degree
    /// `k ≥ 1` that is **irreducible** over `GF(p)`. Returns `None` if `p` is not
    /// prime, if the modulus is not monic of positive degree, or if it is
    /// reducible.
    pub fn new(p: Int, modulus: &[Int]) -> Option<GaloisField> {
        if !p.is_prime_bpsw() {
            return None;
        }
        let reduced: Vec<Int> = modulus.iter().map(|c| c.rem_euclid(&p)).collect();
        if reduced.len() < 2 {
            return None; // degree must be ≥ 1
        }
        let k = reduced.len() - 1;
        if !reduced[k].is_one() {
            return None; // not monic (leading coefficient must be 1)
        }
        if !is_irreducible(&reduced, &p) {
            return None;
        }
        Some(GaloisField {
            data: Rc::new(FieldData::build(p, k, reduced)),
        })
    }

    /// Builds `GF(pᵏ)`, automatically finding a monic irreducible modulus of
    /// degree `k` over `GF(p)`.
    ///
    /// Candidate monic polynomials `xᵏ + c_{k-1}xᵏ⁻¹ + … + c₀` are enumerated
    /// deterministically (the low coefficients as a base-`p` counter) and each is
    /// tested with Rabin's irreducibility test; the first irreducible one is
    /// used. Returns `None` if `p` is not prime or `k == 0`.
    pub fn create(p: Int, k: usize) -> Option<GaloisField> {
        if k == 0 || !p.is_prime_bpsw() {
            return None;
        }
        let mut low = alloc::vec![Int::ZERO; k];
        loop {
            let mut modulus = low.clone();
            modulus.push(Int::ONE); // monic degree-k candidate
            if is_irreducible(&modulus, &p) {
                return Some(GaloisField {
                    data: Rc::new(FieldData::build(p, k, modulus)),
                });
            }
            if !incr_base_p(&mut low, &p) {
                return None; // exhausted (unreachable: irreducibles always exist)
            }
        }
    }

    /// The characteristic `p`.
    #[inline]
    pub fn characteristic(&self) -> Int {
        self.data.p.clone()
    }

    /// The extension degree `k`.
    #[inline]
    pub fn degree(&self) -> usize {
        self.data.k
    }

    /// The field order `pᵏ` (the number of elements).
    #[inline]
    pub fn order(&self) -> Int {
        self.data.p.pow(self.data.k as u32)
    }

    /// The modulus polynomial coefficients, low-to-high (length `k+1`, monic).
    #[inline]
    pub fn modulus(&self) -> &[Int] {
        &self.data.modulus
    }

    /// Builds an element from polynomial coefficients (low-to-high); the input is
    /// reduced modulo `p` and modulo the field modulus, so any length is accepted.
    pub fn element(&self, coeffs: &[Int]) -> GfElement {
        let reduced: Vec<Int> = coeffs.iter().map(|c| c.rem_euclid(&self.data.p)).collect();
        let residue = poly_rem(&reduced, &self.data.modulus, &self.data.p);
        GfElement {
            field: self.data.clone(),
            coeffs: pad_to(residue, self.data.k),
        }
    }

    /// The additive identity `0`.
    #[inline]
    pub fn zero(&self) -> GfElement {
        GfElement {
            field: self.data.clone(),
            coeffs: alloc::vec![Int::ZERO; self.data.k],
        }
    }

    /// The multiplicative identity `1`.
    pub fn one(&self) -> GfElement {
        self.from_int(&Int::ONE)
    }

    /// The generator `a`, the residue class of `x` (`[0, 1, 0, …]`).
    pub fn generator(&self) -> GfElement {
        self.element(&[Int::ZERO, Int::ONE])
    }

    /// The constant element embedding a `GF(p)` value `c` (`c mod p`).
    pub fn from_int(&self, c: &Int) -> GfElement {
        let mut coeffs = alloc::vec![Int::ZERO; self.data.k];
        coeffs[0] = c.rem_euclid(&self.data.p);
        GfElement {
            field: self.data.clone(),
            coeffs,
        }
    }

    /// The Frobenius endomorphism `elem ↦ elemᵖ`. It fixes `GF(p)`, satisfies
    /// `(a+b)ᵖ = aᵖ + bᵖ`, and `elem^(pᵏ) = elem` for every element. Panics if
    /// `elem` belongs to a different field.
    pub fn frobenius(&self, elem: &GfElement) -> GfElement {
        assert!(
            same_field(&self.data, &elem.field),
            "GaloisField::frobenius: element from a different field"
        );
        elem.pow(&self.data.p)
    }

    /// A primitive element: a generator of the (cyclic) multiplicative group
    /// `GF(pᵏ)*`, i.e. an element of multiplicative order `pᵏ − 1`.
    ///
    /// Candidates are enumerated deterministically and each is checked with the
    /// order test `g^{(pᵏ−1)/q} ≠ 1` for every prime `q | pᵏ − 1` (reusing
    /// integer factorization). Such an element always exists.
    pub fn primitive_element(&self) -> GfElement {
        let one = self.one();
        let n = self.order().sub(&Int::ONE); // pᵏ − 1
        // Distinct prime divisors of the group order.
        let mut primes: Vec<Int> = Vec::new();
        for q in n.factorize() {
            if primes.last() != Some(&q) {
                primes.push(q);
            }
        }
        let mut digits = alloc::vec![Int::ZERO; self.data.k];
        while incr_base_p(&mut digits, &self.data.p) {
            let cand = self.element(&digits);
            if cand.is_zero() {
                continue;
            }
            let is_primitive = primes.iter().all(|q| cand.pow(&n.div_exact(q)) != one);
            if is_primitive {
                return cand;
            }
        }
        unreachable!("a primitive element always exists in a finite field");
    }
}

/// Whether two field contexts are the same field (pointer or value equality).
fn same_field(a: &Rc<FieldData>, b: &Rc<FieldData>) -> bool {
    Rc::ptr_eq(a, b) || (a.p == b.p && a.modulus == b.modulus)
}

impl GfElement {
    /// The field this element belongs to.
    pub fn field(&self) -> GaloisField {
        GaloisField {
            data: self.field.clone(),
        }
    }

    /// The canonical residue coefficients, low-to-high (length `k`, each in
    /// `[0, p)`).
    #[inline]
    pub fn to_coefficients(&self) -> &[Int] {
        &self.coeffs
    }

    /// Returns `true` if this is the zero element.
    pub fn is_zero(&self) -> bool {
        self.coeffs.iter().all(Int::is_zero)
    }

    /// Returns `true` if this is the multiplicative identity `1`.
    pub fn is_one(&self) -> bool {
        self.coeffs[0].is_one() && self.coeffs[1..].iter().all(Int::is_zero)
    }

    fn same_field(&self, other: &GfElement) {
        assert!(
            same_field(&self.field, &other.field),
            "GfElement: operands from different fields"
        );
    }

    #[inline]
    fn wrap(&self, coeffs: Vec<Int>) -> GfElement {
        GfElement {
            field: self.field.clone(),
            coeffs,
        }
    }

    /// Returns `self + rhs`.
    pub fn add(&self, rhs: &GfElement) -> GfElement {
        self.same_field(rhs);
        let p = &self.field.p;
        let coeffs = (0..self.field.k)
            .map(|i| self.coeffs[i].add(&rhs.coeffs[i]).rem_euclid(p))
            .collect();
        self.wrap(coeffs)
    }

    /// Returns `self − rhs`.
    pub fn sub(&self, rhs: &GfElement) -> GfElement {
        self.same_field(rhs);
        let p = &self.field.p;
        let coeffs = (0..self.field.k)
            .map(|i| self.coeffs[i].sub(&rhs.coeffs[i]).rem_euclid(p))
            .collect();
        self.wrap(coeffs)
    }

    /// Returns `−self`.
    pub fn neg(&self) -> GfElement {
        let p = &self.field.p;
        let coeffs = self.coeffs.iter().map(|c| c.neg().rem_euclid(p)).collect();
        self.wrap(coeffs)
    }

    /// Returns `self · rhs`.
    pub fn mul(&self, rhs: &GfElement) -> GfElement {
        self.same_field(rhs);
        let coeffs = self.field.mul_coeffs(&self.coeffs, &rhs.coeffs);
        self.wrap(coeffs)
    }

    /// Returns the multiplicative inverse `self⁻¹`, or `None` if `self` is zero.
    ///
    /// Computed by the extended Euclidean algorithm in `GF(p)[x]`; cross-checks
    /// against `self^(pᵏ−2)`.
    pub fn inv(&self) -> Option<GfElement> {
        let p = &self.field.p;
        let u = poly_inv(&self.coeffs, &self.field.modulus, p)?;
        Some(self.wrap(pad_to(u, self.field.k)))
    }

    /// Returns `self / rhs = self · rhs⁻¹`. Panics if `rhs` is zero.
    pub fn div(&self, rhs: &GfElement) -> GfElement {
        self.mul(
            &rhs.inv()
                .expect("GfElement::div: divisor is not invertible"),
        )
    }

    /// Returns `self` raised to `exp` (negative exponents invert first).
    pub fn pow(&self, exp: &Int) -> GfElement {
        if exp.is_negative() {
            return self
                .inv()
                .expect("GfElement::pow: base not invertible for a negative exponent")
                .pow(&exp.abs());
        }
        // Square-and-multiply, low bit first.
        let mut result = GfElement {
            field: self.field.clone(),
            coeffs: {
                let mut c = alloc::vec![Int::ZERO; self.field.k];
                c[0] = Int::ONE;
                c
            },
        };
        let mut base = self.clone();
        for i in 0..exp.bit_len() {
            if exp.bit(i) {
                result = result.mul(&base);
            }
            base = base.mul(&base);
        }
        result
    }
}

impl PartialEq for GfElement {
    fn eq(&self, other: &Self) -> bool {
        same_field(&self.field, &other.field) && self.coeffs == other.coeffs
    }
}

impl Eq for GfElement {}

impl fmt::Display for GfElement {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        if self.is_zero() {
            return f.write_str("0");
        }
        let mut first = true;
        for (i, c) in self.coeffs.iter().enumerate().rev() {
            if c.is_zero() {
                continue;
            }
            if !first {
                f.write_str(" + ")?;
            }
            first = false;
            match i {
                0 => write!(f, "{c}")?,
                1 if c.is_one() => f.write_str("a")?,
                1 => write!(f, "{c}·a")?,
                _ if c.is_one() => write!(f, "a^{i}")?,
                _ => write!(f, "{c}·a^{i}")?,
            }
        }
        Ok(())
    }
}

impl fmt::Debug for GfElement {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "GfElement({} in GF({}^{}))",
            self, self.field.p, self.field.k
        )
    }
}

impl fmt::Debug for GaloisField {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "GaloisField(GF({}^{}))", self.data.p, self.data.k)
    }
}

macro_rules! gf_binop {
    ($tr:ident, $m:ident, $atr:ident, $am:ident) => {
        impl core::ops::$tr for GfElement {
            type Output = GfElement;
            #[inline]
            fn $m(self, rhs: GfElement) -> GfElement {
                GfElement::$m(&self, &rhs)
            }
        }
        impl core::ops::$tr<&GfElement> for GfElement {
            type Output = GfElement;
            #[inline]
            fn $m(self, rhs: &GfElement) -> GfElement {
                GfElement::$m(&self, rhs)
            }
        }
        impl core::ops::$tr<GfElement> for &GfElement {
            type Output = GfElement;
            #[inline]
            fn $m(self, rhs: GfElement) -> GfElement {
                GfElement::$m(self, &rhs)
            }
        }
        impl core::ops::$tr<&GfElement> for &GfElement {
            type Output = GfElement;
            #[inline]
            fn $m(self, rhs: &GfElement) -> GfElement {
                GfElement::$m(self, rhs)
            }
        }
        impl core::ops::$atr<GfElement> for GfElement {
            #[inline]
            fn $am(&mut self, rhs: GfElement) {
                *self = GfElement::$m(self, &rhs);
            }
        }
        impl core::ops::$atr<&GfElement> for GfElement {
            #[inline]
            fn $am(&mut self, rhs: &GfElement) {
                *self = GfElement::$m(self, rhs);
            }
        }
    };
}

gf_binop!(Add, add, AddAssign, add_assign);
gf_binop!(Sub, sub, SubAssign, sub_assign);
gf_binop!(Mul, mul, MulAssign, mul_assign);
gf_binop!(Div, div, DivAssign, div_assign);

impl core::ops::Neg for GfElement {
    type Output = GfElement;
    #[inline]
    fn neg(self) -> GfElement {
        GfElement::neg(&self)
    }
}
impl core::ops::Neg for &GfElement {
    type Output = GfElement;
    #[inline]
    fn neg(self) -> GfElement {
        GfElement::neg(self)
    }
}

// ===========================================================================
// Nested (2-D) Kronecker-substitution multiplication for `Poly<GfElement>`.
//
// A `Poly<GfElement>` is an *outer* polynomial `A(y) = Σ aᵢ yⁱ` whose
// coefficients `aᵢ ∈ GF(pᵏ)` are themselves *inner* polynomials `aᵢ(t) =
// Σ_{s<k} aᵢ,ₛ tˢ` over `GF(p)` (residues in `[0, p)`). Treating both operands as
// bivariate integer polynomials `A(y, t)`, `B(y, t)` (coefficients in `[0, p)`),
// their product `P(y, t) = A·B` has
//
//     P_{l,s} = Σ_{i+j=l, s'+s''=s} aᵢ,ₛ' · bⱼ,ₛ''      (an integer),
//
// and the field product is `cₗ = (Σ P_{l,·}) reduced mod (f, p)`. We recover the
// integer `P_{l,s}` with a single big-integer multiply by a two-level Kronecker
// substitution `t = 2^{k_in}`, `y = 2^{k_out}`, then reduce each outer
// coefficient back into `GF(pᵏ)`.
//
// No-overlap (exactness) — the subtle part, proven at BOTH levels:
//
//   Inner slot `k_in`. A single output coefficient `P_{l,s}` is a sum of at most
//   `C = min(nₐ, n_b)·k` products of two residues `< p` (≤ `min(nₐ,n_b)` pairs
//   `i+j=l`, ≤ `k` pairs `s'+s''=s`), so `P_{l,s} ≤ C·(p−1)²`. With `bp =
//   bit_len(p)` we have `(p−1)² < 2^{2bp}`, and `clog_in = bit_len(C)` gives
//   `C < 2^{clog_in}`, hence `P_{l,s} < 2^{2bp+clog_in} = 2^{k_in}`. So every
//   inner slot holds its coefficient with no carry into its neighbour — and,
//   crucially, `k_in` already accounts for the *outer* accumulation (the `min`
//   factor), because the single big multiply sums all `i+j=l` inner blocks into
//   one outer slot.
//
//   Outer slot `k_out`. The packed inner block of the l-th outer coefficient is
//   `P̂ₗ = Σ_{s=0}^{2k-2} P_{l,s} 2^{s·k_in} ≤ 2^{(2k-1)·k_in} − 1 <
//   2^{(2k-1)·k_in}`. Choosing `k_out = 2k·k_in + clog_out ≥ (2k-1)·k_in` (with
//   `clog_out = bit_len(min(nₐ,n_b))` a harmless margin) keeps `P̂ₗ < 2^{k_out}`,
//   so the outer slots never overlap either.
//
// Reduction mod `(f, p)` is a ring homomorphism `GF(p)[t] → GF(pᵏ)`, so reducing
// the summed raw products once (`reduce_raw`) yields the same *canonical* residue
// as the Karatsuba path (reduce each product, then sum). The recovered
// `Vec<GfElement>` is therefore bit-identical.
// ===========================================================================

/// Smaller-operand coefficient count at or above which `Poly<GfElement>::mul`
/// routes through nested Kronecker substitution. Below it the generic Karatsuba
/// path wins (and stays the differential reference). Chosen from benchmarking
/// (see `examples/gf_kronecker_bench.rs`).
#[cfg(feature = "poly")]
const GF_KRONECKER_THRESHOLD: usize = 24;

/// Nested-Kronecker multiply hook for `Poly<GfElement>` (see
/// [`Ring::poly_mul`](crate::ring::Ring::poly_mul)). Returns `None` below the
/// degree threshold (and `Some(empty)` for an empty operand) so small products
/// stay on the Karatsuba path.
#[cfg(feature = "poly")]
pub(crate) fn gf_kronecker_mul(a: &[GfElement], b: &[GfElement]) -> Option<Vec<GfElement>> {
    if a.is_empty() || b.is_empty() {
        return Some(Vec::new());
    }
    if a.len().min(b.len()) < GF_KRONECKER_THRESHOLD {
        return None;
    }
    Some(gf_kronecker_convolve(a, b))
}

/// Core nested-Kronecker product of two nonempty `GfElement` coefficient slices
/// over one field `GF(pᵏ)`. Returns the product coefficients (length
/// `nₐ + n_b − 1`), each a canonical field residue — exactly the Karatsuba
/// convolution. See the module-level no-overlap argument above.
#[cfg(feature = "poly")]
pub(crate) fn gf_kronecker_convolve(a: &[GfElement], b: &[GfElement]) -> Vec<GfElement> {
    use crate::poly::{kronecker_pack, kronecker_unpack};

    let field = a[0].field.clone();
    let k = field.k;
    let bp = field.p.bit_len();
    let na = a.len();
    let nb = b.len();
    let min_len = na.min(nb);
    let num_outer = na + nb - 1;

    // Inner slot: 2·bit_len(p) plus room for the ≤ min_len·k accumulated terms.
    let inner_count = (min_len as u64) * (k as u64);
    let clog_in = 64 - inner_count.leading_zeros(); // = bit_len(inner_count) ≥ 1
    let k_in = 2 * bp + clog_in;

    // Outer slot: holds the (2k−1)-slot inner block of each outer coefficient.
    let clog_out = 64 - (min_len as u64).leading_zeros(); // = bit_len(min_len)
    let k_out_u64 = 2u64 * (k as u64) * u64::from(k_in) + u64::from(clog_out);
    let k_out =
        u32::try_from(k_out_u64).expect("gf_kronecker_convolve: outer slot exceeds 2³² bits");

    // Pack each element into a k·k_in-bit inner block, then pack the outer
    // polynomial at stride k_out; one big multiply carries out both convolutions.
    let blocks_a: Vec<Int> = a.iter().map(|e| kronecker_pack(&e.coeffs, k_in)).collect();
    let blocks_b: Vec<Int> = b.iter().map(|e| kronecker_pack(&e.coeffs, k_in)).collect();
    let prod = kronecker_pack(&blocks_a, k_out).mul(&kronecker_pack(&blocks_b, k_out));

    // Unpack outer slots, then the 2k−1 inner slots of each, and reduce mod (f,p).
    let inner_slots = 2 * k - 1;
    kronecker_unpack(&prod, k_out, num_outer)
        .into_iter()
        .map(|block| {
            let raw = kronecker_unpack(&block, k_in, inner_slots);
            GfElement {
                field: field.clone(),
                coeffs: field.reduce_raw(&raw),
            }
        })
        .collect()
}