alkahest-cas 3.5.1

//! Polynomial reduction (division algorithm) for Gröbner bases.

use crate::poly::groebner::ideal::GbPoly;
use crate::poly::groebner::monomial_order::MonomialOrder;

/// 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 {
    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 { lt_exp.iter().sum() } 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) {
                // Degree-skip: for graded orders, LM(g) with total degree > lt_deg
                // cannot divide lt, so skip without checking all exponents.
                if is_graded && lg_exp.iter().sum::<u32>() > 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
}

/// 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)
}

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

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

    #[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());
    }
}