alkahest-cas 2.0.3

//! 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);

    '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 g in gs {
            if let Some((lg_exp, lg_coeff)) = g.leading_term(order) {
                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);
                    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(rug::Rational::from(1) / &lf_coeff);
    let coeff_g = rug::Rational::from(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());
    }
}