alkahest-cas 3.12.0

//! Polynomial reduction (division algorithm) for Gröbner bases.
//!
//! The division algorithm runs on a private working representation (see
//! `Divisors`) rather than on [`GbPoly`] values:
//!
//! * the dividend lives in a `BTreeMap` keyed by an *order-preserving encoding*
//!   of each exponent vector, so its leading term is the map's last entry —
//!   `O(log n)` instead of a full scan per step — and each reduction step
//!   updates it in place instead of rebuilding it through `sub(neg(..))`;
//! * each divisor's leading monomial, leading coefficient, total degree and
//!   tail are extracted once, not rescanned on every step; a basis held by the
//!   Buchberger loop is prepared incrementally, element by element;
//! * `interreduce` skips the element being reduced by index instead of cloning
//!   the rest of the basis once per element;
//! * where only the remainder's *direction* matters — Buchberger and
//!   `interreduce` make every remainder monic — the `*_scaled` variants divide
//!   fraction-free over ℤ (primitive integer divisors, periodic content
//!   removal), which avoids the gcd canonicalisation `mpq` arithmetic pays on
//!   every operation.
//!
//! The choice of divisor at every step (first divisible leading monomial in
//! cyclic order from the last one used) and the exact rational arithmetic are
//! the same as in the straightforward algorithm, so the remainders — and hence
//! the Gröbner bases built from them — are identical term for term.  (Scaling a
//! dividend by a non-zero constant leaves its support unchanged, so the
//! fraction-free runs make the same choices and end at a non-zero multiple of
//! the same remainder, which `make_monic` maps to the same polynomial.)  Inputs the
//! fast path does not model (exponent vectors of mixed length, explicitly
//! stored zero coefficients) are routed to that straightforward algorithm,
//! `reduce_reference`, unchanged.

use std::collections::BTreeMap;

use rug::ops::NegAssign;
use rug::{Assign, Integer, Rational};

use crate::poly::groebner::ideal::GbPoly;
use crate::poly::groebner::monomial_order::MonomialOrder;
// Exponent sums panic rather than wrap: a wrapped exponent is a different
// monomial, i.e. a wrong basis (see `poly::exponent`).
use crate::poly::exponent::{add_or_panic, total_degree_or_panic};

/// Divide `f` by the ordered list `gs` and return the remainder.
///
/// This is the multivariate division algorithm. The remainder is zero if and
/// only if `f` is in the ideal generated by `gs` (when `gs` is a Gröbner
/// basis under the given order).
pub fn reduce(f: &GbPoly, gs: &[GbPoly], order: MonomialOrder) -> GbPoly {
    Divisors::new(gs, order).reduce(f, None)
}

/// Compute the S-polynomial of two polynomials f and g under the given order.
pub fn s_polynomial(f: &GbPoly, g: &GbPoly, order: MonomialOrder) -> GbPoly {
    let lf_exp = match f.leading_exp(order) {
        Some(e) => e,
        None => return GbPoly::zero(f.n_vars),
    };
    let lg_exp = match g.leading_exp(order) {
        Some(e) => e,
        None => return GbPoly::zero(f.n_vars),
    };
    let lf_coeff = f.leading_coeff(order).unwrap();
    let lg_coeff = g.leading_coeff(order).unwrap();

    // lcm(lm(f), lm(g))
    let lcm_exp: Vec<u32> = lf_exp
        .iter()
        .zip(lg_exp.iter())
        .map(|(&a, &b)| a.max(b))
        .collect();

    let shift_f: Vec<u32> = lcm_exp
        .iter()
        .zip(lf_exp.iter())
        .map(|(l, a)| l - a)
        .collect();
    let shift_g: Vec<u32> = lcm_exp
        .iter()
        .zip(lg_exp.iter())
        .map(|(l, b)| l - b)
        .collect();

    // S(f, g) = (1/lc(f)) * x^shift_f * f - (1/lc(g)) * x^shift_g * g
    let coeff_f = rug::Rational::from(1) / &lf_coeff;
    let coeff_g = rug::Rational::from(1) / &lg_coeff;

    let term_f = f.mul_monomial(&shift_f, &coeff_f);
    let term_g = g.mul_monomial(&shift_g, &coeff_g);

    term_f.sub(&term_g)
}

// ---------------------------------------------------------------------------
// In-place reduction engine
// ---------------------------------------------------------------------------

/// Encode `e + shift` so that the natural (lexicographic) order on the encoded
/// vectors coincides with `order` on exponent vectors of one fixed length.
///
/// * Lex: the exponents themselves.
/// * GrLex: `[deg, e_0, …, e_{n-1}]`.
/// * GRevLex: `[deg, MAX - e_{n-1}, …, MAX - e_0]` — among equal degrees the
///   *smaller* exponent of the *last* variable ranks higher.
///
/// For the graded orders `key[0]` is the total degree.
#[inline]
fn encode_shifted(order: MonomialOrder, e: &[u32], shift: &[u32], out: &mut Vec<u32>) {
    out.clear();
    match order {
        MonomialOrder::Lex => out.extend(e.iter().zip(shift).map(|(&a, &b)| add_or_panic(a, b))),
        MonomialOrder::GrLex => {
            out.push(0);
            let mut deg = 0u32;
            for (&a, &b) in e.iter().zip(shift) {
                let x = add_or_panic(a, b);
                deg = add_or_panic(deg, x);
                out.push(x);
            }
            out[0] = deg;
        }
        MonomialOrder::GRevLex => {
            out.push(0);
            let mut deg = 0u32;
            for (&a, &b) in e.iter().zip(shift).rev() {
                let x = add_or_panic(a, b);
                deg = add_or_panic(deg, x);
                out.push(u32::MAX - x);
            }
            out[0] = deg;
        }
    }
}

/// Inverse of [`encode_shifted`] (with a zero shift).
#[inline]
fn decode(order: MonomialOrder, key: &[u32], out: &mut Vec<u32>) {
    out.clear();
    match order {
        MonomialOrder::Lex => out.extend_from_slice(key),
        MonomialOrder::GrLex => out.extend_from_slice(&key[1..]),
        MonomialOrder::GRevLex => out.extend(key[1..].iter().rev().map(|k| u32::MAX - k)),
    }
}

/// Add `±c·q` (`q = None` meaning 1; minus when `negate`) into the working map
/// at `key`, dropping the entry if it cancels.  `tmp` is scratch space whose
/// allocation is handed to the map when a new term is created.
#[inline]
fn accumulate(
    p: &mut BTreeMap<Vec<u32>, Rational>,
    key: &[u32],
    c: &Rational,
    q: Option<&Rational>,
    negate: bool,
    tmp: &mut Rational,
) {
    match q {
        Some(q) => tmp.assign(c * q),
        None => tmp.assign(c),
    }
    match p.get_mut(key) {
        Some(v) => {
            if negate {
                *v -= &*tmp;
            } else {
                *v += &*tmp;
            }
            if v.cmp0() == std::cmp::Ordering::Equal {
                p.remove(key);
            }
        }
        None => {
            if negate {
                tmp.neg_assign();
            }
            p.insert(key.to_vec(), std::mem::take(tmp));
        }
    }
}

/// A polynomial prepared for use as a divisor: its leading data extracted once.
#[derive(Debug, Clone)]
pub(crate) struct Divisor {
    lead: Vec<u32>,
    lead_deg: u32,
    lc: Rational,
    lc_is_one: bool,
    /// Every term but the leading one, in the polynomial's own exponent form.
    tail: Vec<(Vec<u32>, Rational)>,
    n_vars: usize,
    /// `Some(w)` when every exponent vector has length `w` and no coefficient
    /// is an explicitly stored zero — the inputs the fast path models exactly.
    width: Option<usize>,
    /// Primitive integer multiple, for the fraction-free path (when enabled).
    int: Option<IntForm>,
}

/// A polynomial's primitive integer multiple `c·g` (`c ∈ ℚ*`): leading
/// coefficient and tail coefficients, aligned with [`Divisor::tail`].
#[derive(Debug, Clone)]
struct IntForm {
    lc: Integer,
    tail: Vec<Integer>,
}

impl IntForm {
    fn new(lc: &Rational, tail: &[(Vec<u32>, Rational)]) -> Self {
        let mut ints =
            primitive_integer_multiple(std::iter::once(lc).chain(tail.iter().map(|(_, c)| c)));
        let tail = ints.split_off(1);
        let lc = ints.pop().expect("leading coefficient");
        IntForm { lc, tail }
    }
}

/// The coefficients of the primitive integer multiple of a coefficient list:
/// clear denominators by their lcm, then divide by the gcd of the numerators.
fn primitive_integer_multiple<'a>(
    coeffs: impl Iterator<Item = &'a Rational> + Clone,
) -> Vec<Integer> {
    let mut den = Integer::from(1);
    for c in coeffs.clone() {
        den.lcm_mut(c.denom());
    }
    let mut ints: Vec<Integer> = coeffs
        .map(|c| {
            let mut m = Integer::from(den.div_exact_ref(c.denom()));
            m *= c.numer();
            m
        })
        .collect();
    let mut content = Integer::new();
    for c in &ints {
        content.gcd_mut(c);
        if content == 1 {
            return ints;
        }
    }
    if content != 0 {
        for c in &mut ints {
            c.div_exact_mut(&content);
        }
    }
    ints
}

/// Reductions that rescale the dividend between two content removals in the
/// fraction-free loop.
const CONTENT_EVERY: u32 = 2;

/// Integer analogue of [`accumulate`]: add `±c·q` into the map at `key`.
#[inline]
fn accumulate_int(
    p: &mut BTreeMap<Vec<u32>, Integer>,
    key: &[u32],
    c: &Integer,
    q: &Integer,
    negate: bool,
    tmp: &mut Integer,
) {
    tmp.assign(c * q);
    match p.get_mut(key) {
        Some(v) => {
            if negate {
                *v -= &*tmp;
            } else {
                *v += &*tmp;
            }
            if v.cmp0() == std::cmp::Ordering::Equal {
                p.remove(key);
            }
        }
        None => {
            if negate {
                tmp.neg_assign();
            }
            p.insert(key.to_vec(), std::mem::take(tmp));
        }
    }
}

/// Divide the working dividend and the remainder so far by their common
/// integer content (a positive scalar, so the pair stays proportional).
fn remove_content(
    p: &mut BTreeMap<Vec<u32>, Integer>,
    rem: &mut [(Vec<u32>, Integer)],
    g: &mut Integer,
) {
    g.assign(0);
    for c in p.values().chain(rem.iter().map(|(_, c)| c)) {
        g.gcd_mut(c);
        if *g == 1 {
            return;
        }
    }
    if *g == 0 {
        return;
    }
    for c in p.values_mut().chain(rem.iter_mut().map(|(_, c)| c)) {
        c.div_exact_mut(g);
    }
}

impl Divisor {
    /// `None` for the zero polynomial (which never divides anything).
    fn new(g: &GbPoly, order: MonomialOrder, fraction_free: bool) -> Option<Self> {
        let (lead, lc) = g.leading_term(order)?;
        let w = lead.len();
        let regular = g.terms.iter().all(|(e, c)| e.len() == w && *c != 0);
        let tail: Vec<(Vec<u32>, Rational)> = g
            .terms
            .iter()
            .filter(|(e, _)| *e != lead)
            .map(|(e, c)| (e.clone(), c.clone()))
            .collect();
        let int = (fraction_free && regular).then(|| IntForm::new(lc, &tail));
        Some(Divisor {
            int,
            lead: lead.clone(),
            lead_deg: total_degree_or_panic(lead),
            lc_is_one: *lc == 1,
            lc: lc.clone(),
            tail,
            n_vars: g.n_vars,
            width: regular.then_some(w),
        })
    }

    /// Rebuild the polynomial (for the reference fallback).
    fn to_poly(&self) -> GbPoly {
        let mut terms: BTreeMap<Vec<u32>, Rational> = self.tail.iter().cloned().collect();
        terms.insert(self.lead.clone(), self.lc.clone());
        GbPoly {
            terms,
            n_vars: self.n_vars,
        }
    }
}

/// An ordered list of divisors (a basis, typically) prepared for repeated
/// reduction under one monomial order.  Indices mirror the caller's list; zero
/// polynomials occupy a slot but never divide.
#[derive(Debug, Clone)]
pub(crate) struct Divisors {
    order: MonomialOrder,
    divs: Vec<Option<Divisor>>,
    /// `n_vars` of each slot, zero polynomials included.
    n_vars: Vec<usize>,
    /// Also prepare integer forms, for [`Self::reduce_s_pair_scaled`] and
    /// [`Self::reduce_scaled`].
    fraction_free: bool,
}

impl Divisors {
    pub(crate) fn new(gs: &[GbPoly], order: MonomialOrder) -> Self {
        Divisors {
            order,
            divs: gs.iter().map(|g| Divisor::new(g, order, false)).collect(),
            n_vars: gs.iter().map(|g| g.n_vars).collect(),
            fraction_free: false,
        }
    }

    /// Like [`Self::new`], additionally preparing each divisor's primitive
    /// integer form so the `*_scaled` (fraction-free) reductions can be used.
    pub(crate) fn new_fraction_free(gs: &[GbPoly], order: MonomialOrder) -> Self {
        Divisors {
            order,
            divs: gs.iter().map(|g| Divisor::new(g, order, true)).collect(),
            n_vars: gs.iter().map(|g| g.n_vars).collect(),
            fraction_free: true,
        }
    }

    pub(crate) fn push(&mut self, g: &GbPoly) {
        self.divs
            .push(Divisor::new(g, self.order, self.fraction_free));
        self.n_vars.push(g.n_vars);
    }

    pub(crate) fn set(&mut self, i: usize, g: &GbPoly) {
        self.divs[i] = Divisor::new(g, self.order, self.fraction_free);
        self.n_vars[i] = g.n_vars;
    }

    pub(crate) fn remove(&mut self, i: usize) {
        self.divs.remove(i);
        self.n_vars.remove(i);
    }

    /// The divisors other than `skip`, as polynomials (reference fallback only).
    fn polys_except(&self, skip: Option<usize>) -> Vec<GbPoly> {
        self.divs
            .iter()
            .enumerate()
            .filter(|&(k, _)| Some(k) != skip)
            .map(|(_, d)| d.as_ref().map_or_else(|| GbPoly::zero(0), Divisor::to_poly))
            .collect()
    }

    /// True when every non-zero divisor other than `skip` is regular of width `w`.
    fn all_width(&self, w: usize, skip: Option<usize>) -> bool {
        self.divs
            .iter()
            .enumerate()
            .all(|(k, d)| Some(k) == skip || d.as_ref().map_or(true, |d| d.width == Some(w)))
    }

    /// Remainder of `f` on division by every divisor except `skip`.
    ///
    /// Equal to `reduce(f, gs_without_skip, order)`: removing one element from
    /// a cyclic sequence keeps the cyclic order of the rest, so the
    /// first-divisor-from-the-last-one search visits candidates in the same
    /// order either way.
    pub(crate) fn reduce(&self, f: &GbPoly, skip: Option<usize>) -> GbPoly {
        let order = self.order;
        let Some(first) = f.terms.keys().next() else {
            return GbPoly::zero(f.n_vars);
        };
        let w = first.len();
        let regular = f.terms.iter().all(|(e, c)| e.len() == w && *c != 0);
        if !regular || !self.all_width(w, skip) {
            return reduce_reference(f, &self.polys_except(skip), order);
        }
        let zero_shift = vec![0u32; w];
        let mut key = Vec::with_capacity(w + 1);
        let mut p = BTreeMap::new();
        for (e, c) in &f.terms {
            encode_shifted(order, e, &zero_shift, &mut key);
            p.insert(key.clone(), c.clone());
        }
        self.reduce_work(p, f.n_vars, skip)
    }

    /// Remainder of the S-polynomial of divisors `i` and `j` — that is,
    /// `reduce(&s_polynomial(g_i, g_j, order), gs, order)` — with the
    /// S-polynomial assembled directly in the working representation.
    pub(crate) fn reduce_s_pair(&self, i: usize, j: usize) -> GbPoly {
        let order = self.order;
        let (Some(f), Some(g)) = (&self.divs[i], &self.divs[j]) else {
            // An S-polynomial involving a zero polynomial is zero.
            return GbPoly::zero(self.n_vars[i]);
        };
        let w = f.lead.len();
        if f.width != Some(w) || !self.all_width(w, None) {
            let sp = s_polynomial(&f.to_poly(), &g.to_poly(), order);
            return reduce_reference(&sp, &self.polys_except(None), order);
        }
        let lcm: Vec<u32> = f.lead.iter().zip(&g.lead).map(|(a, b)| *a.max(b)).collect();
        let shift_f: Vec<u32> = lcm.iter().zip(&f.lead).map(|(l, a)| l - a).collect();
        let shift_g: Vec<u32> = lcm.iter().zip(&g.lead).map(|(l, b)| l - b).collect();
        let inv_f = (!f.lc_is_one).then(|| Rational::from(1) / &f.lc);
        let inv_g = (!g.lc_is_one).then(|| Rational::from(1) / &g.lc);
        // S = x^sf·f/lc(f) − x^sg·g/lc(g).  The two leading terms cancel exactly
        // and every tail term lies strictly below the lcm, so only tails matter.
        let mut p = BTreeMap::new();
        let mut key = Vec::with_capacity(w + 1);
        let mut tmp = Rational::new();
        for (e, c) in &f.tail {
            encode_shifted(order, e, &shift_f, &mut key);
            accumulate(&mut p, &key, c, inv_f.as_ref(), false, &mut tmp);
        }
        for (e, c) in &g.tail {
            encode_shifted(order, e, &shift_g, &mut key);
            accumulate(&mut p, &key, c, inv_g.as_ref(), true, &mut tmp);
        }
        self.reduce_work(p, f.n_vars, None)
    }

    /// A non-zero rational multiple of [`Self::reduce`]`(f, skip)`, computed
    /// fraction-free over ℤ.  Callers that normalise the result (`make_monic`)
    /// get exactly what they would from `reduce`.
    ///
    /// Scaling the dividend by a non-zero constant never changes which terms
    /// are present, so every step picks the same divisor as the rational
    /// algorithm and the two runs stay proportional throughout.
    pub(crate) fn reduce_scaled(&self, f: &GbPoly, skip: Option<usize>) -> GbPoly {
        let order = self.order;
        let Some(first) = f.terms.keys().next() else {
            return GbPoly::zero(f.n_vars);
        };
        let w = first.len();
        let regular = f.terms.iter().all(|(e, c)| e.len() == w && *c != 0);
        if !self.fraction_free || !regular || !self.all_width(w, skip) {
            return self.reduce(f, skip);
        }
        let ints = primitive_integer_multiple(f.terms.values());
        let zero_shift = vec![0u32; w];
        let mut key = Vec::with_capacity(w + 1);
        let mut p = BTreeMap::new();
        for (e, c) in f.terms.keys().zip(ints) {
            encode_shifted(order, e, &zero_shift, &mut key);
            p.insert(key.clone(), c);
        }
        self.reduce_work_int(p, f.n_vars, skip)
    }

    /// A non-zero rational multiple of [`Self::reduce_s_pair`]`(i, j)`,
    /// computed fraction-free over ℤ (see [`Self::reduce_scaled`]).
    pub(crate) fn reduce_s_pair_scaled(&self, i: usize, j: usize) -> GbPoly {
        let order = self.order;
        let (Some(f), Some(g)) = (&self.divs[i], &self.divs[j]) else {
            return GbPoly::zero(self.n_vars[i]);
        };
        let w = f.lead.len();
        let (Some(fi), Some(gi)) = (&f.int, &g.int) else {
            return self.reduce_s_pair(i, j);
        };
        if f.width != Some(w) || !self.all_width(w, None) {
            return self.reduce_s_pair(i, j);
        }
        let lcm: Vec<u32> = f.lead.iter().zip(&g.lead).map(|(a, b)| *a.max(b)).collect();
        let shift_f: Vec<u32> = lcm.iter().zip(&f.lead).map(|(l, a)| l - a).collect();
        let shift_g: Vec<u32> = lcm.iter().zip(&g.lead).map(|(l, b)| l - b).collect();
        // (lc(G)/h)·x^sf·F − (lc(F)/h)·x^sg·G with h = gcd(lc(F), lc(G)): the
        // leading terms cancel, and the result is a multiple of S(f, g).
        let h = Integer::from(fi.lc.gcd_ref(&gi.lc));
        let cf = Integer::from(gi.lc.div_exact_ref(&h));
        let cg = Integer::from(fi.lc.div_exact_ref(&h));
        let mut p = BTreeMap::new();
        let mut key = Vec::with_capacity(w + 1);
        let mut tmp = Integer::new();
        for ((e, _), c) in f.tail.iter().zip(&fi.tail) {
            encode_shifted(order, e, &shift_f, &mut key);
            accumulate_int(&mut p, &key, c, &cf, false, &mut tmp);
        }
        for ((e, _), c) in g.tail.iter().zip(&gi.tail) {
            encode_shifted(order, e, &shift_g, &mut key);
            accumulate_int(&mut p, &key, c, &cg, true, &mut tmp);
        }
        self.reduce_work_int(p, f.n_vars, None)
    }

    /// Fraction-free division loop: `P ← (lc(G)/h)·P − (lc(P)/h)·x^s·G`, with
    /// `h = gcd(lc(P), lc(G))`; the remainder terms are scaled along with `P`
    /// so that the pair stays a multiple of the rational computation.  The
    /// integer content is divided out periodically to curb coefficient growth.
    fn reduce_work_int(
        &self,
        mut p: BTreeMap<Vec<u32>, Integer>,
        n_vars: usize,
        skip: Option<usize>,
    ) -> GbPoly {
        let order = self.order;
        let is_graded = order.is_graded();
        let mut rem: Vec<(Vec<u32>, Integer)> = Vec::new();
        let mut last_divisor = 0usize;
        let mut lt_exp: Vec<u32> = Vec::new();
        let mut shift: Vec<u32> = Vec::new();
        let mut key: Vec<u32> = Vec::new();
        let mut tmp = Integer::new();
        let mut h = Integer::new();
        let mut mp = Integer::new();
        let mut scalings = 0u32;

        while let Some((lt_key, mut mg)) = p.pop_last() {
            decode(order, &lt_key, &mut lt_exp);
            let lt_deg = if is_graded { lt_key[0] } else { 0 };
            let Some((idx, d)) = self.find_divisor(&lt_exp, lt_deg, last_divisor, skip) else {
                rem.push((lt_exp.clone(), mg));
                continue;
            };
            let di = d.int.as_ref().expect("fraction-free divisor");
            h.assign(mg.gcd_ref(&di.lc));
            mp.assign(di.lc.div_exact_ref(&h));
            mg.div_exact_mut(&h);
            if mp != 1 {
                for v in p.values_mut() {
                    *v *= &mp;
                }
                for (_, v) in &mut rem {
                    *v *= &mp;
                }
                scalings += 1;
            }
            shift.clear();
            shift.extend(lt_exp.iter().zip(&d.lead).map(|(a, b)| a - b));
            for ((e, _), c) in d.tail.iter().zip(&di.tail) {
                encode_shifted(order, e, &shift, &mut key);
                accumulate_int(&mut p, &key, c, &mg, true, &mut tmp);
            }
            last_divisor = idx;
            if scalings >= CONTENT_EVERY {
                scalings = 0;
                remove_content(&mut p, &mut rem, &mut tmp);
            }
        }

        GbPoly {
            terms: rem
                .into_iter()
                .map(|(e, c)| (e, Rational::from(c)))
                .collect(),
            n_vars,
        }
    }

    /// The first divisor, in cyclic order from `last` and ignoring `skip`,
    /// whose leading monomial divides `lt_exp` (of total degree `lt_deg`,
    /// consulted only for graded orders).
    #[inline]
    fn find_divisor(
        &self,
        lt_exp: &[u32],
        lt_deg: u32,
        last: usize,
        skip: Option<usize>,
    ) -> Option<(usize, &Divisor)> {
        let is_graded = self.order.is_graded();
        let n_divs = self.divs.len();
        for offset in 0..n_divs {
            let idx = (last + offset) % n_divs;
            if Some(idx) == skip {
                continue;
            }
            let Some(d) = &self.divs[idx] else {
                continue;
            };
            if is_graded && d.lead_deg > lt_deg {
                continue;
            }
            if lt_exp.iter().zip(&d.lead).all(|(a, b)| a >= b) {
                return Some((idx, d));
            }
        }
        None
    }

    /// The division loop proper, on an encoded working map.
    fn reduce_work(
        &self,
        mut p: BTreeMap<Vec<u32>, Rational>,
        n_vars: usize,
        skip: Option<usize>,
    ) -> GbPoly {
        let order = self.order;
        let is_graded = order.is_graded();
        let mut rem: Vec<(Vec<u32>, Rational)> = Vec::new();
        let mut last_divisor = 0usize;
        let mut lt_exp: Vec<u32> = Vec::new();
        let mut shift: Vec<u32> = Vec::new();
        let mut key: Vec<u32> = Vec::new();
        let mut tmp = Rational::new();

        while let Some((lt_key, mut q)) = p.pop_last() {
            decode(order, &lt_key, &mut lt_exp);
            let lt_deg = if is_graded { lt_key[0] } else { 0 };
            let Some((idx, d)) = self.find_divisor(&lt_exp, lt_deg, last_divisor, skip) else {
                // No divisor: the leading term moves to the remainder.
                rem.push((lt_exp.clone(), q));
                continue;
            };
            // p ← p − (lt(p)/lt(d))·x^shift·d; the leading terms cancel, and
            // p's has already been popped.
            if !d.lc_is_one {
                q /= &d.lc;
            }
            shift.clear();
            shift.extend(lt_exp.iter().zip(&d.lead).map(|(a, b)| a - b));
            for (e, c) in &d.tail {
                encode_shifted(order, e, &shift, &mut key);
                accumulate(&mut p, &key, c, Some(&q), true, &mut tmp);
            }
            last_divisor = idx;
        }

        GbPoly {
            terms: rem.into_iter().collect(),
            n_vars,
        }
    }
}

/// The straightforward division algorithm on [`GbPoly`] values: the
/// implementation [`reduce`] had before the in-place engine.  Used for inputs
/// outside the engine's model (mixed-length exponent vectors, explicitly stored
/// zero coefficients) and as the oracle of the differential tests.
pub(crate) fn reduce_reference(f: &GbPoly, gs: &[GbPoly], order: MonomialOrder) -> GbPoly {
    let mut p = f.clone();
    let mut r = GbPoly::zero(f.n_vars);

    // Cache the index of the last successful divisor.  When consecutive leading
    // terms of `p` are divisible by the same basis element, this avoids scanning
    // from the start of `gs` on every step.
    let mut last_divisor: usize = 0;

    let is_graded = order.is_graded();

    'outer: while !p.is_zero() {
        let (lt_exp, lt_coeff) = match p.leading_term(order) {
            Some((e, c)) => (e.clone(), c.clone()),
            None => break,
        };
        // For graded orders, precompute the total degree of the current leading term.
        // Any basis element whose LM has higher total degree cannot divide it.
        let lt_deg: u32 = if is_graded {
            total_degree_or_panic(&lt_exp)
        } else {
            0
        };

        // Try gs[last_divisor], then wrap around through all of gs.
        for offset in 0..gs.len() {
            let idx = (last_divisor + offset) % gs.len();
            let g = &gs[idx];
            if let Some((lg_exp, lg_coeff)) = g.leading_term(order) {
                if is_graded && total_degree_or_panic(lg_exp) > lt_deg {
                    continue;
                }
                if lt_exp.len() == lg_exp.len()
                    && lt_exp.iter().zip(lg_exp.iter()).all(|(a, b)| a >= b)
                {
                    let shift: Vec<u32> = lt_exp
                        .iter()
                        .zip(lg_exp.iter())
                        .map(|(a, b)| a - b)
                        .collect();
                    let coeff = rug::Rational::from(&lt_coeff / lg_coeff);
                    let subtrahend = g.mul_monomial(&shift, &coeff);
                    p = p.sub(&subtrahend);
                    last_divisor = idx;
                    continue 'outer;
                }
            }
        }

        // No divisor found — move leading term to remainder.
        let lt = GbPoly::monomial(lt_exp.clone(), lt_coeff.clone());
        r = r.add(&lt);
        let mut p_terms = p.terms.clone();
        p_terms.remove(&lt_exp);
        p.terms = p_terms;
    }

    r
}

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

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

    const ORDERS: [MonomialOrder; 3] = [
        MonomialOrder::Lex,
        MonomialOrder::GrLex,
        MonomialOrder::GRevLex,
    ];

    fn assert_same(a: &GbPoly, b: &GbPoly) {
        assert_eq!(a.terms, b.terms);
        assert_eq!(a.n_vars, b.n_vars);
    }

    /// A polynomial in `n` variables with ≤ `max_terms` terms of degree ≤ 3 and
    /// small rational (not necessarily integral) coefficients.
    fn arb_poly(n: usize, max_terms: usize) -> impl Strategy<Value = GbPoly> {
        let term = (proptest::collection::vec(0u32..=3, n), -4i64..=4, 1i64..=3);
        proptest::collection::vec(term, 0..=max_terms).prop_map(move |ts| {
            let mut p = GbPoly::zero(n);
            for (e, num, den) in ts {
                p = p.add(&GbPoly::monomial(e, rat(num, den)));
            }
            p
        })
    }

    fn arb_problem() -> impl Strategy<Value = (GbPoly, Vec<GbPoly>, usize)> {
        (1usize..=3).prop_flat_map(|n| {
            (
                arb_poly(n, 8),
                proptest::collection::vec(arb_poly(n, 4), 0..=4),
                0usize..3,
            )
        })
    }

    #[test]
    fn encoding_preserves_every_order() {
        let mut exps = vec![];
        for a in 0..3u32 {
            for b in 0..3u32 {
                for c in 0..3u32 {
                    exps.push(vec![a, b, c]);
                }
            }
        }
        let zero = [0u32; 3];
        for order in ORDERS {
            for x in &exps {
                let (mut kx, mut back) = (vec![], vec![]);
                encode_shifted(order, x, &zero, &mut kx);
                decode(order, &kx, &mut back);
                assert_eq!(&back, x);
                for y in &exps {
                    let mut ky = vec![];
                    encode_shifted(order, y, &zero, &mut ky);
                    assert_eq!(kx.cmp(&ky), order.cmp(x, y), "{order:?} {x:?} {y:?}");
                }
            }
        }
    }

    #[test]
    fn irregular_inputs_take_the_reference_path() {
        // Mixed exponent widths among the divisors.
        let f = GbPoly::monomial(vec![2, 1], rat(3, 2));
        let gs = vec![
            GbPoly::monomial(vec![1], rat(1, 1)),
            GbPoly {
                terms: [(vec![1, 0], rat(2, 1)), (vec![0, 1], rat(1, 1))]
                    .into_iter()
                    .collect(),
                n_vars: 2,
            },
        ];
        for order in ORDERS {
            assert_same(&reduce(&f, &gs, order), &reduce_reference(&f, &gs, order));
        }
        // An explicitly stored zero coefficient in the dividend.
        let f = GbPoly {
            terms: [(vec![0, 3], rat(0, 1)), (vec![2, 0], rat(1, 1))]
                .into_iter()
                .collect(),
            n_vars: 2,
        };
        let gs = vec![GbPoly::monomial(vec![1, 0], rat(1, 1))];
        for order in ORDERS {
            assert_same(&reduce(&f, &gs, order), &reduce_reference(&f, &gs, order));
        }
    }

    proptest! {
        #![proptest_config(ProptestConfig::with_cases(256))]

        /// `reduce` equals the reference division algorithm term for term,
        /// for arbitrary (non-Gröbner) divisor lists.
        #[test]
        fn reduce_matches_reference((f, gs, ord) in arb_problem()) {
            let order = ORDERS[ord];
            let new = reduce(&f, &gs, order);
            let old = reduce_reference(&f, &gs, order);
            prop_assert_eq!(&new.terms, &old.terms);
            prop_assert_eq!(new.n_vars, old.n_vars);
        }

        /// Skipping index `k` equals dividing by the list with `gs[k]` removed.
        #[test]
        fn reduce_skip_matches_removal((f, gs, ord) in arb_problem(), k in 0usize..4) {
            prop_assume!(k < gs.len());
            let order = ORDERS[ord];
            let new = Divisors::new(&gs, order).reduce(&f, Some(k));
            let mut others = gs.clone();
            others.remove(k);
            let old = reduce_reference(&f, &others, order);
            prop_assert_eq!(&new.terms, &old.terms);
        }

        /// The fused S-pair path equals `reduce_reference(s_polynomial(..))`.
        #[test]
        fn reduce_s_pair_matches_reference((_f, gs, ord) in arb_problem(), i in 0usize..12, j in 0usize..12) {
            prop_assume!(gs.len() >= 2);
            let (i, j) = (i % gs.len(), j % gs.len());
            prop_assume!(i != j);
            let order = ORDERS[ord];
            let new = Divisors::new(&gs, order).reduce_s_pair(i, j);
            let sp = s_polynomial(&gs[i], &gs[j], order);
            let old = reduce_reference(&sp, &gs, order);
            prop_assert_eq!(&new.terms, &old.terms);
            prop_assert_eq!(new.n_vars, old.n_vars);
        }

        /// The fraction-free remainder is a non-zero multiple of the rational
        /// one: identical once made monic (and zero exactly when it is).
        #[test]
        fn reduce_scaled_is_a_multiple((f, gs, ord) in arb_problem(), k in 0usize..5) {
            let order = ORDERS[ord];
            let skip = (k < gs.len()).then_some(k);
            let divs = Divisors::new_fraction_free(&gs, order);
            let scaled = divs.reduce_scaled(&f, skip);
            let exact = divs.reduce(&f, skip);
            prop_assert_eq!(scaled.n_vars, exact.n_vars);
            prop_assert_eq!(&scaled.make_monic(order).terms, &exact.make_monic(order).terms);
        }

        #[test]
        fn reduce_s_pair_scaled_is_a_multiple((_f, gs, ord) in arb_problem(), i in 0usize..12, j in 0usize..12) {
            prop_assume!(gs.len() >= 2);
            let (i, j) = (i % gs.len(), j % gs.len());
            prop_assume!(i != j);
            let order = ORDERS[ord];
            let divs = Divisors::new_fraction_free(&gs, order);
            let scaled = divs.reduce_s_pair_scaled(i, j);
            let exact = divs.reduce_s_pair(i, j);
            prop_assert_eq!(scaled.n_vars, exact.n_vars);
            prop_assert_eq!(&scaled.make_monic(order).terms, &exact.make_monic(order).terms);
        }
    }

    #[test]
    fn reduce_to_remainder() {
        // Reduce x^2 by [x^2 - 1] → remainder = 1
        let f = GbPoly::monomial(vec![2, 0], rat(1, 1));
        let g = GbPoly {
            terms: [(vec![2, 0], rat(1, 1)), (vec![0, 0], rat(-1, 1))]
                .into_iter()
                .collect(),
            n_vars: 2,
        };
        let r = reduce(&f, &[g], MonomialOrder::Lex);
        // f = x^2, g = x^2 - 1 → f - g = 1
        assert!(!r.is_zero());
    }

    #[test]
    fn s_polynomial_example() {
        // S(x^2 - 1, x - 1) under Lex over 1 variable
        // lcm(x^2, x) = x^2, shift_f = 1, shift_g = x
        // S = (x^2 - 1) - x*(x-1) = x^2 - 1 - x^2 + x = x - 1
        let f = GbPoly {
            terms: [(vec![2], rat(1, 1)), (vec![0], rat(-1, 1))]
                .into_iter()
                .collect(),
            n_vars: 1,
        };
        let g = GbPoly {
            terms: [(vec![1], rat(1, 1)), (vec![0], rat(-1, 1))]
                .into_iter()
                .collect(),
            n_vars: 1,
        };
        let s = s_polynomial(&f, &g, MonomialOrder::Lex);
        assert!(!s.is_zero());
    }
}