alkahest-cas 3.10.0

//! Buchberger's algorithm for Gröbner basis computation over ℚ.
//!
//! Implements the sequential Buchberger algorithm with:
//! - Gebauer-Möller criteria M and F to prune S-pairs
//! - Sugar selection strategy: process pair with minimum sugar degree first, break ties by lcm degree
//! - Incremental basis update: each new element is added before selecting the next pair
//!
//! The pair machinery itself is in `super::pairs` (private — internal to this
//! module), shared with the `Q(params)` engine in [`super::parametric`].
//!
//! Reference: Becker & Weispfenning (1993) "Gröbner Bases", Algorithm 6.5 (GROEBNERNEWS2),
//! Gebauer & Möller (1988) "On an Installation of Buchberger's Algorithm",
//! and Giovini et al. (1991) "One Sugar Cube, Please" for the sugar selection strategy.

use std::collections::BinaryHeap;

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

// ---------------------------------------------------------------------------
// Main algorithm
// ---------------------------------------------------------------------------

/// Compute a Gröbner basis for the ideal generated by `generators` under `order`.
///
/// Uses sequential Buchberger with Gebauer-Möller pair management and
/// normal selection (process the pair with minimum lcm degree first).
pub fn compute_buchberger_basis(generators: Vec<GbPoly>, order: MonomialOrder) -> Vec<GbPoly> {
    let initial: Vec<GbPoly> = generators
        .into_iter()
        .filter(|g| !g.is_zero())
        .map(|g| g.make_monic(order))
        .collect();

    if initial.is_empty() {
        return initial;
    }

    let mut basis: Vec<GbPoly> = Vec::with_capacity(initial.len() * 2);
    let mut basis_sugar: Vec<u32> = Vec::with_capacity(initial.len() * 2);
    // Leading exponents, kept in step with `basis`.  Every element pushed here
    // is non-zero (the generators were filtered and reductions are only added
    // when they do not vanish), so `leading_exp` always yields a value.
    let mut basis_lead: Vec<Vec<u32>> = Vec::with_capacity(initial.len() * 2);
    let mut pair_vec: Vec<CriticalPair> = Vec::new();

    // Add initial generators one by one, applying GM update after each.
    for gen in initial {
        let sugar = gen.sugar();
        let Some(lead) = gen.leading_exp(order) else {
            continue;
        };
        let new_idx = basis.len();
        basis.push(gen);
        basis_sugar.push(sugar);
        basis_lead.push(lead);
        update_pairs(&basis_lead, &basis_sugar, &mut pair_vec, new_idx);
    }

    // Build min-heap (CriticalPair::Ord is reversed for min-heap behaviour).
    let mut heap: BinaryHeap<CriticalPair> = BinaryHeap::from(pair_vec);

    while let Some(pair) = heap.pop() {
        let sp = s_polynomial(&basis[pair.i], &basis[pair.j], order);
        let r = reduce(&sp, &basis, order);

        if !r.is_zero() {
            let r = r.make_monic(order);
            let sugar = r.sugar();
            let Some(lead) = r.leading_exp(order) else {
                continue;
            };
            let new_idx = basis.len();
            basis.push(r);
            basis_sugar.push(sugar);
            basis_lead.push(lead);

            // Flatten heap → apply GM update → rebuild heap.
            let mut pv: Vec<CriticalPair> = heap.into_vec();
            update_pairs(&basis_lead, &basis_sugar, &mut pv, new_idx);
            heap = BinaryHeap::from(pv);
        }
    }

    interreduce(basis, order)
}

// ---------------------------------------------------------------------------
// Interreduction
// ---------------------------------------------------------------------------

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

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[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() {
        let f = poly(&[(&[2], 1), (&[0], -1)]);
        let g = poly(&[(&[1], 1), (&[0], -1)]);
        let basis = compute_buchberger_basis(vec![f, g], MonomialOrder::Lex);
        assert_eq!(basis.len(), 1);
        assert!(basis[0].terms.contains_key(&vec![1]));
    }

    #[test]
    fn groebner_trivial_ideal() {
        let f = poly(&[(&[1, 0], 1)]);
        let g = poly(&[(&[0, 1], 1)]);
        let basis = compute_buchberger_basis(vec![f, g], MonomialOrder::GRevLex);
        assert_eq!(basis.len(), 2);
    }

    #[test]
    fn groebner_linear_system() {
        let f = poly(&[(&[1, 0], 1), (&[0, 1], 1), (&[0, 0], -1)]);
        let g = poly(&[(&[1, 0], 1), (&[0, 1], -1)]);
        let basis = compute_buchberger_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() {
        let f = poly(&[(&[1, 0], 1)]);
        let basis = compute_buchberger_basis(vec![f], MonomialOrder::Lex);
        let zero = GbPoly::zero(2);
        let r = reduce(&zero, &basis, MonomialOrder::Lex);
        assert!(r.is_zero());
    }

    #[test]
    fn circle_parabola_grevlex() {
        // x^2 + y^2 - 4, y - x^2 + 1
        let x2_y2_m4 = poly(&[(&[2, 0], 1), (&[0, 2], 1), (&[0, 0], -4)]);
        let y_mx2_p1 = poly(&[(&[0, 1], 1), (&[2, 0], -1), (&[0, 0], 1)]);
        let basis =
            compute_buchberger_basis(vec![x2_y2_m4, y_mx2_p1.clone()], MonomialOrder::GRevLex);
        // Verify generators reduce to 0
        assert!(!basis.is_empty());
        let r = reduce(&y_mx2_p1, &basis, MonomialOrder::GRevLex);
        assert!(r.is_zero());
    }
}