alkahest-cas 3.8.0

//! Primary decomposition and radicals over ℚ\[x₁,…,xₙ\] (Gianni–Trager–Zacharias
//! fragment).
//!
//! # What is certified, and what is refused
//!
//! There is no general primary-decomposition algorithm here, and none of the
//! routines pretend otherwise. Every component this module returns is primary
//! for a reason it can state, and every radical it returns is radical for a
//! reason it can state:
//!
//! * **Monomial ideals** decompose completely, by the coprime split
//!   `⟨J, u·v⟩ = ⟨J, u⟩ ∩ ⟨J, v⟩` for coprime monomials `u, v`, into irreducible
//!   monomial ideals `⟨x_{j₁}^{a₁}, …, x_{j_k}^{a_k}⟩`, each primary with
//!   associated prime `⟨x_{j₁}, …, x_{j_k}⟩`. Their radical is generated by the
//!   square-free part of each generating monomial.
//! * **Principal ideals** decompose as `⟨∏ pᵢ^{eᵢ}⟩ = ∩ ⟨pᵢ^{eᵢ}⟩`, since
//!   ℚ\[x₁,…,xₙ\] is a UFD and `⟨pᵉ⟩` is `⟨p⟩`-primary for a prime element `p`.
//!   Their radical is `⟨∏ pᵢ⟩`.
//! * **Zero-dimensional ideals** have `√I = I + ⟨sqfree(pᵢ) : i⟩` where
//!   `pᵢ` generates `I ∩ ℚ[xᵢ]` (Seidenberg; ℚ is perfect), and a
//!   shape-position basis `⟨x₀ − g₀(t), …, x_{n−2} − g_{n−2}(t), h(t)⟩` with
//!   `h = c·pᵉ` and `p` irreducible is primary, because its quotient is the
//!   local Artinian ring `ℚ[t]/⟨pᵉ⟩`.
//! * **Splits** — `I = (I : x_i^∞) ∩ (I + ⟨x_i⟩)` when the intersection checks
//!   out, and `I = ∩_j (I + ⟨p_j^{e_j}⟩)` by CRT when `I` contains a univariate
//!   with several distinct irreducible factors — reduce an ideal to smaller
//!   ones, which must themselves land in one of the certified classes.
//!
//! Anything else **refuses**: see [`IdealRefusal`]. Returning `I` itself as
//! though it were its own radical, or as though it were primary, is the failure
//! this module is written to avoid — a caller who reads a field named
//! `associated_prime` is entitled to assume it names a prime.

use crate::errors::AlkahestError;
use crate::flint::mpoly::{FlintMPoly, FlintMPolyCtx, FlintMPolyFactor};
use crate::flint::FlintPoly;
use crate::poly::groebner::ideal::GbPoly;
use crate::poly::groebner::monomial_order::MonomialOrder;
use crate::poly::groebner::{is_zero_dimensional, GroebnerBasis};
use std::cell::RefCell;
use std::collections::BTreeMap;
use std::fmt;
use std::sync::Arc;

const MAX_SPLIT_DEPTH: usize = 48;

/// Ceiling on the number of irreducible components a monomial split may produce
/// before the routine gives up and refuses. The split is exponential in the
/// worst case (`⟨x₁y₁, …, x_ky_k⟩` has `2^k` components) and this gate runs on
/// every pull request.
const MAX_MONOMIAL_COMPONENTS: usize = 256;

/// One primary component together with its associated prime (√Q).
///
/// Both fields carry a guarantee: `primary` is primary and `associated_prime`
/// is prime, for one of the reasons listed in the module documentation. A
/// component that cannot be certified is never returned — the whole call
/// refuses instead.
#[derive(Clone, Debug)]
pub struct PrimaryComponent {
    pub primary: GroebnerBasis,
    pub associated_prime: GroebnerBasis,
}

/// Primary-decomposition failures (inconsistent input, depth limit, FLINT).
///
/// `Factorization` doubles as the carrier for [`IdealRefusal`] — see that type
/// for why the refusal cannot be a variant of its own.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum PrimaryDecompositionError {
    EmptyGenerators,
    InconsistentNvars,
    RecursionDepth,
    Factorization(&'static str),
}

impl fmt::Display for PrimaryDecompositionError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            PrimaryDecompositionError::EmptyGenerators => {
                write!(f, "ideal generators must be non-empty")
            }
            PrimaryDecompositionError::InconsistentNvars => {
                write!(f, "inconsistent n_vars across generators")
            }
            PrimaryDecompositionError::RecursionDepth => {
                write!(f, "primary decomposition exceeded recursion depth")
            }
            // Reworded from "univariate factorization failed: {msg}": this variant
            // is also the carrier for `IdealRefusal`, whose message is a full
            // sentence of its own.  Every in-module construction site now spells
            // out what failed.
            PrimaryDecompositionError::Factorization(msg) => write!(f, "{msg}"),
        }
    }
}

impl std::error::Error for PrimaryDecompositionError {}

impl AlkahestError for PrimaryDecompositionError {
    fn code(&self) -> &'static str {
        match self {
            PrimaryDecompositionError::EmptyGenerators => "E-IDEAL-001",
            PrimaryDecompositionError::InconsistentNvars => "E-IDEAL-002",
            PrimaryDecompositionError::RecursionDepth => "E-IDEAL-003",
            PrimaryDecompositionError::Factorization(_) => "E-IDEAL-004",
        }
    }

    fn remediation(&self) -> Option<&'static str> {
        match self {
            PrimaryDecompositionError::EmptyGenerators => Some("pass at least one generator"),
            PrimaryDecompositionError::InconsistentNvars => {
                Some("all generators must be polynomials in the same variable list")
            }
            PrimaryDecompositionError::RecursionDepth => Some(
                "the saturation split recursed past its depth limit; simplify the \
                 generating set",
            ),
            PrimaryDecompositionError::Factorization(_) => {
                Some("report the generating set as a minimal failing example")
            }
        }
    }
}

// ---------------------------------------------------------------------------
// Refusals, reported out of band
// ---------------------------------------------------------------------------

/// Which routine declined, which fixes the stable code the refusal carries.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum RefusalSite {
    /// `√I` for an ideal outside the certified classes — `E-IDEAL-005`.
    Radical,
    /// A primary decomposition with a component that cannot be certified
    /// primary — `E-IDEAL-006`.
    Decomposition,
}

/// An ideal-theoretic question this module cannot answer, with the code it
/// carries.
///
/// # Why this is not an error variant
///
/// [`PrimaryDecompositionError`] is a public *exhaustive* enum, so growing it a
/// `NotCertifiable` variant is a major semver break — and so is marking it
/// `#[non_exhaustive]` to allow one later. A correctness fix inside a patch
/// release cannot spend a major version, so the refusal travels out of band:
/// the refusing routine returns `PrimaryDecompositionError::Factorization` with
/// this type's message, and the real code is recorded here for
/// [`take_ideal_refusal`] to hand to the bindings.
///
/// This is the pattern [`crate::matrix::take_zero_test_refusal`] already uses
/// for undecided zero tests inside `LinearAlgebraError::UnsupportedField`, and
/// [`crate::calculus::limits::last_budget_trip`] for budget trips inside
/// `LimitError::DepthExceeded`.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct IdealRefusal {
    site: RefusalSite,
    message: &'static str,
}

impl fmt::Display for IdealRefusal {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str(self.message)
    }
}

impl std::error::Error for IdealRefusal {}

impl AlkahestError for IdealRefusal {
    fn code(&self) -> &'static str {
        match self.site {
            RefusalSite::Radical => "E-IDEAL-005",
            RefusalSite::Decomposition => "E-IDEAL-006",
        }
    }

    fn remediation(&self) -> Option<&'static str> {
        match self.site {
            RefusalSite::Radical => Some(
                "radical is certified for monomial, principal and zero-dimensional \
                 ideals; intersect the associated primes of a primary decomposition \
                 if one is available",
            ),
            RefusalSite::Decomposition => Some(
                "primary decomposition is certified for monomial and principal ideals, \
                 for saturation/CRT splits of them, and for shape-position \
                 zero-dimensional ideals; no general algorithm is implemented",
            ),
        }
    }
}

thread_local! {
    /// The refusal behind the `Factorization` error the current thread is about
    /// to return, when that variant is a carrier rather than what it usually
    /// means.
    static LAST_IDEAL_REFUSAL: RefCell<Option<IdealRefusal>> = const { RefCell::new(None) };
}

/// Drop any recorded refusal, so a genuine FLINT failure can never be
/// re-attributed to a refusal left behind by an earlier call on this thread.
fn forget_ideal_refusal() {
    LAST_IDEAL_REFUSAL.with(|c| *c.borrow_mut() = None);
}

fn refuse(site: RefusalSite, message: &'static str) -> PrimaryDecompositionError {
    LAST_IDEAL_REFUSAL.with(|c| *c.borrow_mut() = Some(IdealRefusal { site, message }));
    PrimaryDecompositionError::Factorization(message)
}

/// Take the refusal behind the error that just came back, if there was one.
///
/// Bindings call this when [`radical`] or [`primary_decomposition`] returns
/// `PrimaryDecompositionError::Factorization` and raise the refusal's own
/// `E-IDEAL-005` / `E-IDEAL-006` when it is present, so the caller still gets
/// the specific code. `Some` means *this* error is a refusal; `None` means the
/// variant means what it usually means — a FLINT failure.
///
/// Consuming, so one refusal is reported once and cannot leak into a later
/// unrelated error. Thread-local.
pub fn take_ideal_refusal() -> Option<IdealRefusal> {
    LAST_IDEAL_REFUSAL.with(|c| c.borrow_mut().take())
}

const RADICAL_NOT_CERTIFIABLE: &str = "cannot certify √I for this ideal: it is neither \
     monomial nor principal nor zero-dimensional, and no primary decomposition of it \
     could be certified either. The general radical needs Gianni–Trager–Zacharias or a \
     characteristic-set method, which is not implemented — refusing rather than \
     returning I unchanged as though it were already radical";

const DECOMPOSITION_NOT_CERTIFIABLE: &str = "cannot certify a primary decomposition of \
     this ideal: it splits no further and is not one of the classes whose primarity can \
     be established (monomial, principal, or a shape-position zero-dimensional basis). \
     Refusing rather than reporting the ideal itself as a primary component with an \
     unjustified associated prime";

const MONOMIAL_TOO_LARGE: &str = "the irreducible decomposition of this monomial ideal \
     exceeded the component ceiling; refusing rather than returning a truncated \
     intersection that is not equal to the input ideal";

const FACTORIZATION_FAILED: &str = "univariate factorization failed: FLINT could not \
     factor a generator over ℤ";

const MULTIVARIATE_FACTORIZATION_FAILED: &str = "multivariate factorization failed: \
     FLINT could not factor the generator of a principal ideal over ℤ";

fn lcm_rational_denoms(coeffs: &[rug::Rational]) -> rug::Integer {
    let mut m = rug::Integer::from(1);
    for c in coeffs {
        m = m.lcm(c.denom());
    }
    m
}

/// Clear denominators and strip integer content → primitive `FlintPoly`.
fn primitive_flint_from_rational_asc(coeffs: &[rug::Rational]) -> Option<FlintPoly> {
    if coeffs.is_empty() {
        return None;
    }
    let mut hi = coeffs.len();
    while hi > 0 && coeffs[hi - 1] == 0 {
        hi -= 1;
    }
    if hi == 0 {
        return None;
    }
    let coeffs = &coeffs[..hi];
    let lcm = lcm_rational_denoms(coeffs);
    let mut ints: Vec<rug::Integer> = Vec::with_capacity(coeffs.len());
    for c in coeffs {
        let t = c * rug::Rational::from((lcm.clone(), 1));
        let (n, d) = t.into_numer_denom();
        debug_assert_eq!(d, rug::Integer::from(1));
        ints.push(n);
    }
    let mut g = ints[0].clone();
    for a in ints.iter().skip(1) {
        g = g.gcd(a);
    }
    if g != 0 {
        for a in &mut ints {
            let (q, r) = a.clone().div_rem(g.clone());
            debug_assert_eq!(r, 0);
            *a = q;
        }
    }
    Some(FlintPoly::from_rug_coefficients(&ints))
}

/// Radical √I, or a refusal.
///
/// Certified for monomial, principal and zero-dimensional ideals, and for any
/// ideal whose [`primary_decomposition`] can be certified (√I is then the
/// intersection of the associated primes). Everything else raises
/// `E-IDEAL-005`; see [`IdealRefusal`].
///
/// Before 3.8 this returned the input basis unchanged whenever no univariate
/// generator was available, which asserted `√I = I` with nothing behind it:
/// `radical([(x−y)²])` did not contain `x−y`, although
/// `√⟨(x−y)²⟩ = ⟨x−y⟩` exactly.
pub fn radical(
    gens: Vec<GbPoly>,
    order: MonomialOrder,
) -> Result<GroebnerBasis, PrimaryDecompositionError> {
    validate_gens(&gens)?;
    forget_ideal_refusal();
    let gb = GroebnerBasis::compute(gens, order);
    radical_from_basis(&gb, order)
        .ok_or_else(|| refuse(RefusalSite::Radical, RADICAL_NOT_CERTIFIABLE))
}

/// Irredundant primary decomposition, or a refusal.
///
/// Every returned component is primary and carries a prime `associated_prime`,
/// for one of the reasons in the module documentation; the intersection of the
/// components is the input ideal. When some component cannot be certified the
/// whole call raises `E-IDEAL-006` rather than reporting an unjustified one.
pub fn primary_decomposition(
    gens: Vec<GbPoly>,
    order: MonomialOrder,
) -> Result<Vec<PrimaryComponent>, PrimaryDecompositionError> {
    validate_gens(&gens)?;
    forget_ideal_refusal();
    let gb = GroebnerBasis::compute(gens, order);
    if is_unit_ideal(&gb) {
        return Ok(vec![]);
    }
    let mut raw = decompose_recursive(gb, order, 0)?;
    dedup_components(&mut raw);
    drop_redundant_components(&mut raw);
    Ok(raw)
}

fn validate_gens(gens: &[GbPoly]) -> Result<(), PrimaryDecompositionError> {
    if gens.is_empty() {
        return Err(PrimaryDecompositionError::EmptyGenerators);
    }
    let n = gens[0].n_vars;
    if gens.iter().any(|g| g.n_vars != n) {
        return Err(PrimaryDecompositionError::InconsistentNvars);
    }
    Ok(())
}

fn is_unit_ideal(gb: &GroebnerBasis) -> bool {
    gb.generators().iter().any(|g| {
        g.terms.len() == 1
            && g.terms
                .keys()
                .next()
                .is_some_and(|e| e.iter().all(|&x| x == 0))
            && g.terms.values().next().is_some_and(|c| *c != 0)
    })
}

fn ideals_equal(a: &GroebnerBasis, b: &GroebnerBasis) -> bool {
    for g in a.generators() {
        if !b.contains(g) {
            return false;
        }
    }
    for g in b.generators() {
        if !a.contains(g) {
            return false;
        }
    }
    true
}

/// Embed `p ∈ R` into `R[t]` with `t` the first exponent (`t^0 · …`).
fn embed_add_t_front(p: &GbPoly) -> GbPoly {
    let n = p.n_vars + 1;
    let mut terms = BTreeMap::new();
    for (e, c) in &p.terms {
        let mut ne = Vec::with_capacity(n);
        ne.push(0u32);
        ne.extend_from_slice(e);
        terms.insert(ne, c.clone());
    }
    GbPoly { terms, n_vars: n }
}

/// `1 - t·f` in `R[t,…]` (first variable is `t`).
fn one_minus_t_times_f(f: &GbPoly) -> GbPoly {
    let fe = embed_add_t_front(f);
    let mut terms: BTreeMap<Vec<u32>, rug::Rational> = BTreeMap::new();
    for (e, c) in &fe.terms {
        let mut ne = e.clone();
        ne[0] += 1;
        let entry = terms.entry(ne).or_insert_with(|| rug::Rational::from(0));
        *entry -= c;
    }
    let zero = vec![0u32; fe.n_vars];
    let one = terms.entry(zero).or_insert_with(|| rug::Rational::from(0));
    *one += 1;
    terms.retain(|_, v| *v != 0);
    GbPoly {
        terms,
        n_vars: fe.n_vars,
    }
}

/// Saturation `I : ⟨f⟩^∞` as `(I, 1 - t·f) ∩ R` (new variable `t` first in Lex).
fn saturate_ideal(generators: &[GbPoly], f: &GbPoly, order: MonomialOrder) -> GroebnerBasis {
    let mut ext: Vec<GbPoly> = Vec::with_capacity(generators.len() + 1);
    for g in generators {
        ext.push(embed_add_t_front(g));
    }
    ext.push(one_minus_t_times_f(f));
    let gb_ext = GroebnerBasis::compute(ext, order);
    let elim = gb_ext.eliminate(&[0]);
    let stripped: Vec<GbPoly> = elim.generators().iter().map(strip_first_var).collect();
    GroebnerBasis::compute(stripped, order)
}

fn strip_first_var(p: &GbPoly) -> GbPoly {
    let old_n = p.n_vars;
    if old_n == 0 {
        return GbPoly::zero(0);
    }
    let n = old_n - 1;
    let mut terms = BTreeMap::new();
    for (e, c) in &p.terms {
        if e.len() != old_n || e[0] != 0 {
            continue;
        }
        if n == 0 {
            terms.insert(vec![], c.clone());
        } else {
            terms.insert(e[1..].to_vec(), c.clone());
        }
    }
    GbPoly { terms, n_vars: n }
}

/// `I ∩ J` from `t·I + (t-1)·J` (first variable is `t`).
fn ideal_intersection(i: &[GbPoly], j: &[GbPoly], order: MonomialOrder) -> GroebnerBasis {
    let mut ext = Vec::with_capacity(i.len() + j.len());
    for g in i {
        let ge = embed_add_t_front(g);
        ext.push(mul_t(&ge));
    }
    for h in j {
        let he = embed_add_t_front(h);
        let th = mul_t(&he);
        ext.push(th.sub(&he));
    }
    let gb_ext = GroebnerBasis::compute(ext, order);
    let elim = gb_ext.eliminate(&[0]);
    let stripped: Vec<GbPoly> = elim.generators().iter().map(strip_first_var).collect();
    GroebnerBasis::compute(stripped, order)
}

fn mul_t(p: &GbPoly) -> GbPoly {
    let mut terms = BTreeMap::new();
    for (e, c) in &p.terms {
        let mut ne = e.clone();
        if ne.is_empty() {
            continue;
        }
        ne[0] += 1;
        terms.insert(ne, c.clone());
    }
    GbPoly {
        terms,
        n_vars: p.n_vars,
    }
}

fn var_monomial(n_vars: usize, idx: usize) -> GbPoly {
    let mut exp = vec![0u32; n_vars];
    exp[idx] = 1;
    GbPoly::monomial(exp, rug::Rational::from(1))
}

fn decompose_recursive(
    gb: GroebnerBasis,
    order: MonomialOrder,
    depth: usize,
) -> Result<Vec<PrimaryComponent>, PrimaryDecompositionError> {
    if depth > MAX_SPLIT_DEPTH {
        return Err(PrimaryDecompositionError::RecursionDepth);
    }
    if is_unit_ideal(&gb) {
        return Ok(vec![]);
    }
    if gb.generators().iter().all(|g| g.is_zero()) {
        // ⟨0⟩ is prime in the domain ℚ[x₁,…,xₙ], hence primary, and is its own
        // radical.
        return Ok(vec![PrimaryComponent {
            primary: gb.clone(),
            associated_prime: gb,
        }]);
    }

    let n_vars = gb.generators()[0].n_vars;

    // Monomial ideals decompose completely and every piece is certified
    // irreducible, so there is nothing left for the general machinery to do.
    if let Some(mons) = monomial_ideal_generators(gb.generators()) {
        return match decompose_monomial_ideal(&mons, n_vars, order) {
            Some(comps) => Ok(comps),
            None => Err(refuse(RefusalSite::Decomposition, MONOMIAL_TOO_LARGE)),
        };
    }

    // Principal ideals decompose completely in the UFD ℚ[x₁,…,xₙ].
    if gb.generators().len() == 1 {
        return match principal_components(&gb.generators()[0], order) {
            Some(comps) => Ok(comps),
            None => Err(refuse(
                RefusalSite::Decomposition,
                MULTIVARIATE_FACTORIZATION_FAILED,
            )),
        };
    }

    for i in 0..n_vars {
        let f = var_monomial(n_vars, i);
        let sat_gb = saturate_ideal(gb.generators(), &f, order);
        let mut sum_gens = gb.generators().to_vec();
        sum_gens.push(var_monomial(n_vars, i));
        let sum_gb = GroebnerBasis::compute(sum_gens, order);

        if is_unit_ideal(&sat_gb) || is_unit_ideal(&sum_gb) {
            continue;
        }
        if ideals_equal(&sat_gb, &gb) || ideals_equal(&sum_gb, &gb) {
            continue;
        }
        let inter = ideal_intersection(sat_gb.generators(), sum_gb.generators(), order);
        if !ideals_equal(&inter, &gb) {
            continue;
        }

        let left = decompose_recursive(sat_gb, order, depth + 1)?;
        let right = decompose_recursive(sum_gb, order, depth + 1)?;
        let mut out = left;
        out.extend(right);
        return Ok(out);
    }

    if let Some(pieces) = try_univariate_factor_split(gb.generators(), n_vars)? {
        let mut acc = Vec::new();
        for piece_gens in pieces {
            let piece = GroebnerBasis::compute(piece_gens, order);
            acc.extend(decompose_recursive(piece, order, depth + 1)?);
        }
        return Ok(acc);
    }

    if let Some(component) = certify_shape_position(gb.generators(), n_vars, order)? {
        return Ok(vec![component]);
    }

    // An ideal whose radical is *maximal* is primary: `R/I` is then local with
    // its maximal ideal equal to the nilradical `√I/I`, so everything outside
    // that ideal is a unit and everything inside is nilpotent — no zero divisor
    // is left over.  This is what certifies `⟨x² + y², xy⟩`, whose radical is
    // `⟨x, y⟩`.
    if let Some(rad) = radical_direct(&gb, order) {
        if let Some(certified) = certify_shape_position(rad.generators(), n_vars, order)? {
            if ideals_equal(&certified.primary, &certified.associated_prime) {
                // A shape-position basis whose eliminant is irreducible *and*
                // square-free has the field ℚ[t]/⟨p⟩ as its quotient, i.e. it is
                // maximal.
                return Ok(vec![PrimaryComponent {
                    primary: gb,
                    associated_prime: rad,
                }]);
            }
        }
    }

    Err(refuse(
        RefusalSite::Decomposition,
        DECOMPOSITION_NOT_CERTIFIABLE,
    ))
}

/// If the basis contains a univariate polynomial in *some* variable whose
/// factorization over ℚ has several distinct irreducible factors, split along
/// them as `⟨I, p_j^{e_j}⟩`.
///
/// The pieces are pairwise comaximal, because `gcd(p_i^{e_i}, p_j^{e_j}) = 1` in
/// the univariate ring ℚ\[x\] gives a Bézout identity, and their product lies in
/// `I`, so CRT makes `I = ∩_j (I + ⟨p_j^{e_j}⟩)` an equality rather than a
/// containment.
///
/// Until 3.8 this looked only in variable 0, so a zero-dimensional ideal in
/// shape position — whose eliminant lives in the *last* variable — was never
/// split at all.
fn try_univariate_factor_split(
    gens: &[GbPoly],
    n_vars: usize,
) -> Result<Option<Vec<Vec<GbPoly>>>, PrimaryDecompositionError> {
    for var in 0..n_vars {
        let u = match find_any_univariate(gens, var) {
            Some(u) => u,
            None => continue,
        };
        let facs = factor_univariate_q_monic(&u, var, n_vars)?;
        if facs.len() <= 1 {
            continue;
        }
        let mut out = Vec::with_capacity(facs.len());
        for (p, e) in facs {
            let mut g = gens.to_vec();
            g.push(gbpoly_pow(&p, e));
            out.push(g);
        }
        return Ok(Some(out));
    }
    Ok(None)
}

fn gbpoly_pow(p: &GbPoly, e: u32) -> GbPoly {
    let mut acc = GbPoly::constant(rug::Rational::from(1), p.n_vars);
    for _ in 0..e {
        acc = acc.mul(p);
    }
    acc
}

fn flint_monic_to_gbpoly_in_var(fz: &FlintPoly, var: usize, n_vars: usize) -> GbPoly {
    let deg = fz.degree();
    if deg < 0 {
        return GbPoly::zero(n_vars);
    }
    let lc = fz.get_coeff_flint(deg as usize).to_rug();
    let mut terms = BTreeMap::new();
    for d in 0..=deg as usize {
        let cz = fz.get_coeff_flint(d).to_rug();
        if cz == 0 {
            continue;
        }
        let rq = rug::Rational::from((cz.clone(), lc.clone()));
        let mut expv = vec![0u32; n_vars];
        expv[var] = d as u32;
        terms.insert(expv, rq);
    }
    GbPoly { terms, n_vars }
}

/// Factor a univariate `p(x_var)`; returns **monic** irreducible factors over ℚ.
fn factor_univariate_q_monic(
    p: &GbPoly,
    var: usize,
    n_vars: usize,
) -> Result<Vec<(GbPoly, u32)>, PrimaryDecompositionError> {
    if var >= n_vars || !is_univariate_in_var(p, var) {
        return Err(PrimaryDecompositionError::Factorization(
            "internal: expected a univariate polynomial in the requested variable",
        ));
    }
    let mut coeff_map: BTreeMap<u32, rug::Rational> = BTreeMap::new();
    for (e, c) in &p.terms {
        coeff_map.insert(e[var], c.clone());
    }
    let deg = *coeff_map.keys().max().unwrap_or(&0);
    let coeffs_r: Vec<rug::Rational> = (0..=deg)
        .map(|d| {
            coeff_map
                .get(&d)
                .cloned()
                .unwrap_or_else(|| rug::Rational::from(0))
        })
        .collect();
    let fp = primitive_flint_from_rational_asc(&coeffs_r).ok_or(
        PrimaryDecompositionError::Factorization(
            "univariate factorization failed: could not build an integer model of the \
             eliminant",
        ),
    )?;
    let (_unit, facs) = fp
        .factor_over_z()
        .map_err(|_| PrimaryDecompositionError::Factorization(FACTORIZATION_FAILED))?;
    let mut pairs = Vec::new();
    for (fz, exp) in facs {
        if fz.degree() < 1 {
            continue;
        }
        let g = flint_monic_to_gbpoly_in_var(&fz, var, n_vars);
        pairs.push((g, exp));
    }
    Ok(pairs)
}

fn is_univariate_in_var(p: &GbPoly, var: usize) -> bool {
    p.terms
        .keys()
        .all(|e| e.len() == p.n_vars && e.iter().enumerate().all(|(i, &v)| i == var || v == 0))
}

// ---------------------------------------------------------------------------
// Radical
// ---------------------------------------------------------------------------

/// `√I` when it can be certified, `None` when it cannot.
///
/// The order of the attempts is cheapest-first; every one of them is exact.
/// Before 3.8 this function's fallback was `I` itself, which is what made
/// `radical` and `PrimaryComponent::associated_prime` able to be wrong.
fn radical_from_basis(gb: &GroebnerBasis, order: MonomialOrder) -> Option<GroebnerBasis> {
    if let Some(r) = radical_direct(gb, order) {
        return Some(r);
    }

    // Last resort: √I is the intersection of the associated primes of any
    // primary decomposition (the embedded ones contain a minimal one, so they
    // do not change the intersection).  This only succeeds where the
    // decomposition itself is certified.
    let comps = decompose_recursive(gb.clone(), order, 0).ok()?;
    if comps.is_empty() {
        return Some(gb.clone());
    }
    let mut acc = comps[0].associated_prime.clone();
    for c in &comps[1..] {
        acc = ideal_intersection(acc.generators(), c.associated_prime.generators(), order);
    }
    Some(acc)
}

/// The part of [`radical_from_basis`] that does not go through
/// [`decompose_recursive`] — the decomposition consults it, so it must not
/// consult the decomposition.
fn radical_direct(gb: &GroebnerBasis, order: MonomialOrder) -> Option<GroebnerBasis> {
    let gens = gb.generators();
    if gens.is_empty() || gens.iter().all(|g| g.is_zero()) {
        // ⟨0⟩ is prime in a domain, so it is its own radical.
        return Some(gb.clone());
    }
    if is_unit_ideal(gb) {
        return Some(gb.clone());
    }
    let n = gens[0].n_vars;

    // √⟨monomials⟩ is generated by the square-free part of each monomial: that
    // ideal J is radical (a monomial ideal generated by square-free monomials
    // is), and each of its generators has a power in I, so I ⊆ J ⊆ √I ⊆ √J = J.
    if let Some(mons) = monomial_ideal_generators(gens) {
        let sq: Vec<GbPoly> = mons
            .iter()
            .map(|m| monomial_poly(&squarefree_exp(m)))
            .collect();
        return Some(GroebnerBasis::compute(sq, order));
    }

    // √⟨f⟩ = ⟨∏ pᵢ⟩ over the distinct irreducible factors of f: ℚ[x₁,…,xₙ] is a
    // UFD, so gⁿ ∈ ⟨f⟩ for some n iff every pᵢ divides g.
    if gens.len() == 1 {
        let facs = factor_gbpoly_q(&gens[0], order)?;
        let mut prod = GbPoly::constant(rug::Rational::from(1), n);
        for (p, _) in &facs {
            prod = prod.mul(p);
        }
        return Some(GroebnerBasis::compute(vec![prod], order));
    }

    radical_zero_dimensional(gb, n, order)
}

/// Seidenberg: over a perfect field, an ideal containing a square-free
/// univariate polynomial in *every* variable is radical.
///
/// So when `I` is zero-dimensional — equivalently, `I ∩ ℚ[xᵢ] ≠ 0` for every
/// `i` — the ideal `J = I + ⟨sqfree(pᵢ) : i⟩` is radical, and `I ⊆ J ⊆ √I`
/// because `pᵢ | sqfree(pᵢ)^{deg pᵢ}`; hence `√I = √J = J`.
///
/// The old code only used the eliminants that happened to appear *verbatim* in
/// the basis. Under Lex that is normally just the last variable, so the result
/// was an ideal between `I` and `√I` reported as `√I`.
fn radical_zero_dimensional(
    gb: &GroebnerBasis,
    n: usize,
    order: MonomialOrder,
) -> Option<GroebnerBasis> {
    if n == 0 {
        return None;
    }
    // `is_zero_dimensional` reads GRevLex leading monomials, so it needs a
    // GRevLex basis — cheaper than discovering the answer through n failed
    // eliminations.
    let grevlex = GroebnerBasis::compute(gb.generators().to_vec(), MonomialOrder::GRevLex);
    if !is_zero_dimensional(grevlex.generators(), n) {
        return None;
    }
    let mut gens = gb.generators().to_vec();
    for i in 0..n {
        let u = eliminant_in_var(gb.generators(), n, i)?;
        gens.push(univariate_squarefree_part(&u, i, order));
    }
    Some(GroebnerBasis::compute(gens, order))
}

/// A non-zero element of `I ∩ ℚ[x_i]`, or `None` when there is none.
fn eliminant_in_var(gens: &[GbPoly], n: usize, i: usize) -> Option<GbPoly> {
    if let Some(u) = find_any_univariate(gens, i) {
        return Some(u);
    }
    // Re-rank the variables so that x_i is lex-least, then the elimination
    // theorem says the basis elements free of every other variable generate
    // `I ∩ ℚ[x_i]`.
    let mut perm: Vec<usize> = (0..n).filter(|&j| j != i).collect();
    perm.push(i);
    let permuted: Vec<GbPoly> = gens.iter().map(|g| permute_vars(g, &perm)).collect();
    let pgb = GroebnerBasis::compute(permuted, MonomialOrder::Lex);
    let u = find_any_univariate(pgb.generators(), n - 1)?;
    let mut terms = BTreeMap::new();
    for (e, c) in &u.terms {
        let mut ne = vec![0u32; n];
        ne[i] = e[n - 1];
        terms.insert(ne, c.clone());
    }
    Some(GbPoly { terms, n_vars: n })
}

/// Re-index the exponent vectors: `new[k] = old[perm[k]]`.
fn permute_vars(p: &GbPoly, perm: &[usize]) -> GbPoly {
    let n = perm.len();
    let mut terms: BTreeMap<Vec<u32>, rug::Rational> = BTreeMap::new();
    for (e, c) in &p.terms {
        let ne: Vec<u32> = perm
            .iter()
            .map(|&j| e.get(j).copied().unwrap_or(0))
            .collect();
        terms.insert(ne, c.clone());
    }
    GbPoly { terms, n_vars: n }
}

// ---------------------------------------------------------------------------
// Monomial ideals
// ---------------------------------------------------------------------------

/// The exponent vectors of `gens` when every generator is a single monomial.
fn monomial_ideal_generators(gens: &[GbPoly]) -> Option<Vec<Vec<u32>>> {
    let mut out = Vec::with_capacity(gens.len());
    for g in gens {
        if g.is_zero() {
            continue;
        }
        if g.terms.len() != 1 {
            return None;
        }
        let (e, c) = g.terms.iter().next()?;
        if *c == 0 {
            return None;
        }
        out.push(e.clone());
    }
    if out.is_empty() {
        None
    } else {
        Some(out)
    }
}

fn monomial_poly(exp: &[u32]) -> GbPoly {
    GbPoly::monomial(exp.to_vec(), rug::Rational::from(1))
}

fn squarefree_exp(m: &[u32]) -> Vec<u32> {
    m.iter().map(|&e| e.min(1)).collect()
}

fn monomial_divides(a: &[u32], b: &[u32]) -> bool {
    a.iter().zip(b.iter()).all(|(x, y)| x <= y)
}

/// The minimal generating set of a monomial ideal: drop every generator that a
/// different one divides.
fn minimal_monomial_generators(mut mons: Vec<Vec<u32>>) -> Vec<Vec<u32>> {
    mons.sort();
    mons.dedup();
    mons.iter()
        .filter(|m| !mons.iter().any(|o| o != *m && monomial_divides(o, m)))
        .cloned()
        .collect()
}

/// Irreducible decomposition of a monomial ideal, by the coprime split.
///
/// For monomial ideals `⟨J, u·v⟩ = ⟨J, u⟩ ∩ ⟨J, v⟩` whenever `u` and `v` are
/// coprime monomials: a monomial in both right-hand ideals that is not in `J`
/// is divisible by `u` and by `v`, hence by `uv`. Splitting the first variable
/// off each generator that mentions more than one strictly decreases
/// `Σ (|supp m| − 1)`, so the recursion terminates, and what it terminates on is
/// a set of pure prime powers `⟨x_{j₁}^{a₁}, …⟩` — an irreducible monomial
/// ideal, primary with associated prime `⟨x_{j₁}, …⟩`.
fn collect_irreducible_monomial_components(
    mons: Vec<Vec<u32>>,
    n: usize,
    out: &mut Vec<Vec<Vec<u32>>>,
) -> bool {
    if out.len() >= MAX_MONOMIAL_COMPONENTS {
        return false;
    }
    let mons = minimal_monomial_generators(mons);
    let split_at = mons
        .iter()
        .position(|m| m.iter().filter(|&&e| e > 0).count() >= 2);
    match split_at {
        None => {
            out.push(mons);
            true
        }
        Some(i) => {
            let m = mons[i].clone();
            let j = m
                .iter()
                .position(|&e| e > 0)
                .expect("a generator with ≥2 variables has support");
            let mut u = vec![0u32; n];
            u[j] = m[j];
            let mut v = m;
            v[j] = 0;
            let mut left = mons.clone();
            left[i] = u;
            let mut right = mons;
            right[i] = v;
            collect_irreducible_monomial_components(left, n, out)
                && collect_irreducible_monomial_components(right, n, out)
        }
    }
}

fn decompose_monomial_ideal(
    mons: &[Vec<u32>],
    n: usize,
    order: MonomialOrder,
) -> Option<Vec<PrimaryComponent>> {
    let mut pieces: Vec<Vec<Vec<u32>>> = Vec::new();
    if !collect_irreducible_monomial_components(mons.to_vec(), n, &mut pieces) {
        return None;
    }
    Some(
        pieces
            .into_iter()
            .map(|pure| {
                let primary =
                    GroebnerBasis::compute(pure.iter().map(|m| monomial_poly(m)).collect(), order);
                let associated_prime = GroebnerBasis::compute(
                    pure.iter()
                        .map(|m| monomial_poly(&squarefree_exp(m)))
                        .collect(),
                    order,
                );
                PrimaryComponent {
                    primary,
                    associated_prime,
                }
            })
            .collect(),
    )
}

// ---------------------------------------------------------------------------
// Principal ideals
// ---------------------------------------------------------------------------

/// `⟨∏ pᵢ^{eᵢ}⟩ = ∩ ⟨pᵢ^{eᵢ}⟩`, each factor primary.
///
/// ℚ\[x₁,…,xₙ\] is a UFD, so the `pᵢ^{eᵢ}` are pairwise coprime and their
/// intersection is generated by their product; and `⟨pᵉ⟩` is `⟨p⟩`-primary for
/// a prime element `p`, because `ab ∈ ⟨pᵉ⟩` with `p ∤ a` forces `pᵉ | b`.
fn principal_components(f: &GbPoly, order: MonomialOrder) -> Option<Vec<PrimaryComponent>> {
    let facs = factor_gbpoly_q(f, order)?;
    if facs.is_empty() {
        return None;
    }
    Some(
        facs.into_iter()
            .map(|(p, e)| PrimaryComponent {
                primary: GroebnerBasis::compute(vec![gbpoly_pow(&p, e)], order),
                associated_prime: GroebnerBasis::compute(vec![p], order),
            })
            .collect(),
    )
}

/// Non-constant irreducible factors of `p` over ℚ, made monic, with
/// multiplicity. `None` when FLINT declines.
fn factor_gbpoly_q(p: &GbPoly, order: MonomialOrder) -> Option<Vec<(GbPoly, u32)>> {
    let n = p.n_vars;
    if n == 0 || p.is_zero() {
        return None;
    }
    let ctx = FlintMPolyCtx::new(n);
    let fp = gbpoly_to_flint(p, &ctx)?;
    let mut fac = FlintMPolyFactor::new(Arc::clone(&ctx));
    if !fac.factor(&fp) || !fac.constant_den_is_one() {
        return None;
    }
    let mut out = Vec::with_capacity(fac.len());
    for i in 0..fac.len() {
        let base = flint_to_gbpoly(&fac.base_at(i), n);
        let exp = fac.exp_at(i);
        if base.is_zero() || is_constant(&base) {
            continue;
        }
        out.push((base.make_monic(order), exp));
    }
    Some(out)
}

fn is_constant(p: &GbPoly) -> bool {
    p.terms.keys().all(|e| e.iter().all(|&v| v == 0))
}

/// Clear denominators and hand the integer model to FLINT.
fn gbpoly_to_flint(p: &GbPoly, ctx: &Arc<FlintMPolyCtx>) -> Option<FlintMPoly> {
    if p.is_zero() {
        return None;
    }
    let coeffs: Vec<rug::Rational> = p.terms.values().cloned().collect();
    let lcm = lcm_rational_denoms(&coeffs);
    let nv = ctx.nvars();
    let mut fp = FlintMPoly::new(Arc::clone(ctx));
    for (e, c) in &p.terms {
        let scaled = c.clone() * rug::Rational::from((lcm.clone(), 1));
        let (num, den) = scaled.into_numer_denom();
        debug_assert_eq!(den, rug::Integer::from(1));
        let mut exp = vec![0u64; nv];
        for (i, &v) in e.iter().enumerate() {
            if i < nv {
                exp[i] = u64::from(v);
            }
        }
        fp.push_term(&num, &exp);
    }
    fp.finish();
    Some(fp)
}

fn flint_to_gbpoly(f: &FlintMPoly, n_vars: usize) -> GbPoly {
    let mut terms: BTreeMap<Vec<u32>, rug::Rational> = BTreeMap::new();
    for (e, c) in f.terms() {
        if c == 0 {
            continue;
        }
        let mut exp = vec![0u32; n_vars];
        for (i, &v) in e.iter().enumerate() {
            if i < n_vars {
                exp[i] = v;
            }
        }
        terms.insert(exp, rug::Rational::from((c, 1)));
    }
    GbPoly { terms, n_vars }
}

// ---------------------------------------------------------------------------
// Shape-position certificate
// ---------------------------------------------------------------------------

/// `⟨x₀ − g₀(t), …, x_{n−2} − g_{n−2}(t), h(t)⟩` with `t = x_{n−1}`, `h = c·pᵉ`
/// and `p` irreducible over ℚ.
///
/// The quotient is `ℚ[t]/⟨pᵉ⟩` — a local Artinian ring, whose only zero
/// divisors are nilpotent — so the ideal is primary, and replacing `h` by `p`
/// gives the associated prime `ℚ[t]/⟨p⟩`, a field.
///
/// This is the certificate the zero-dimensional leaves of the saturation and
/// CRT splits land on: `⟨x − 1, y⟩` and `⟨x + 1, y⟩` from `⟨x² − 1, y⟩`, for
/// instance.
fn certify_shape_position(
    gens: &[GbPoly],
    n: usize,
    order: MonomialOrder,
) -> Result<Option<PrimaryComponent>, PrimaryDecompositionError> {
    if n == 0 || gens.len() != n {
        return Ok(None);
    }
    let t = n - 1;
    let mut eliminant: Option<GbPoly> = None;
    let mut solved: Vec<(usize, GbPoly)> = Vec::new();
    for g in gens {
        if g.is_zero() {
            return Ok(None);
        }
        if is_univariate_in_var(g, t) {
            if eliminant.is_some() {
                return Ok(None);
            }
            eliminant = Some(g.clone());
        } else if let Some(j) = solved_variable_over_last(g, t) {
            solved.push((j, g.clone()));
        } else {
            return Ok(None);
        }
    }
    let Some(h) = eliminant else {
        return Ok(None);
    };
    let mut seen: Vec<usize> = solved.iter().map(|(j, _)| *j).collect();
    seen.sort_unstable();
    if seen != (0..t).collect::<Vec<_>>() {
        return Ok(None);
    }

    let facs = factor_univariate_q_monic(&h, t, n)?;
    if facs.len() != 1 {
        // Several distinct factors: the CRT split should have taken this ideal
        // apart already, so this is not a leaf.
        return Ok(None);
    }
    let (p, _e) = facs.into_iter().next().expect("exactly one factor");
    let mut prime_gens: Vec<GbPoly> = solved.into_iter().map(|(_, g)| g).collect();
    prime_gens.push(p);
    Ok(Some(PrimaryComponent {
        primary: GroebnerBasis::compute(gens.to_vec(), order),
        associated_prime: GroebnerBasis::compute(prime_gens, order),
    }))
}

/// `Some(j)` when `g = c·x_j + (a polynomial in x_t alone)` with `j ≠ t`.
fn solved_variable_over_last(g: &GbPoly, t: usize) -> Option<usize> {
    let mut lead: Option<usize> = None;
    for e in g.terms.keys() {
        let support: Vec<usize> = e
            .iter()
            .enumerate()
            .filter(|(_, &v)| v > 0)
            .map(|(i, _)| i)
            .collect();
        if support.iter().all(|&i| i == t) {
            continue;
        }
        if support.len() == 1 && support[0] != t && e[support[0]] == 1 {
            if lead.is_some() {
                return None;
            }
            lead = Some(support[0]);
        } else {
            return None;
        }
    }
    lead
}

fn find_any_univariate(gens: &[GbPoly], var: usize) -> Option<GbPoly> {
    for g in gens {
        if g.is_zero() || !is_univariate_in_var(g, var) {
            continue;
        }
        return Some(g.clone());
    }
    None
}

fn univariate_squarefree_part(u: &GbPoly, var: usize, order: MonomialOrder) -> GbPoly {
    let n = u.n_vars;
    let mut coeff_map: BTreeMap<u32, rug::Rational> = BTreeMap::new();
    for (e, c) in &u.terms {
        coeff_map.insert(e[var], c.clone());
    }
    let deg = *coeff_map.keys().max().unwrap_or(&0);
    let coeffs_r: Vec<rug::Rational> = (0..=deg)
        .map(|d| {
            coeff_map
                .get(&d)
                .cloned()
                .unwrap_or_else(|| rug::Rational::from(0))
        })
        .collect();
    let fp = match primitive_flint_from_rational_asc(&coeffs_r) {
        Some(p) => p,
        None => return GbPoly::zero(n),
    };
    let der = fp.derivative();
    let g = fp.gcd(&der);
    let sf = fp.div_exact(&g);
    let mut terms = BTreeMap::new();
    for d in 0..=sf.degree() {
        let cz = sf.get_coeff_flint(d as usize).to_rug();
        if cz == 0 {
            continue;
        }
        let mut expv = vec![0u32; n];
        expv[var] = d as u32;
        terms.insert(expv, rug::Rational::from((cz, 1)));
    }
    GbPoly { terms, n_vars: n }.make_monic(order)
}

/// Drop every component that contains another one.
///
/// `Q_j ⊆ Q_i` makes `Q_i ∩ Q_j = Q_j`, so `Q_i` contributes nothing to the
/// intersection and the decomposition stays exact without it. This is what
/// removes the `⟨x, z⟩` that `⟨xz, yz⟩` used to report alongside `⟨z⟩` — the
/// ideal is square-free monomial, hence radical, so exactly its two minimal
/// primes are associated.
fn drop_redundant_components(v: &mut Vec<PrimaryComponent>) {
    let mut i = 0;
    while i < v.len() {
        let redundant = v
            .iter()
            .enumerate()
            .any(|(j, other)| j != i && ideal_contains_ideal(&v[i].primary, &other.primary));
        if redundant {
            v.remove(i);
        } else {
            i += 1;
        }
    }
}

/// `inner ⊆ outer`.
fn ideal_contains_ideal(outer: &GroebnerBasis, inner: &GroebnerBasis) -> bool {
    inner.generators().iter().all(|g| outer.contains(g))
}

fn dedup_components(v: &mut Vec<PrimaryComponent>) {
    let mut i = 0;
    while i < v.len() {
        let mut dup = false;
        for j in 0..i {
            if ideals_equal(&v[i].primary, &v[j].primary) {
                dup = true;
                break;
            }
        }
        if dup {
            v.remove(i);
        } else {
            i += 1;
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn rat(n: i64, d: i64) -> rug::Rational {
        rug::Rational::from((n, d))
    }

    #[test]
    fn intersection_xy_xz() {
        let n = 3usize;
        let xy = GbPoly {
            terms: [(vec![1, 1, 0], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let xz = GbPoly {
            terms: [(vec![1, 0, 1], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let gb_i = GroebnerBasis::compute(vec![xy, xz], MonomialOrder::Lex);
        let f = var_monomial(n, 0);
        let a = saturate_ideal(gb_i.generators(), &f, MonomialOrder::Lex);
        let mut sg = gb_i.generators().to_vec();
        sg.push(f);
        let b = GroebnerBasis::compute(sg, MonomialOrder::Lex);
        let inter = ideal_intersection(a.generators(), b.generators(), MonomialOrder::Lex);
        assert!(ideals_equal(&inter, &gb_i));
    }

    #[test]
    fn primary_xy_xz() {
        let n = 3usize;
        let xy = GbPoly {
            terms: [(vec![1, 1, 0], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let xz = GbPoly {
            terms: [(vec![1, 0, 1], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let dec = primary_decomposition(vec![xy, xz], MonomialOrder::Lex).unwrap();
        assert_eq!(dec.len(), 2);
    }

    #[test]
    fn primary_x2_xy_embedded() {
        let n = 2usize;
        let x2 = GbPoly {
            terms: [(vec![2, 0], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let xy_ = GbPoly {
            terms: [(vec![1, 1], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let gens = vec![x2.clone(), xy_.clone()];
        let dec = primary_decomposition(gens.clone(), MonomialOrder::Lex).unwrap();
        assert_eq!(dec.len(), 2);
        let r = radical(gens, MonomialOrder::Lex).unwrap();
        let one_x = GbPoly {
            terms: [(vec![1, 0], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        assert!(r.contains(&one_x));
    }

    #[test]
    fn factor_split_x2_minus_one() {
        let n = 2usize;
        let xm1 = GbPoly {
            terms: [(vec![2, 0], rat(1, 1)), (vec![0, 0], rat(-1, 1))]
                .into_iter()
                .collect(),
            n_vars: n,
        };
        let y = GbPoly {
            terms: [(vec![0, 1], rat(1, 1))].into_iter().collect(),
            n_vars: n,
        };
        let dec = primary_decomposition(vec![xm1, y], MonomialOrder::Lex).unwrap();
        assert_eq!(dec.len(), 2);
    }

    // -----------------------------------------------------------------------
    // 3.8 — silent errors in the radical and the associated primes
    // -----------------------------------------------------------------------

    /// `c₀·x^a·y^b + c₁·x^c·y^d + …` from `[(a, b, c₀), …]`.
    fn poly2(terms: &[(u32, u32, i64)]) -> GbPoly {
        GbPoly {
            terms: terms
                .iter()
                .map(|&(a, b, c)| (vec![a, b], rat(c, 1)))
                .collect(),
            n_vars: 2,
        }
    }

    fn poly3(terms: &[(u32, u32, u32, i64)]) -> GbPoly {
        GbPoly {
            terms: terms
                .iter()
                .map(|&(a, b, c, k)| (vec![a, b, c], rat(k, 1)))
                .collect(),
            n_vars: 3,
        }
    }

    #[test]
    fn radical_of_a_square_contains_its_base() {
        // √⟨(x−y)²⟩ = ⟨x−y⟩ exactly: ℚ[x,y]/(x−y) ≅ ℚ[y] is a domain, so ⟨x−y⟩
        // is prime, and it contains (x−y)².  This returned ⟨(x−y)²⟩ before 3.8,
        // so `contains(x−y)` was False while `contains((x−y)²)` was True.
        let sq = poly2(&[(2, 0, 1), (1, 1, -2), (0, 2, 1)]);
        let r = radical(vec![sq.clone()], MonomialOrder::Lex).unwrap();
        let base = poly2(&[(1, 0, 1), (0, 1, -1)]);
        assert!(r.contains(&base), "√⟨(x−y)²⟩ must contain x−y");
        assert!(r.contains(&sq));
        // …and no more than that: y ∉ ⟨x−y⟩.
        assert!(!r.contains(&poly2(&[(0, 1, 1)])));
    }

    #[test]
    fn decomposition_of_a_difference_of_squares_is_two_primes() {
        // x² − y² = (x−y)(x+y), both irreducible and non-associate, so
        // ⟨x²−y²⟩ = ⟨x−y⟩ ∩ ⟨x+y⟩ and neither factor ideal is the whole thing.
        // ⟨x²−y²⟩ itself is *not* primary: (x−y)(x+y) ∈ I, x−y ∉ I by degree,
        // and no power of x+y is divisible by x²−y².
        let f = poly2(&[(2, 0, 1), (0, 2, -1)]);
        let dec = primary_decomposition(vec![f], MonomialOrder::Lex).unwrap();
        assert_eq!(dec.len(), 2);
        let minus = poly2(&[(1, 0, 1), (0, 1, -1)]);
        let plus = poly2(&[(1, 0, 1), (0, 1, 1)]);
        assert!(dec.iter().any(|c| ideals_equal(
            &c.primary,
            &GroebnerBasis::compute(vec![minus.clone()], MonomialOrder::Lex)
        )));
        assert!(dec.iter().any(|c| ideals_equal(
            &c.primary,
            &GroebnerBasis::compute(vec![plus.clone()], MonomialOrder::Lex)
        )));
        for c in &dec {
            assert!(
                ideals_equal(&c.primary, &c.associated_prime),
                "both are prime"
            );
        }
    }

    #[test]
    fn squarefree_monomial_ideal_has_exactly_its_minimal_primes() {
        // ⟨xz, yz⟩ is square-free monomial, hence radical, so its associated
        // primes are exactly its minimal primes: ⟨z⟩ and ⟨x, y⟩.  The third
        // component this used to report, ⟨x, z⟩, contains ⟨z⟩ and is therefore
        // redundant.
        let xz = poly3(&[(1, 0, 1, 1)]);
        let yz = poly3(&[(0, 1, 1, 1)]);
        let dec = primary_decomposition(vec![xz, yz], MonomialOrder::Lex).unwrap();
        assert_eq!(dec.len(), 2);
        let z_only = GroebnerBasis::compute(vec![poly3(&[(0, 0, 1, 1)])], MonomialOrder::Lex);
        let x_and_y = GroebnerBasis::compute(
            vec![poly3(&[(1, 0, 0, 1)]), poly3(&[(0, 1, 0, 1)])],
            MonomialOrder::Lex,
        );
        assert!(dec.iter().any(|c| ideals_equal(&c.primary, &z_only)));
        assert!(dec.iter().any(|c| ideals_equal(&c.primary, &x_and_y)));
        for c in &dec {
            assert!(
                ideals_equal(&c.primary, &c.associated_prime),
                "radical ideal"
            );
        }
    }

    #[test]
    fn embedded_monomial_component_keeps_its_multiplicity() {
        // ⟨x², xy⟩ = ⟨x⟩ ∩ ⟨x², y⟩.  The second is ⟨x,y⟩-primary and *not*
        // prime, so a decomposition that reported prime = primary here would be
        // wrong in the other direction.
        let dec = primary_decomposition(
            vec![poly2(&[(2, 0, 1)]), poly2(&[(1, 1, 1)])],
            MonomialOrder::Lex,
        )
        .unwrap();
        assert_eq!(dec.len(), 2);
        let embedded = dec
            .iter()
            .find(|c| !ideals_equal(&c.primary, &c.associated_prime))
            .expect("⟨x², y⟩ is primary but not prime");
        assert!(embedded.associated_prime.contains(&poly2(&[(1, 0, 1)])));
        assert!(embedded.associated_prime.contains(&poly2(&[(0, 1, 1)])));
        assert!(
            !embedded.primary.contains(&poly2(&[(1, 0, 1)])),
            "x ∉ ⟨x², y⟩"
        );
    }

    #[test]
    fn primary_when_the_radical_is_maximal() {
        // √⟨x² + y², xy⟩ = ⟨x, y⟩: y(x²+y²) − x(xy) = y³ and x(x²+y²) − y(xy) = x³
        // are both in the ideal.  A maximal radical makes the ideal primary.
        let gens = vec![poly2(&[(2, 0, 1), (0, 2, 1)]), poly2(&[(1, 1, 1)])];
        let r = radical(gens.clone(), MonomialOrder::Lex).unwrap();
        assert!(r.contains(&poly2(&[(1, 0, 1)])));
        assert!(r.contains(&poly2(&[(0, 1, 1)])));
        let dec = primary_decomposition(gens, MonomialOrder::Lex).unwrap();
        assert_eq!(dec.len(), 1);
        assert!(dec[0].associated_prime.contains(&poly2(&[(1, 0, 1)])));
        assert!(!dec[0].primary.contains(&poly2(&[(1, 0, 1)])));
    }

    #[test]
    fn radical_refuses_rather_than_return_its_input() {
        // The twisted cubic ⟨y − x², z − x³⟩ is prime and one-dimensional; it is
        // neither monomial nor principal nor zero-dimensional, and none of the
        // splits reduce it.  √I = I here, but nothing in this module can *show*
        // that, so it must refuse instead of guessing right by accident.
        let gens = vec![
            poly3(&[(0, 1, 0, 1), (2, 0, 0, -1)]),
            poly3(&[(0, 0, 1, 1), (3, 0, 0, -1)]),
        ];
        let err = radical(gens.clone(), MonomialOrder::Lex)
            .expect_err("must refuse rather than assert √I = I");
        assert!(matches!(err, PrimaryDecompositionError::Factorization(_)));
        let refusal = take_ideal_refusal().expect("refusal recorded out of band");
        assert_eq!(refusal.code(), "E-IDEAL-005");
        assert_eq!(take_ideal_refusal(), None, "consuming");

        let err = primary_decomposition(gens, MonomialOrder::Lex).expect_err("must refuse");
        assert!(matches!(err, PrimaryDecompositionError::Factorization(_)));
        assert_eq!(
            take_ideal_refusal().expect("refusal recorded").code(),
            "E-IDEAL-006"
        );
    }

    #[test]
    fn every_reported_associated_prime_is_radical() {
        // The weakest check that would have caught the old behaviour on every
        // ideal at once: √P = P for a prime P, so re-radicalising a reported
        // associated prime must be a no-op.
        let cases: Vec<Vec<GbPoly>> = vec![
            vec![poly2(&[(2, 0, 1), (0, 2, -1)])],
            vec![poly2(&[(2, 0, 1)]), poly2(&[(1, 1, 1)])],
            vec![poly3(&[(1, 0, 1, 1)]), poly3(&[(0, 1, 1, 1)])],
            vec![poly2(&[(2, 0, 1), (0, 0, -1)]), poly2(&[(0, 1, 1)])],
            vec![poly2(&[(2, 0, 1), (0, 2, 1)]), poly2(&[(1, 1, 1)])],
        ];
        for gens in cases {
            let dec = primary_decomposition(gens.clone(), MonomialOrder::Lex).unwrap();
            assert!(!dec.is_empty());
            for c in &dec {
                let again = radical(c.associated_prime.generators().to_vec(), MonomialOrder::Lex)
                    .expect("the radical of a certified prime is computable");
                assert!(
                    ideals_equal(&again, &c.associated_prime),
                    "√P ≠ P — the reported associated prime is not prime"
                );
                for g in c.primary.generators() {
                    assert!(
                        c.associated_prime.contains(g),
                        "Q ⊄ √Q — the component does not lie in its own prime"
                    );
                }
            }
        }
    }
}