alkahest-cas 2.0.3

//! F4-style Gröbner basis computation over ℚ.
//!
//! Implements Buchberger's algorithm with:
//! - Parallel S-polynomial reduction via Rayon (when `parallel` feature is enabled)
//! - Buchberger's product criterion to prune S-pairs
//! - Interreduction of the final basis

use crate::poly::groebner::ideal::GbPoly;
use crate::poly::groebner::monomial_order::MonomialOrder;
use crate::poly::groebner::reduce::{reduce, s_polynomial};

/// Check Buchberger's product criterion: if lcm(lm(f), lm(g)) = lm(f)*lm(g),
/// i.e., leading monomials are coprime, the S-polynomial reduces to 0.
fn product_criterion(f: &GbPoly, g: &GbPoly, order: MonomialOrder) -> bool {
    let lf = match f.leading_exp(order) {
        Some(e) => e,
        None => return true,
    };
    let lg = match g.leading_exp(order) {
        Some(e) => e,
        None => return true,
    };
    lf.iter().zip(lg.iter()).all(|(&a, &b)| a == 0 || b == 0)
}

/// Compute a Gröbner basis for the ideal generated by `generators` under `order`.
///
/// Uses Buchberger's algorithm with the product criterion for pruning.
/// S-polynomial reductions are parallelized when the `parallel` feature is enabled.
pub fn compute_groebner_basis(generators: Vec<GbPoly>, order: MonomialOrder) -> Vec<GbPoly> {
    let mut basis: Vec<GbPoly> = generators
        .into_iter()
        .filter(|g| !g.is_zero())
        .map(|g| g.make_monic(order))
        .collect();

    if basis.is_empty() {
        return basis;
    }

    let mut pairs: Vec<(usize, usize)> = vec![];
    for i in 0..basis.len() {
        for j in (i + 1)..basis.len() {
            if !product_criterion(&basis[i], &basis[j], order) {
                pairs.push((i, j));
            }
        }
    }

    while !pairs.is_empty() {
        let s_polys: Vec<GbPoly> = {
            #[cfg(feature = "parallel")]
            {
                use rayon::prelude::*;
                pairs
                    .par_iter()
                    .map(|&(i, j)| s_polynomial(&basis[i], &basis[j], order))
                    .collect()
            }
            #[cfg(not(feature = "parallel"))]
            {
                pairs
                    .iter()
                    .map(|&(i, j)| s_polynomial(&basis[i], &basis[j], order))
                    .collect()
            }
        };

        let reduced: Vec<GbPoly> = {
            #[cfg(feature = "parallel")]
            {
                use rayon::prelude::*;
                s_polys
                    .par_iter()
                    .map(|sp| reduce(sp, &basis, order))
                    .collect()
            }
            #[cfg(not(feature = "parallel"))]
            {
                s_polys.iter().map(|sp| reduce(sp, &basis, order)).collect()
            }
        };

        let new_start = basis.len();
        for r in reduced {
            if !r.is_zero() {
                basis.push(r.make_monic(order));
            }
        }

        pairs.clear();
        for i in 0..basis.len() {
            for j in (i + 1)..basis.len() {
                if !product_criterion(&basis[i], &basis[j], order) {
                    if i >= new_start || j >= new_start {
                        pairs.push((i, j));
                    }
                }
            }
        }
    }

    interreduce(basis, order)
}

/// Interreduce a Gröbner basis: reduce each element by all others and remove
/// elements whose leading term is divisible by another's.
pub(crate) fn interreduce(mut basis: Vec<GbPoly>, order: MonomialOrder) -> Vec<GbPoly> {
    let mut i = 0;
    while i < basis.len() {
        let others: Vec<GbPoly> = basis
            .iter()
            .enumerate()
            .filter(|&(j, _)| j != i)
            .map(|(_, g)| g.clone())
            .collect();
        let reduced = reduce(&basis[i], &others, order);
        if reduced.is_zero() {
            basis.remove(i);
        } else {
            basis[i] = reduced.make_monic(order);
            i += 1;
        }
    }
    basis
}

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

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

    fn poly(terms: &[(&[u32], i64)]) -> GbPoly {
        let n_vars = terms.first().map(|(e, _)| e.len()).unwrap_or(1);
        GbPoly {
            terms: terms
                .iter()
                .map(|(e, c)| (e.to_vec(), rat(*c, 1)))
                .collect(),
            n_vars,
        }
    }

    #[test]
    fn groebner_x_squared_minus_1_and_x_minus_1() {
        // Ideal: (x^2 - 1, x - 1) over 1 variable
        // Gröbner basis under Lex should be {x - 1}
        let f = poly(&[(&[2], 1), (&[0], -1)]); // x^2 - 1
        let g = poly(&[(&[1], 1), (&[0], -1)]); // x - 1
        let basis = compute_groebner_basis(vec![f, g], MonomialOrder::Lex);
        assert_eq!(basis.len(), 1);
        assert!(basis[0].terms.contains_key(&vec![1]));
    }

    #[test]
    fn groebner_trivial_ideal() {
        // Ideal: (x, y) over 2 variables — basis should be {x, y}
        let f = poly(&[(&[1, 0], 1)]); // x
        let g = poly(&[(&[0, 1], 1)]); // y
        let basis = compute_groebner_basis(vec![f, g], MonomialOrder::GRevLex);
        assert_eq!(basis.len(), 2);
    }

    #[test]
    fn groebner_linear_system() {
        // Ideal: (x + y - 1, x - y) over 2 variables
        // Solutions: x = 1/2, y = 1/2
        let f = poly(&[(&[1, 0], 1), (&[0, 1], 1), (&[0, 0], -1)]); // x + y - 1
        let g = poly(&[(&[1, 0], 1), (&[0, 1], -1)]); // x - y
        let basis = compute_groebner_basis(vec![f, g], MonomialOrder::Lex);
        assert!(!basis.is_empty());
        let orig_f = poly(&[(&[1, 0], 1), (&[0, 1], 1), (&[0, 0], -1)]);
        let orig_g = poly(&[(&[1, 0], 1), (&[0, 1], -1)]);
        let r_f = reduce(&orig_f, &basis, MonomialOrder::Lex);
        let r_g = reduce(&orig_g, &basis, MonomialOrder::Lex);
        assert!(r_f.is_zero(), "f does not reduce to 0 mod basis: {:?}", r_f);
        assert!(r_g.is_zero(), "g does not reduce to 0 mod basis: {:?}", r_g);
    }

    #[test]
    fn contains_check() {
        // 0 is in every ideal
        let f = poly(&[(&[1, 0], 1)]); // x
        let basis = compute_groebner_basis(vec![f], MonomialOrder::Lex);
        let zero = GbPoly::zero(2);
        let r = reduce(&zero, &basis, MonomialOrder::Lex);
        assert!(r.is_zero());
    }
}