ocas-poly 0.26.0

//! Multi-modular Gröbner basis computation.
//!
//! Computes Gröbner bases of ℚ-ideals by:
//!
//! 1. clearing denominators to a primitive ℤ ideal;
//! 2. computing F5 bases of the mod-p images for several primes
//!    ([`crate::groebner::f5::f5`], which takes the native i64 fast path
//!    for `p < 2³¹`) in parallel, keeping only "lucky" images whose
//!    leading-monomial set matches the current reference;
//! 3. reconstructing rational coefficients by CRT + rational
//!    reconstruction ([`crate::rational_reconstruction`]);
//! 4. verifying the candidate exactly over ℚ (S-pair criterion plus
//!    membership of the original generators);
//! 5. falling back to a p-adic Hensel lift ([`hensel_lift_groebner`])
//!    once enough primes have accumulated, and ultimately to the exact
//!    ℚ F5 algorithm.
//!
//! The driver mirrors the multi-prime pattern of the modular GCD
//! ([`crate::multivariate_gcd::modular_gcd_x`]): lucky-prime filtering by
//! degree comparison there becomes leading-monomial-ideal inclusion here,
//! and the reconstruction/verification loop is identical in spirit.
//!
//! All primes are drawn from [2³⁰, 2³¹) so the ℤ_p fast paths apply.

use num_bigint::BigInt;
use ocas_core::FastHashMap as HashMap;
use ocas_domain::number_theory::{crt, primes_from, symmetric_mod};
use ocas_domain::{Domain, FiniteField, Integer, IntegerDomain, Rational, RationalDomain};
use rayon::prelude::*;
use smallvec::SmallVec;

use super::f4::{echelonize_fp, norm_mod};
use super::f5::f5;
use super::GroebnerBasis;
use crate::multivariate_gcd::qmpoly_to_primitive_zmpoly;
use crate::rational_reconstruction::rational_reconstruction;
use crate::sparse::{MonomialOrder, SparseMultivariatePolynomial, monomial_divides, monomial_lcm};

/// Primes in [2³⁰, 2³¹): large enough to reconstruct small rational
/// coefficients from a single image, small enough for the ℤ_p F4/F5 fast
/// paths (which require `p < 2³¹`).
fn gb_primes() -> impl Iterator<Item = i64> {
    primes_from(&Integer::from(1i64 << 30))
        .filter(|p| *p < Integer::from(1i64 << 31))
        .map(|p| p.to_i64().expect("prime below 2^31 fits i64"))
}

/// Reduce a primitive ℤ ideal modulo `p`.
///
/// Returns `None` when some leading coefficient vanishes mod `p`: such a
/// prime is "bad" — the mod-p image would change the leading term and
/// cannot be used for reconstruction.
fn ideal_mod_p<O: MonomialOrder>(
    ideal_z: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    p: i64,
) -> Option<Vec<SparseMultivariatePolynomial<FiniteField, O>>> {
    let field = FiniteField::new(BigInt::from(p));
    let p_int = Integer::from(p);
    let mut out = Vec::with_capacity(ideal_z.len());
    for poly in ideal_z {
        if let Some(lc) = poly.leading_coeff()
            && lc.mod_floor(&p_int).is_zero()
        {
            return None;
        }
        let terms: Vec<(Vec<usize>, _)> = poly
            .terms_ref()
            .iter()
            .map(|(e, c)| (e.to_vec(), field.element(c.to_bigint())))
            .collect();
        out.push(SparseMultivariatePolynomial::from_terms(
            field.clone(),
            poly.n_vars(),
            terms,
        ));
    }
    Some(out)
}

/// The sorted leading-monomial set of a Gröbner basis.
fn lm_set<D: Domain, O: MonomialOrder>(gb: &GroebnerBasis<D, O>) -> Vec<SmallVec<[usize; 4]>> {
    let mut lms: Vec<SmallVec<[usize; 4]>> = gb
        .basis
        .iter()
        .filter_map(|p| p.leading_monomial().cloned())
        .collect();
    lms.sort();
    lms
}

/// True iff the monomial ideal generated by `a` is contained in that of
/// `b`: every monomial of `a` is divisible by some monomial of `b`.
///
/// Note: `monomial_divides(x, y)` is true iff `y` divides `x`.
fn lm_subset(a: &[SmallVec<[usize; 4]>], b: &[SmallVec<[usize; 4]>]) -> bool {
    a.iter()
        .all(|ma| b.iter().any(|mb| monomial_divides(ma, mb)))
}

/// Convert a sparse integer polynomial to a sparse rational polynomial.
fn zmpoly_to_qmpoly<O: MonomialOrder>(
    f: &SparseMultivariatePolynomial<IntegerDomain, O>,
) -> SparseMultivariatePolynomial<RationalDomain, O> {
    SparseMultivariatePolynomial::from_terms(
        RationalDomain,
        f.n_vars(),
        f.terms_ref()
            .iter()
            .map(|(e, c)| (e.to_vec(), Rational::from_integer(c.clone())))
            .collect(),
    )
}

// =========================================================================
//  CRT accumulation + rational reconstruction
// =========================================================================

/// Incremental CRT state per basis element: monomial → (residue, modulus).
///
/// Each accepted prime image folds in its coefficient residues; a
/// monomial absent from an image contributes a 0 residue (its coefficient
/// is 0 mod that prime).
/// Monomial → (accumulated residue, accumulated modulus) for one basis element.
type CrtCoeffs = HashMap<SmallVec<[usize; 4]>, (Integer, Integer)>;

struct CrtAcc {
    /// Per basis element: monomial → (accumulated residue, accumulated modulus).
    coeffs: Vec<CrtCoeffs>,
}

impl CrtAcc {
    fn new(n_elements: usize) -> Self {
        Self {
            coeffs: (0..n_elements).map(|_| HashMap::default()).collect(),
        }
    }

    /// Fold one prime image into the accumulator.
    fn fold<O: MonomialOrder>(&mut self, p: i64, gb: &GroebnerBasis<FiniteField, O>) {
        let p_int = Integer::from(p);
        for (i, poly) in gb.basis.iter().enumerate() {
            let map = &mut self.coeffs[i];
            let mut seen: Vec<SmallVec<[usize; 4]>> = Vec::new();
            for (exp, c) in poly.terms_ref() {
                seen.push(exp.clone());
                let residue = Integer::from(c.value().clone());
                let entry = map
                    .entry(exp.clone())
                    .or_insert_with(|| (Integer::from(0), Integer::from(1)));
                if let Some((r, m)) = crt(&entry.0, &entry.1, &residue, &p_int) {
                    *entry = (r, m);
                }
            }
            // Monomials known from earlier images but absent here get a
            // 0 residue for this prime.
            let known: Vec<SmallVec<[usize; 4]>> = map.keys().cloned().collect();
            for exp in known {
                if !seen.contains(&exp) {
                    let entry = map.get_mut(&exp).unwrap();
                    if let Some((r, m)) = crt(&entry.0, &entry.1, &Integer::from(0), &p_int) {
                        *entry = (r, m);
                    }
                }
            }
        }
    }
}

/// Attempt to assemble a rational candidate basis from the accumulated
/// CRT state. Returns `None` when some coefficient cannot be rationally
/// reconstructed with the current modulus.
fn reconstruct_q<O: MonomialOrder>(
    acc: &CrtAcc,
    n_vars: usize,
) -> Option<Vec<SparseMultivariatePolynomial<RationalDomain, O>>> {
    let mut basis: Vec<SparseMultivariatePolynomial<RationalDomain, O>> =
        Vec::with_capacity(acc.coeffs.len());
    for map in &acc.coeffs {
        let mut terms: Vec<(Vec<usize>, Rational)> = Vec::with_capacity(map.len());
        for (exp, (r, m)) in map {
            let s = symmetric_mod(r, m);
            let (n, d) = rational_reconstruction(&s, m)?;
            terms.push((
                exp.to_vec(),
                Rational::from_bigints(n.to_bigint(), d.to_bigint()),
            ));
        }
        basis.push(SparseMultivariatePolynomial::from_terms(
            RationalDomain,
            n_vars,
            terms,
        ));
    }
    Some(basis)
}

/// Exact verification of a candidate ℚ basis against the primitive ℤ
/// ideal: the candidate must be a Gröbner basis (all S-polynomials reduce
/// to zero) and every generator must reduce to zero against it.
fn verify_candidate<O: MonomialOrder>(
    ideal_z: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    cand: &[SparseMultivariatePolynomial<RationalDomain, O>],
) -> bool {
    if cand.is_empty() {
        return ideal_z.iter().all(|p| p.is_zero());
    }
    let gb = GroebnerBasis {
        basis: cand.to_vec(),
    }
    .minimize()
    .auto_reduce();
    if !gb.is_groebner_basis() {
        return false;
    }
    for f in ideal_z {
        let fq = zmpoly_to_qmpoly(f);
        if !fq.reduce(&gb.basis).is_zero() {
            return false;
        }
    }
    true
}

// =========================================================================
//  Driver
// =========================================================================

/// Compute a Gröbner basis of a ℚ-ideal by the multi-modular method.
///
/// Parallel F5 over batches of lucky primes, CRT + rational
/// reconstruction, exact ℚ verification, a Hensel-lift shortcut once
/// enough primes have accumulated, and a final fallback to the exact ℚ
/// F5 algorithm. The function is total: it always returns a valid
/// reduced Gröbner basis.
#[allow(clippy::too_many_lines)]
pub fn groebner_basis_multi_modular<O: MonomialOrder + Send + Sync>(
    ideal: &[SparseMultivariatePolynomial<RationalDomain, O>],
) -> GroebnerBasis<RationalDomain, O> {
    // 1. Clear denominators: the ideal as primitive integer polynomials.
    let f_z: Vec<SparseMultivariatePolynomial<IntegerDomain, O>> = ideal
        .iter()
        .filter(|p| !p.is_zero())
        .map(qmpoly_to_primitive_zmpoly)
        .collect();
    if f_z.is_empty() {
        return GroebnerBasis { basis: vec![] };
    }
    let n_vars = f_z[0].n_vars();

    let mut primes = gb_primes();

    // 2. The first valid image establishes the reference leading-monomial set.
    type LmRef = Vec<SmallVec<[usize; 4]>>;
    type Accepted<O2> = Vec<(i64, GroebnerBasis<FiniteField, O2>)>;
    let (mut lm_ref, mut accepted): (LmRef, Accepted<O>) = loop {
        let Some(p) = primes.next() else {
            return f5(ideal);
        };
        if let Some(gens) = ideal_mod_p(&f_z, p) {
            let gb = f5(&gens);
            break (lm_set(&gb), vec![(p, gb)]);
        }
    };
    let mut acc = CrtAcc::new(lm_ref.len());
    acc.fold(accepted[0].0, &accepted[0].1);

    let batch_size = rayon::current_num_threads().max(1);
    let mut hensel_tried = false;

    loop {
        // 3. Batch of primes, F5 images computed in parallel.
        let batch: Vec<i64> = primes.by_ref().take(batch_size).collect();
        if batch.is_empty() {
            break;
        }
        let results: Vec<(i64, Option<GroebnerBasis<FiniteField, O>>)> = batch
            .par_iter()
            .map(|&p| (p, ideal_mod_p(&f_z, p).map(|gens| f5(&gens))))
            .collect();

        // 4. Lucky-prime filter by leading-monomial ideal inclusion.
        for (p, img) in results {
            let Some(gb) = img else {
                continue; // bad prime (leading coefficient vanished)
            };
            let lms = lm_set(&gb);
            if lms == lm_ref {
                accepted.push((p, gb.clone()));
                acc.fold(p, &gb);
            } else if lm_subset(&lms, &lm_ref) {
                // Strictly smaller monomial ideal: closer to the true
                // basis — restart accumulation around this image.
                lm_ref = lms;
                accepted.clear();
                accepted.push((p, gb.clone()));
                acc = CrtAcc::new(lm_ref.len());
                acc.fold(p, &gb);
            }
            // Strictly larger or incomparable: discard the image.
        }

        // 5. CRT + rational reconstruction after each batch.
        if let Some(cand) = reconstruct_q(&acc, n_vars)
            && verify_candidate(&f_z, &cand)
        {
            return GroebnerBasis { basis: cand }.minimize().auto_reduce();
        }

        // 7. Hensel shortcut: once enough primes have accumulated, lift the
        // first lucky image p-adically; on failure continue with CRT.
        if !hensel_tried && accepted.len() > 16 {
            hensel_tried = true;
            let (p1, gb1) = &accepted[0];
            if let Some(lifted) = hensel_lift_groebner(&f_z, &gb1.basis, *p1, 64) {
                let lifted_q: Vec<SparseMultivariatePolynomial<RationalDomain, O>> =
                    lifted.iter().map(zmpoly_to_qmpoly).collect();
                if verify_candidate(&f_z, &lifted_q) {
                    return GroebnerBasis { basis: lifted_q }.minimize().auto_reduce();
                }
            }
        }
        if accepted.len() > 64 {
            break;
        }
    }

    // Fallback: exact Gröbner basis over ℚ.
    f5(ideal)
}

/// Multi-modular Gröbner basis for any domain, gated by an `Any` check on
/// the coefficients.
///
/// Returns `None` when the coefficients are not [`Rational`] (e.g. finite
/// fields), in which case the caller should fall back to its usual path.
/// For ℚ ideals this always returns `Some`.
pub(crate) fn groebner_basis_mm<D: Domain + 'static, O: MonomialOrder + Send + Sync>(
    ideal: &[SparseMultivariatePolynomial<D, O>],
) -> Option<GroebnerBasis<D, O>> {
    if ideal.is_empty() {
        return Some(GroebnerBasis { basis: vec![] });
    }
    let mut q_ideal: Vec<SparseMultivariatePolynomial<RationalDomain, O>> =
        Vec::with_capacity(ideal.len());
    for poly in ideal {
        let mut terms = Vec::with_capacity(poly.n_terms());
        for (exp, c) in poly.terms_ref() {
            let r = (c as &dyn std::any::Any)
                .downcast_ref::<Rational>()?
                .clone();
            terms.push((exp.to_vec(), r));
        }
        q_ideal.push(SparseMultivariatePolynomial::from_terms(
            RationalDomain,
            poly.n_vars(),
            terms,
        ));
    }
    let gb_q = groebner_basis_multi_modular(&q_ideal);
    let domain = ideal[0].domain().clone();
    let mut basis: Vec<SparseMultivariatePolynomial<D, O>> =
        Vec::with_capacity(gb_q.basis.len());
    for poly in &gb_q.basis {
        let mut terms = Vec::with_capacity(poly.n_terms());
        for (exp, c) in poly.terms_ref() {
            let boxed: Box<dyn std::any::Any> = Box::new(c.clone());
            let elem = *boxed.downcast::<D::Element>().ok()?;
            terms.push((exp.to_vec(), elem));
        }
        basis.push(SparseMultivariatePolynomial::from_terms(
            domain.clone(),
            poly.n_vars(),
            terms,
        ));
    }
    Some(GroebnerBasis { basis })
}

// =========================================================================
//  Hensel lifting (trace-free linear p-adic lift)
// =========================================================================

/// Lift the mod-p reduced Gröbner basis `g0` of the primitive ℤ ideal
/// `ideal_z` to a ℚ Gröbner basis by trace-free linear p-adic lifting.
///
/// `g0` must be the reduced monic Gröbner basis of `ideal_mod_p(ideal_z, p)`
/// (i.e. the output of F5 over 𝔽_p). The support of the basis is fixed to
/// that of `g0` (a lucky-prime assumption). Returns the primitive integer
/// polynomials of the verified ℚ basis, or `None` when `p` is unsuitable
/// (a residual fails the required congruence) or `max_digits` digits are
/// exhausted.
///
/// Algorithm (each step adds one p-adic digit):
/// 1. Compute all residuals mod `p^(k+1)`: the normal forms of the
///    generators and of every S-polynomial against the lifted basis. All
///    residual coefficients must be ≡ 0 (mod `p^k`), otherwise `p` is
///    rejected.
/// 2. Linearize: the residual is affine in the unknown tail-coefficient
///    corrections; each column is obtained by perturbing one unknown by
///    `p^k` and re-computing all residuals.
/// 3. Solve the resulting 𝔽_p system with [`echelonize_fp`] (free
///    variables set to zero); no solution → `None`.
/// 4. `G_{k+1} = G_k + p^k·γ`, symmetrically normalized mod `p^(k+1)`.
/// 5. From the second step on, try convergence: rational-reconstruct
///    every coefficient; when the candidate matches the previous step's
///    and passes the ℚ verification, return it.
///
/// All division is exact because the basis stays monic (corrections touch
/// only tail coefficients), so no modular inverses are needed — this is
/// what makes the lift "trace-free".
pub fn hensel_lift_groebner<O: MonomialOrder>(
    ideal_z: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    g0: &[SparseMultivariatePolynomial<FiniteField, O>],
    p: i64,
    max_digits: usize,
) -> Option<Vec<SparseMultivariatePolynomial<IntegerDomain, O>>> {
    let n_vars = g0[0].n_vars();
    let p_int = Integer::from(p);

    // Lift g0 to an integer basis with coefficients in [0, p).
    let mut g: Vec<SparseMultivariatePolynomial<IntegerDomain, O>> = g0
        .iter()
        .map(|poly| {
            let terms: Vec<(Vec<usize>, Integer)> = poly
                .terms_ref()
                .iter()
                .map(|(e, c)| (e.to_vec(), Integer::from(c.value().clone())))
                .collect();
            SparseMultivariatePolynomial::from_terms(IntegerDomain, n_vars, terms)
        })
        .collect();

    // Unknowns: the tail coefficients of each basis element.
    let mut unknowns: Vec<(usize, SmallVec<[usize; 4]>)> = Vec::new();
    for (i, poly) in g0.iter().enumerate() {
        let lm = poly.leading_monomial()?.clone();
        for exp in poly.terms_ref().keys() {
            if *exp != lm {
                unknowns.push((i, exp.clone()));
            }
        }
    }
    let n_unknowns = unknowns.len();

    let mut prev_candidate: Option<Vec<SparseMultivariatePolynomial<RationalDomain, O>>> = None;
    let mut pk = p_int.clone(); // p^k at the start of iteration k (k = 1: p^1)
    for k in 1..=max_digits {
        let modulus = &pk * &p_int; // p^(k+1)

        // (a) Baseline residuals mod p^(k+1); every coefficient must be
        // divisible by p^k or the previous lift step did not converge.
        let rho0 = residuals(&g, ideal_z, &modulus)?;
        for r in &rho0 {
            for c in r.terms_ref().values() {
                if !c.mod_floor(&pk).is_zero() {
                    return None;
                }
            }
        }

        // (b) Jacobian columns: perturb each unknown by p^k and recompute
        // all residuals; the difference / p^k is the column (mod p).
        let mut residuals_all: Vec<Vec<SparseMultivariatePolynomial<IntegerDomain, O>>> =
            Vec::with_capacity(1 + n_unknowns);
        residuals_all.push(rho0);
        for (gi, m) in &unknowns {
            let mut g_pert = g.clone();
            let mut coeff = g_pert[*gi].coeff(m);
            coeff += &pk;
            g_pert[*gi].set_term_external(m.to_vec(), coeff);
            let rho_e = residuals(&g_pert, ideal_z, &modulus)?;
            for r in &rho_e {
                for c in r.terms_ref().values() {
                    if !c.mod_floor(&pk).is_zero() {
                        return None;
                    }
                }
            }
            residuals_all.push(rho_e);
        }

        // Equation rows: (residual index within one residuals() result,
        // monomial). Every residuals() call returns the same number of
        // polynomials in the same order, so the inner index is shared
        // across all perturbed runs.
        let n_residuals = residuals_all[0].len();
        let mut row_map: HashMap<(usize, SmallVec<[usize; 4]>), usize> = HashMap::default();
        let mut matrix: Vec<Vec<(i64, usize)>> = Vec::new();
        for res in &residuals_all {
            for (ri, r) in res.iter().enumerate() {
                debug_assert_eq!(res.len(), n_residuals);
                for exp in r.terms_ref().keys() {
                    row_map.entry((ri, exp.clone())).or_insert_with(|| {
                        matrix.push(Vec::new());
                        matrix.len() - 1
                    });
                }
            }
        }

        for (col, _) in unknowns.iter().enumerate() {
            let rho_e = &residuals_all[1 + col];
            let rho0 = &residuals_all[0];
            for (&(ri, ref exp), &row) in &row_map {
                let ce = rho_e[ri].coeff(exp);
                let c0 = rho0[ri].coeff(exp);
                let diff = &ce - &c0;
                let v = (diff / &pk).mod_floor(&p_int).to_i64().unwrap();
                if v != 0 {
                    matrix[row].push((v, col));
                }
            }
        }

        // RHS: −rho0 / p^k (mod p).
        let rho0 = &residuals_all[0];
        for (&(ri, ref exp), &row) in &row_map {
            let c0 = rho0[ri].coeff(exp);
            let v = (c0 / &pk).mod_floor(&p_int).to_i64().unwrap();
            if v != 0 {
                matrix[row].push((p_int.to_i64().unwrap() - v, n_unknowns));
            }
        }

        // (c) Solve over 𝔽_p.
        let x = solve_fp(&mut matrix, n_unknowns, p)?;

        // (d) G_{k+1} = G_k + p^k·γ, symmetrically normalized mod p^(k+1).
        for (col, (gi, m)) in unknowns.iter().enumerate() {
            let gamma = x[col];
            if gamma == 0 {
                continue;
            }
            let delta = &pk * &Integer::from(gamma);
            let mut coeff = g[*gi].coeff(m);
            coeff += &delta;
            coeff = symmetric_mod(&coeff, &modulus);
            g[*gi].set_term_external(m.to_vec(), coeff);
        }
        for poly in &mut g {
            let terms: Vec<(Vec<usize>, Integer)> = poly
                .terms_ref()
                .iter()
                .map(|(e, c)| (e.to_vec(), symmetric_mod(c, &modulus)))
                .collect();
            *poly = SparseMultivariatePolynomial::from_terms(IntegerDomain, n_vars, terms);
        }

        // (e) Convergence: rational-reconstruct every coefficient; when the
        // candidate matches the previous step's, verify exactly over ℚ.
        if k >= 2 {
            match rr_candidate(&g, &modulus) {
                Some(cand) => {
                    if prev_candidate.as_ref() == Some(&cand)
                        && verify_candidate(ideal_z, &cand)
                    {
                        return Some(cand.iter().map(qmpoly_to_primitive_zmpoly).collect());
                    }
                    prev_candidate = Some(cand);
                }
                None => prev_candidate = None,
            }
        }

        pk = modulus;
    }
    None
}

/// Rational-reconstruct every coefficient of the lifted basis.
fn rr_candidate<O: MonomialOrder>(
    g: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    modulus: &Integer,
) -> Option<Vec<SparseMultivariatePolynomial<RationalDomain, O>>> {
    let mut out = Vec::with_capacity(g.len());
    for poly in g {
        let mut terms: Vec<(Vec<usize>, Rational)> = Vec::with_capacity(poly.n_terms());
        for (exp, c) in poly.terms_ref() {
            let s = symmetric_mod(c, modulus);
            let (n, d) = rational_reconstruction(&s, modulus)?;
            terms.push((
                exp.to_vec(),
                Rational::from_bigints(n.to_bigint(), d.to_bigint()),
            ));
        }
        out.push(SparseMultivariatePolynomial::from_terms(
            RationalDomain,
            poly.n_vars(),
            terms,
        ));
    }
    Some(out)
}

/// All residuals: normal forms of the generators and of every S-polynomial
/// of the lifted basis, computed mod `modulus` with coefficients in [0, m).
fn residuals<O: MonomialOrder>(
    g: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    ideal_z: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    modulus: &Integer,
) -> Option<Vec<SparseMultivariatePolynomial<IntegerDomain, O>>> {
    let mut out = Vec::with_capacity(ideal_z.len() + g.len() * (g.len() - 1) / 2);
    for f in ideal_z {
        out.push(reduce_mod_ring(f, g, modulus));
    }
    for i in 0..g.len() {
        for j in (i + 1)..g.len() {
            let s = spoly_mod_ring(&g[i], &g[j], modulus)?;
            out.push(reduce_mod_ring(&s, g, modulus));
        }
    }
    Some(out)
}

/// S-polynomial of two monic polynomials over ℤ/modulus; the leading
/// terms cancel exactly, so the result has degree < max.
fn spoly_mod_ring<O: MonomialOrder>(
    a: &SparseMultivariatePolynomial<IntegerDomain, O>,
    b: &SparseMultivariatePolynomial<IntegerDomain, O>,
    modulus: &Integer,
) -> Option<SparseMultivariatePolynomial<IntegerDomain, O>> {
    let lm_a = a.leading_monomial()?;
    let lm_b = b.leading_monomial()?;
    let lcm = monomial_lcm(lm_a, lm_b);
    let diff_a: SmallVec<[usize; 4]> = lcm.iter().zip(lm_a.iter()).map(|(x, y)| x - y).collect();
    let diff_b: SmallVec<[usize; 4]> = lcm.iter().zip(lm_b.iter()).map(|(x, y)| x - y).collect();
    let sa = a.mul_monomial(&diff_a);
    let sb = b.mul_monomial(&diff_b);
    let mut merged: HashMap<SmallVec<[usize; 4]>, Integer> = HashMap::default();
    for (e, c) in sa.terms_ref() {
        if *e == lcm {
            continue;
        }
        let v = c.mod_floor(modulus);
        let entry = merged.entry(e.clone()).or_insert_with(|| Integer::from(0));
        *entry = (entry.clone() + &v).mod_floor(modulus);
    }
    for (e, c) in sb.terms_ref() {
        if *e == lcm {
            continue;
        }
        let v = (Integer::from(0) - c.mod_floor(modulus)).mod_floor(modulus);
        let entry = merged.entry(e.clone()).or_insert_with(|| Integer::from(0));
        *entry = (entry.clone() + &v).mod_floor(modulus);
    }
    let terms: Vec<(Vec<usize>, Integer)> = merged
        .into_iter()
        .filter(|(_, v)| !v.is_zero())
        .map(|(e, v)| (e.to_vec(), v))
        .collect();
    Some(SparseMultivariatePolynomial::from_terms(
        IntegerDomain,
        a.n_vars(),
        terms,
    ))
}

/// Full reduction of `poly` against the monic basis `basis` over
/// ℤ/modulus. Always reduces the largest reducible monomial first (this
/// makes the sequence of processed monomials strictly decreasing in the
/// monomial order, hence terminating); the reducer is the first basis
/// element whose leading monomial divides it, keeping the process
/// deterministic across perturbed runs.
fn reduce_mod_ring<O: MonomialOrder>(
    poly: &SparseMultivariatePolynomial<IntegerDomain, O>,
    basis: &[SparseMultivariatePolynomial<IntegerDomain, O>],
    modulus: &Integer,
) -> SparseMultivariatePolynomial<IntegerDomain, O> {
    let order = poly.order.clone();
    let n_vars = poly.n_vars();
    let mut coeffs: HashMap<SmallVec<[usize; 4]>, Integer> = HashMap::default();
    for (exp, c) in poly.terms_ref() {
        let v = c.mod_floor(modulus);
        if !v.is_zero() {
            coeffs.insert(exp.clone(), v);
        }
    }
    let basis_lms: Vec<(usize, SmallVec<[usize; 4]>)> = basis
        .iter()
        .enumerate()
        .filter_map(|(i, b)| b.leading_monomial().map(|lm| (i, lm.clone())))
        .collect();
    loop {
        // Largest monomial (per order) divisible by some basis leading
        // monomial.
        type Best = (SmallVec<[usize; 4]>, Integer, usize, SmallVec<[usize; 4]>);
        let mut best: Option<Best> = None;
        for (exp, c) in &coeffs {
            if let Some((bi, lm)) = basis_lms
                .iter()
                .find(|(_, lm)| monomial_divides(exp, lm))
            {
                match &best {
                    Some((be, _, _, _)) if order.cmp(exp, be) != std::cmp::Ordering::Greater => {}
                    _ => best = Some((exp.clone(), c.clone(), *bi, lm.clone())),
                }
            }
        }
        let Some((exp, c, bi, lm)) = best else {
            break;
        };
        let diff: SmallVec<[usize; 4]> = exp.iter().zip(lm.iter()).map(|(a, b)| a - b).collect();
        let mut updates: Vec<(SmallVec<[usize; 4]>, Integer)> = Vec::new();
        for (texp, tc) in basis[bi].terms_ref() {
            let mut mexp: SmallVec<[usize; 4]> = SmallVec::with_capacity(n_vars);
            for v in 0..n_vars {
                mexp.push(diff.get(v).copied().unwrap_or(0) + texp.get(v).copied().unwrap_or(0));
            }
            updates.push((mexp, (&c * tc).mod_floor(modulus)));
        }
        for (mexp, sub) in updates {
            let new_v = match coeffs.get(&mexp) {
                Some(e) => (e - &sub).mod_floor(modulus),
                None => (Integer::from(0) - &sub).mod_floor(modulus),
            };
            if new_v.is_zero() {
                coeffs.remove(&mexp);
            } else {
                coeffs.insert(mexp, new_v);
            }
        }
    }
    let terms: Vec<(Vec<usize>, Integer)> = coeffs
        .into_iter()
        .map(|(e, v)| (e.to_vec(), v))
        .collect();
    SparseMultivariatePolynomial::from_terms(IntegerDomain, n_vars, terms)
}

/// Solve `A·x = b` over 𝔽_p for the augmented sparse system in
/// `matrix` (unknown columns `0..n_unknowns`, RHS in the last column).
///
/// Rows must have columns in ascending order. Free variables are set to
/// zero. Returns `None` when the system is inconsistent.
fn solve_fp(
    matrix: &mut Vec<Vec<(i64, usize)>>,
    n_unknowns: usize,
    p: i64,
) -> Option<Vec<i64>> {
    let ncols = n_unknowns + 1;
    let mut pivots: Vec<Option<usize>> = Vec::new();
    echelonize_fp(matrix, ncols, p, &mut pivots);

    let mut pivot_row_of_col: Vec<Option<usize>> = vec![None; ncols];
    for (r, row) in matrix.iter().enumerate() {
        if let Some(&(_, col)) = row.first() {
            pivot_row_of_col[col] = Some(r);
        }
    }
    // A pivot at the RHS column is an inconsistent equation (0 = nonzero).
    if pivot_row_of_col[n_unknowns].is_some() {
        return None;
    }

    // Gauss-Jordan completion: eliminate each pivot column from every
    // other row so each pivot row carries only its pivot and the RHS.
    let mut scratch: Vec<(i64, usize)> = Vec::new();
    #[allow(clippy::needless_range_loop)]
    for col in 0..n_unknowns {
        let Some(pr) = pivot_row_of_col[col] else {
            continue;
        };
        let pivot = matrix[pr].clone();
        #[allow(clippy::needless_range_loop)]
        for r in 0..matrix.len() {
            if r == pr {
                continue;
            }
            if let Some(pos) = matrix[r].iter().position(|&(_, c)| c == col) {
                let c = matrix[r][pos].0;
                sub_scaled_any_fp(&mut matrix[r], &pivot, c, p, &mut scratch);
            }
        }
    }

    let mut x = vec![0i64; n_unknowns];
    for (col, &pr) in pivot_row_of_col.iter().enumerate().take(n_unknowns) {
        if let Some(pr) = pr {
            let row = &matrix[pr];
            let rhs = row
                .iter()
                .find(|&&(_, c)| c == n_unknowns)
                .map(|&(v, _)| v)
                .unwrap_or(0);
            x[col] = rhs;
        }
    }
    Some(x)
}

/// Sparse `row -= c * pivot` (mod `p`) by merging the two column-ascending
/// rows. General merge: no assumption that the head columns coincide.
fn sub_scaled_any_fp(
    row: &mut Vec<(i64, usize)>,
    pivot: &[(i64, usize)],
    c: i64,
    p: i64,
    scratch: &mut Vec<(i64, usize)>,
) {
    scratch.clear();
    scratch.reserve(row.len() + pivot.len());
    let mut i = 0;
    let mut j = 0;
    while i < row.len() && j < pivot.len() {
        let (rc, rcol) = row[i];
        let (pc, pcol) = pivot[j];
        if rcol < pcol {
            scratch.push((rc, rcol));
            i += 1;
        } else if rcol > pcol {
            let v = norm_mod(-c * pc, p);
            if v != 0 {
                scratch.push((v, pcol));
            }
            j += 1;
        } else {
            let v = norm_mod(rc - c * pc, p);
            if v != 0 {
                scratch.push((v, rcol));
            }
            i += 1;
            j += 1;
        }
    }
    scratch.extend_from_slice(&row[i..]);
    for &(pc, pcol) in &pivot[j..] {
        let v = norm_mod(-c * pc, p);
        if v != 0 {
            scratch.push((v, pcol));
        }
    }
    std::mem::swap(row, scratch);
}

// =========================================================================
//  Tests
// =========================================================================

#[cfg(test)]
mod tests {
    use super::*;
    use crate::groebner::f5::f5;
    use crate::sparse::Lex;

    fn r(n: i64, d: i64) -> Rational {
        Rational::new(n, d)
    }

    fn qpoly(
        terms: Vec<(Vec<usize>, Rational)>,
        n_vars: usize,
    ) -> SparseMultivariatePolynomial<RationalDomain, Lex> {
        SparseMultivariatePolynomial::from_terms(RationalDomain, n_vars, terms)
    }

    fn zpoly(
        terms: Vec<(Vec<usize>, i64)>,
        n_vars: usize,
    ) -> SparseMultivariatePolynomial<IntegerDomain, Lex> {
        SparseMultivariatePolynomial::from_terms(
            IntegerDomain,
            n_vars,
            terms
                .into_iter()
                .map(|(e, c)| (e, Integer::from(c)))
                .collect(),
        )
    }

    fn cyclic_q(n: usize) -> Vec<SparseMultivariatePolynomial<RationalDomain, Lex>> {
        let mut gens = Vec::with_capacity(n);
        for k in 1..n {
            let mut terms = Vec::new();
            for start in 0..n {
                let mut exps = vec![0usize; n];
                for j in 0..k {
                    exps[(start + j) % n] = 1;
                }
                terms.push((exps, r(1, 1)));
            }
            gens.push(qpoly(terms, n));
        }
        let full_exps = vec![1usize; n];
        gens.push(qpoly(
            vec![(full_exps, r(1, 1)), (vec![0usize; n], r(-1, 1))],
            n,
        ));
        gens
    }

    #[test]
    fn crt_rr_reconstruction() {
        // The rational 3/7: 3·7^{-1} mod 11 = 3·8 = 24 ≡ 2; mod 13: 3·2 = 6.
        let m1 = Integer::from(11);
        let r1 = Integer::from(2);
        let m2 = Integer::from(13);
        let r2 = Integer::from(6);
        let (r, m) = crt(&r1, &m1, &r2, &m2).unwrap();
        let s = symmetric_mod(&r, &m);
        let (n, d) = rational_reconstruction(&s, &m).unwrap();
        assert_eq!(n, Integer::from(3));
        assert_eq!(d, Integer::from(7));
    }

    #[test]
    fn bad_prime_skipped() {
        // 2x - 1: the leading coefficient 2 vanishes mod 2.
        let ideal_z = vec![zpoly(vec![(vec![1], 2), (vec![0], -1)], 1)];
        assert!(ideal_mod_p(&ideal_z, 2).is_none());
        assert!(ideal_mod_p(&ideal_z, 3).is_some());
        // The driver still computes the correct basis (x - 1/2).
        let ideal_q = vec![qpoly(vec![(vec![1], r(1, 1)), (vec![0], r(-1, 2))], 1)];
        assert_eq!(
            groebner_basis_multi_modular(&ideal_q),
            f5(&ideal_q),
            "multi-modular must skip bad primes"
        );
    }

    #[test]
    fn mm_matches_f5_cyclic() {
        for n in 2..=4 {
            let ideal = cyclic_q(n);
            assert_eq!(
                groebner_basis_multi_modular(&ideal),
                f5(&ideal),
                "multi-modular != f5 for cyclic-{n}"
            );
        }
    }

    #[test]
    fn mm_matches_f5_rational_coeffs() {
        // GB {x + 1/2, y - x} has a genuine rational coefficient.
        let ideal = vec![
            qpoly(vec![(vec![1, 0], r(1, 1)), (vec![0, 0], r(1, 2))], 2),
            qpoly(vec![(vec![0, 1], r(1, 1)), (vec![1, 0], r(-1, 1))], 2),
        ];
        assert_eq!(
            groebner_basis_multi_modular(&ideal),
            f5(&ideal),
            "multi-modular must reconstruct rational coefficients"
        );
    }

    #[test]
    fn hensel_lift_small_ideal() {
        // {x^2 - 2, y - x} over ℚ: reduced GB is the ideal itself.
        let ideal_z = vec![
            zpoly(vec![(vec![2, 0], 1), (vec![0, 0], -2)], 2),
            zpoly(vec![(vec![0, 1], 1), (vec![1, 0], -1)], 2),
        ];
        let p = gb_primes().next().unwrap();
        let img = ideal_mod_p(&ideal_z, p).unwrap();
        let g0 = f5(&img);
        let lifted =
            hensel_lift_groebner(&ideal_z, &g0.basis, p, 64).expect("lift should succeed");
        let lifted_q: Vec<SparseMultivariatePolynomial<RationalDomain, Lex>> =
            lifted.iter().map(zmpoly_to_qmpoly).collect();
        let expect = f5(&[
            qpoly(vec![(vec![2, 0], r(1, 1)), (vec![0, 0], r(-2, 1))], 2),
            qpoly(vec![(vec![0, 1], r(1, 1)), (vec![1, 0], r(-1, 1))], 2),
        ]);
        assert_eq!(
            GroebnerBasis {
                basis: lifted_q
            }
            .minimize()
            .auto_reduce(),
            expect,
            "lifted basis must equal the direct ℚ basis"
        );
    }

    #[test]
    fn hensel_lift_rational_coeff() {
        // 2x + 1 has GB {x + 1/2}: the lift must recover 1/2 p-adically.
        let ideal_z = vec![zpoly(vec![(vec![1], 2), (vec![0], 1)], 1)];
        let p = gb_primes().next().unwrap();
        let img = ideal_mod_p(&ideal_z, p).unwrap();
        let g0 = f5(&img);
        let lifted =
            hensel_lift_groebner(&ideal_z, &g0.basis, p, 64).expect("lift should succeed");
        assert_eq!(lifted, ideal_z);
    }

    #[test]
    fn hensel_lift_bad_prime_none() {
        // A basis that is NOT the mod-p image of the true GB: the residual
        // of 2x - 1 against {x} is -1, not divisible by p → None.
        let ideal_z = vec![zpoly(vec![(vec![1], 2), (vec![0], -1)], 1)];
        let p = gb_primes().next().unwrap();
        let field = FiniteField::new(BigInt::from(p));
        let wrong_g0 = vec![SparseMultivariatePolynomial::from_terms(
            field.clone(),
            1,
            vec![(vec![1], field.element(1))],
        )];
        assert!(hensel_lift_groebner(&ideal_z, &wrong_g0, p, 64).is_none());
    }
}