add buchberger algorithm

src/poly/buchberger.rs

@@ -0,0 +1,74 @@
+use super::flat::Poly;
+use super::var::Var;
+
+/// Computes a Gröbner basis for the ideal generated by `generators` using
+/// Buchberger's algorithm under lex order.
+///
+/// The returned basis spans the same ideal as the input and satisfies
+/// Buchberger's criterion: every S-polynomial of a pair in the basis
+/// reduces to zero modulo the basis.
+pub fn groebner_basis<V: Var>(generators: Vec<Poly<V>>) -> Vec<Poly<V>> {
+    let mut g: Vec<Poly<V>> = generators.into_iter().filter(|p| !p.is_zero()).collect();
+
+    let mut i = 0;
+    while i < g.len() {
+        let mut j = i + 1;
+        while j < g.len() {
+            let s = g[i].s_poly(&g[j]);
+            let r = reduce(&s, &g);
+            if !r.is_zero() {
+                g.push(r);
+            }
+            j += 1;
+        }
+        i += 1;
+    }
+
+    g
+}
+
+/// Checks whether `basis` satisfies Buchberger's criterion under lex order:
+/// the S-polynomial of every pair reduces to zero modulo the basis.
+///
+/// Returns `true` iff `basis` is a Gröbner basis for the ideal it generates.
+pub fn is_groebner_basis<V: Var>(basis: &[Poly<V>]) -> bool {
+    for i in 0..basis.len() {
+        for j in (i + 1)..basis.len() {
+            let s = basis[i].s_poly(&basis[j]);
+            if !reduce(&s, basis).is_zero() {
+                return false;
+            }
+        }
+    }
+    true
+}
+
+/// Reduces `f` modulo `basis` until no leading term of `f` is divisible
+/// by the leading monomial of any element in `basis`.
+///
+/// Uses the pseudo-division remainder and repeats until stable.
+fn reduce<V: Var>(f: &Poly<V>, basis: &[Poly<V>]) -> Poly<V> {
+    let mut p = f.clone();
+
+    'outer: loop {
+        let Some((lm_p, _)) = p.leading_term_lex() else {
+            break;
+        };
+
+        for g in basis {
+            let Some((lm_g, _)) = g.leading_term_lex() else {
+                continue;
+            };
+
+            if lm_p.contains(&lm_g) {
+                let (_, _, r) = p.clone().div_rem(g);
+                p = r;
+                continue 'outer;
+            }
+        }
+
+        break;
+    }
+
+    p
+}

src/poly/flat.rs

@@ -11,10 +11,18 @@ pub fn lex_cmp<V: Var>(a: &Mono<V>, b: &Mono<V>) -> Ordering {
         match (a_it.peek(), b_it.peek()) {
             (None, None) => return Ordering::Equal,
             (Some((_, a_exp)), None) => {
-                return if *a_exp > 0 { Ordering::Greater } else { Ordering::Equal };
+                return if *a_exp > 0 {
+                    Ordering::Greater
+                } else {
+                    Ordering::Equal
+                };
             }
             (None, Some((_, b_exp))) => {
-                return if *b_exp > 0 { Ordering::Less } else { Ordering::Equal };
+                return if *b_exp > 0 {
+                    Ordering::Less
+                } else {
+                    Ordering::Equal
+                };
             }
             (Some((a_var, a_exp)), Some((b_var, b_exp))) => {
                 if a_var < b_var {

src/poly/mod.rs

@@ -1,3 +1,4 @@
+pub mod buchberger;
 pub mod flat;
 pub mod fmt;
 pub mod ops;

src/poly/ops.rs

@@ -168,9 +168,7 @@ impl<V: Var> Poly<V> {
     ///
     /// Panics if `divisor` is zero.
     pub fn div_rem(self, divisor: &Poly<V>) -> (u32, Poly<V>, Poly<V>) {
-        let (lt_g_mono, lt_g_coeff) = divisor
+        let (lt_g_mono, lt_g_coeff) = divisor.leading_term_lex().expect("divisor must be nonzero");
-            .leading_term_lex()
-            .expect("divisor must be nonzero");
 
         let mut p = self;
         let mut q = Poly::default();

src/poly/tests.rs

@@ -1,4 +1,5 @@
-use super::flat::{lex_cmp, Mono, Poly};
+use super::buchberger::{groebner_basis, is_groebner_basis};
+use super::flat::{Mono, Poly, lex_cmp};
 use super::var::StaticVar;
 
 #[test]
@@ -175,15 +176,21 @@ fn test_s_poly() {
     let f: Poly<StaticVar> = [
         (1i32, Mono::from([("x", 2u32)])),
         (1i32, Mono::from([("y", 1u32)])),
-    ].into_iter().collect();
+    ]
+    .into_iter()
+    .collect();
     let g: Poly<StaticVar> = [
         (1i32, Mono::from([("x", 1u32), ("y", 1u32)])),
         (1i32, Mono::from([("z", 1u32)])),
-    ].into_iter().collect();
+    ]
+    .into_iter()
+    .collect();
     let expected: Poly<StaticVar> = [
         (1i32, Mono::from([("y", 2u32)])),
         (-1i32, Mono::from([("x", 1u32), ("z", 1u32)])),
-    ].into_iter().collect();
+    ]
+    .into_iter()
+    .collect();
     assert_eq!(f.s_poly(&g), expected);
 
     // f = 2x + y, g = 3x + z  (same LM=x, d=gcd(2,3)=1)
@@ -201,10 +208,18 @@ fn test_s_poly() {
 }
 
 fn make_const_poly(c: i32) -> Poly<StaticVar> {
-    Poly { mono: [(Mono { term: vec![] }, c)].into_iter().collect() }
+    Poly {
+        mono: [(Mono { term: vec![] }, c)].into_iter().collect(),
+    }
 }
 
-fn verify_div_rem(f: Poly<StaticVar>, g: &Poly<StaticVar>, d: u32, q: Poly<StaticVar>, r: Poly<StaticVar>) {
+fn verify_div_rem(
+    f: Poly<StaticVar>,
+    g: &Poly<StaticVar>,
+    d: u32,
+    q: Poly<StaticVar>,
+    r: Poly<StaticVar>,
+) {
     // lc(g)^d * f == q * g + r
     let (_, lc_g) = g.leading_term_lex().unwrap();
     let lhs = lc_g.pow(d) * f;
@@ -228,7 +243,9 @@ fn test_div_rem() {
     let f: Poly<StaticVar> = [
         (1i32, Mono::from([("x", 3u32)])),
         (1i32, Mono::from([("x", 2u32), ("y", 1u32)])),
-    ].into_iter().collect();
+    ]
+    .into_iter()
+    .collect();
     let g: Poly<StaticVar> = [(1, [("x", 2)])].into();
     let expected_q: Poly<StaticVar> = [(1, [("x", 1)]), (1, [("y", 1)])].into();
     let (d, q, r) = f.clone().div_rem(&g);
@@ -273,7 +290,9 @@ fn test_div_rem() {
     let f: Poly<StaticVar> = [
         (1i32, Mono::from([("x", 2u32)])),
         (1i32, Mono::from([("x", 1u32), ("y", 1u32)])),
-    ].into_iter().collect();
+    ]
+    .into_iter()
+    .collect();
     let g: Poly<StaticVar> = [(1, [("x", 1)]), (1, [("y", 1)])].into();
     let expected_q: Poly<StaticVar> = [(1, [("x", 1)])].into();
     let (d, q, r) = f.clone().div_rem(&g);
@@ -281,3 +300,80 @@ fn test_div_rem() {
     assert!(r.is_zero());
     verify_div_rem(f, &g, d, q, r);
 }
+
+#[test]
+fn test_groebner() {
+    // Ideal (x²) — already a Gröbner basis
+    let f: Poly<StaticVar> = [(1, [("x", 2)])].into();
+    let _basis = groebner_basis(vec![f]);
+
+    // Ideal (x³ - x², x² - x): gcd is x² - x
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 3u32)])),
+        (-1i32, Mono::from([("x", 2u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let g: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32)])),
+        (-1i32, Mono::from([("x", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let basis = groebner_basis(vec![f, g]);
+    assert!(is_groebner_basis(&basis));
+
+    // Classic example: I = (x²y - x, xy² - y)
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32), ("y", 1u32)])),
+        (-1i32, Mono::from([("x", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let g: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 1u32), ("y", 2u32)])),
+        (-1i32, Mono::from([("y", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let basis = groebner_basis(vec![f, g]);
+    assert!(is_groebner_basis(&basis));
+}
+
+#[test]
+fn test_is_groebner_basis() {
+    // {x} is a GB: only one element, no pairs to check.
+    let f: Poly<StaticVar> = [(1, [("x", 1)])].into();
+    assert!(is_groebner_basis(&[f]));
+
+    // {x, y} is a GB: S(x, y) = y*x - x*y = 0.
+    let x: Poly<StaticVar> = [(1, [("x", 1)])].into();
+    let y: Poly<StaticVar> = [(1, [("y", 1)])].into();
+    assert!(is_groebner_basis(&[x, y]));
+
+    // {x² + y, xy} is NOT a GB:
+    // S(x²+y, xy) = y*(x²+y) - x*(xy) = y² ≠ 0 mod the set.
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32)])),
+        (1i32, Mono::from([("y", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let g: Poly<StaticVar> = [(1i32, Mono::from([("x", 1u32), ("y", 1u32)]))]
+        .into_iter()
+        .collect();
+    assert!(!is_groebner_basis(&[f, g]));
+
+    // After running groebner_basis, the result must always pass.
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32)])),
+        (1i32, Mono::from([("y", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let g: Poly<StaticVar> = [(1i32, Mono::from([("x", 1u32), ("y", 1u32)]))]
+        .into_iter()
+        .collect();
+    let basis = groebner_basis(vec![f, g]);
+    assert!(is_groebner_basis(&basis));
+}