check reduced grobner basis

src/poly/buchberger.rs

@@ -0,0 +1,98 @@
+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;
+    }
+
+    // Minimize: remove any g whose LM is divisible by LM(h) for some other h.
+    let mut i = 0;
+    while i < g.len() {
+        let lm_i = g[i].leading_term_lex().unwrap().0.clone();
+        let redundant =
+            (0..g.len()).any(|j| j != i && lm_i.contains(&g[j].leading_term_lex().unwrap().0));
+        if redundant {
+            g.remove(i);
+        } else {
+            i += 1;
+        }
+    }
+
+    // Interreduce: replace each generator with its normal form modulo the others.
+    for i in 0..g.len() {
+        let others: Vec<Poly<V>> = g
+            .iter()
+            .enumerate()
+            .filter(|(j, _)| *j != i)
+            .map(|(_, p)| p.clone())
+            .collect();
+        g[i] = reduce(&g[i], &others);
+    }
+
+    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.
+pub(crate) 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 {
@@ -116,6 +124,34 @@ impl<V: Var> Mono<V> {
         true
     }
 
+    /// Returns the lcm of self and other (element-wise max of exponents).
+    pub fn lcm(&self, other: &Mono<V>) -> Mono<V> {
+        let mut self_it = self.term.iter().peekable();
+        let mut other_it = other.term.iter().peekable();
+        let mut result: Vec<(V, u32)> = vec![];
+
+        loop {
+            match (self_it.peek(), other_it.peek()) {
+                (None, None) => break,
+                (Some(_), None) => result.push(self_it.next().unwrap().clone()),
+                (None, Some(_)) => result.push(other_it.next().unwrap().clone()),
+                (Some((s_var, _)), Some((o_var, _))) => {
+                    if s_var < o_var {
+                        result.push(self_it.next().unwrap().clone());
+                    } else if s_var > o_var {
+                        result.push(other_it.next().unwrap().clone());
+                    } else {
+                        let (var, s_exp) = self_it.next().unwrap();
+                        let (_, o_exp) = other_it.next().unwrap();
+                        result.push((var.clone(), (*s_exp).max(*o_exp)));
+                    }
+                }
+            }
+        }
+
+        Mono { term: result }
+    }
+
     /// Divides self by other. Assumes `self.contains(other)`.
     pub fn div(self, other: &Mono<V>) -> Mono<V> {
         let mut self_it = self.term.into_iter().peekable();

src/poly/fmt.rs

@@ -4,6 +4,7 @@ use itertools::Itertools;
 
 use crate::fmt::{num_to_subscript, num_to_superscript};
 use crate::poly::flat::{Mono, Poly};
+use crate::poly::ideal::Ideal;
 use crate::poly::var::{StaticVar, Var};
 
 impl Display for StaticVar {
@@ -44,6 +45,20 @@ impl<V: Var> Display for Poly<V> {
     }
 }
 
+impl<V: Var, S> Display for Ideal<V, S> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "<")?;
+        let mut iter = self.generators().iter();
+        if let Some(first) = iter.next() {
+            write!(f, "{first}")?;
+            for p in iter {
+                write!(f, ", {p}")?;
+            }
+        }
+        write!(f, ">")
+    }
+}
+
 impl<V: Var> Display for Mono<V> {
     fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
         write!(

src/poly/ideal.rs

@@ -0,0 +1,66 @@
+use std::marker::PhantomData;
+
+use super::buchberger;
+use super::flat::Poly;
+use super::var::Var;
+
+/// Marker: the ideal's generators are arbitrary polynomials.
+pub struct Generators;
+
+/// Marker: the ideal's generators form a Gröbner basis.
+pub struct GroebnerBasis;
+
+pub struct Ideal<V: Var, S> {
+    generators: Vec<Poly<V>>,
+    _state: PhantomData<S>,
+}
+
+impl<V: Var> Ideal<V, Generators> {
+    pub fn new(generators: Vec<Poly<V>>) -> Self {
+        Ideal {
+            generators,
+            _state: PhantomData,
+        }
+    }
+
+    /// Computes a Gröbner basis for this ideal using Buchberger's algorithm.
+    pub fn groebner_basis(self) -> Ideal<V, GroebnerBasis> {
+        Ideal {
+            generators: buchberger::groebner_basis(self.generators),
+            _state: PhantomData,
+        }
+    }
+}
+
+impl<V: Var> Into<Ideal<V, GroebnerBasis>> for Ideal<V, Generators> {
+    fn into(self) -> Ideal<V, GroebnerBasis> {
+        self.groebner_basis()
+    }
+}
+
+impl<V: Var> FromIterator<Poly<V>> for Ideal<V, Generators> {
+    fn from_iter<T: IntoIterator<Item = Poly<V>>>(iter: T) -> Self {
+        Ideal::new(iter.into_iter().collect())
+    }
+}
+
+impl<V: Var, T: IntoIterator<Item = Poly<V>>> From<T> for Ideal<V, Generators> {
+    fn from(iter: T) -> Self {
+        iter.into_iter().collect()
+    }
+}
+
+impl<V: Var, S> Ideal<V, S> {
+    pub fn generators(&self) -> &[Poly<V>] {
+        &self.generators
+    }
+}
+
+impl<V: Var> Ideal<V, GroebnerBasis> {
+    /// Returns `true` if `p` belongs to this ideal.
+    ///
+    /// Reduces `p` modulo the Gröbner basis; membership holds iff the remainder is zero.
+    pub fn contains(&self, p: &Poly<V>) -> bool {
+        buchberger::reduce(p, &self.generators).is_zero()
+    }
+}

src/poly/mod.rs

@@ -1,5 +1,7 @@
+pub mod buchberger;
 pub mod flat;
 pub mod fmt;
+pub mod ideal;
 pub mod ops;
 pub mod var;
 

src/poly/ops.rs

@@ -128,6 +128,38 @@ impl<V: Var> Sub for Poly<V> {
     }
 }
 
+fn gcd(a: i32, b: i32) -> i32 {
+    let (mut a, mut b) = (a.abs(), b.abs());
+    while b != 0 {
+        let t = b;
+        b = a % b;
+        a = t;
+    }
+    a
+}
+
+impl<V: Var> Poly<V> {
+    /// S-polynomial of `self` and `other` under lex order.
+    ///
+    /// Let `L = lcm(LM(f), LM(g))` and `d = gcd(lc(f), lc(g))`. Returns:
+    /// `(lc(g)/d) * (L/LM(f)) * f - (lc(f)/d) * (L/LM(g)) * g`
+    ///
+    /// The leading terms cancel by construction. Panics if either polynomial is zero.
+    pub fn s_poly(&self, other: &Poly<V>) -> Poly<V> {
+        let (lm_f, lc_f) = self.leading_term_lex().expect("f must be nonzero");
+        let (lm_g, lc_g) = other.leading_term_lex().expect("g must be nonzero");
+
+        let l = lm_f.lcm(&lm_g);
+        let t_f: Poly<V> = 1 * l.clone().div(&lm_f);
+        let t_g: Poly<V> = 1 * l.div(&lm_g);
+
+        let d = gcd(lc_f, lc_g);
+        let left = (lc_g / d) * (t_f * self.clone());
+        let right = (lc_f / d) * (t_g * other.clone());
+        left - right
+    }
+}
+
 impl<V: Var> Poly<V> {
     /// Pseudo-division with remainder.
     ///
@@ -136,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,6 @@
-use super::flat::{lex_cmp, Mono, Poly};
+use super::buchberger::{groebner_basis, is_groebner_basis};
+use super::flat::{Mono, Poly, lex_cmp};
+use super::ideal::{Generators, GroebnerBasis, Ideal};
 use super::var::StaticVar;
 
 #[test]
@@ -145,11 +147,80 @@ fn test_mono_div() {
     assert_eq!(a.div(&b), Mono { term: vec![] });
 }
 
-fn make_const_poly(c: i32) -> Poly<StaticVar> {
+#[test]
-    Poly { mono: [(Mono { term: vec![] }, c)].into_iter().collect() }
+fn test_mono_lcm() {
+    // lcm(x², xy) = x²y
+    let a: Mono<StaticVar> = [("x", 2)].into();
+    let b: Mono<StaticVar> = [("x", 1), ("y", 1)].into();
+    assert_eq!(a.lcm(&b), Mono::from([("x", 2), ("y", 1)]));
+
+    // lcm(x, y) = xy  (disjoint)
+    let a: Mono<StaticVar> = [("x", 1)].into();
+    let b: Mono<StaticVar> = [("y", 1)].into();
+    assert_eq!(a.lcm(&b), Mono::from([("x", 1), ("y", 1)]));
+
+    // lcm(x², x²) = x²  (identical)
+    let a: Mono<StaticVar> = [("x", 2)].into();
+    assert_eq!(a.lcm(&a.clone()), a);
 }
 
-fn verify_div_rem(f: Poly<StaticVar>, g: &Poly<StaticVar>, d: u32, q: Poly<StaticVar>, r: Poly<StaticVar>) {
+#[test]
+fn test_s_poly() {
+    // S(x², xy) = 0: S-poly of two monomials always vanishes
+    let f: Poly<StaticVar> = [(1, [("x", 2)])].into();
+    let g: Poly<StaticVar> = [(1, [("x", 1), ("y", 1)])].into();
+    assert!(f.s_poly(&g).is_zero());
+
+    // f = x² + y, g = xy + z  (lex: LM(f)=x², LM(g)=xy)
+    // L = x²y, t_f = y, t_g = x, d = gcd(1,1) = 1
+    // S = y*(x²+y) - x*(xy+z) = x²y + y² - x²y - xz = y² - xz
+    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)])),
+        (1i32, Mono::from([("z", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    let expected: Poly<StaticVar> = [
+        (1i32, Mono::from([("y", 2u32)])),
+        (-1i32, Mono::from([("x", 1u32), ("z", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+    assert_eq!(f.s_poly(&g), expected);
+
+    // f = 2x + y, g = 3x + z  (same LM=x, d=gcd(2,3)=1)
+    // S = (3/1)*f - (2/1)*g = 3(2x+y) - 2(3x+z) = 3y - 2z
+    let f: Poly<StaticVar> = [(2, [("x", 1)]), (1, [("y", 1)])].into();
+    let g: Poly<StaticVar> = [(3, [("x", 1)]), (1, [("z", 1)])].into();
+    let expected: Poly<StaticVar> = [(3, [("y", 1)]), (-2, [("z", 1)])].into();
+    assert_eq!(f.s_poly(&g), expected);
+
+    // f = 4x, g = 6x  (d = gcd(4,6) = 2)
+    // S = (6/2)*x*4x - (4/2)*x*6x = 3*4x - 2*6x = 12x - 12x = 0
+    let f: Poly<StaticVar> = [(4, [("x", 1)])].into();
+    let g: Poly<StaticVar> = [(6, [("x", 1)])].into();
+    assert!(f.s_poly(&g).is_zero());
+}
+
+fn make_const_poly(c: i32) -> Poly<StaticVar> {
+    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>,
+) {
     // lc(g)^d * f == q * g + r
     let (_, lc_g) = g.leading_term_lex().unwrap();
     let lhs = lc_g.pow(d) * f;
@@ -173,7 +244,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);
@@ -218,7 +291,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);
@@ -226,3 +301,164 @@ 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));
+}
+
+#[test]
+fn test_ideal() {
+    // Construction from Vec and iterator
+    let f: Poly<StaticVar> = [(1, [("x", 2)])].into();
+    let g: Poly<StaticVar> = [(1, [("y", 1)])].into();
+    let ideal = Ideal::new(vec![f.clone(), g.clone()]);
+    assert_eq!(ideal.generators().len(), 2);
+    let ideal: Ideal<StaticVar, Generators> = [f.clone(), g.clone()].into_iter().collect();
+    assert_eq!(ideal.generators().len(), 2);
+
+    // Construction via From / .into()
+    let ideal: Ideal<StaticVar, Generators> = vec![f.clone(), g.clone()].into();
+    assert_eq!(ideal.generators().len(), 2);
+
+    // Display: <x², y>
+    let ideal = Ideal::new(vec![f, g]);
+    assert_eq!(ideal.to_string(), "<x\u{00B2}, y>");
+
+    // groebner_basis transitions state and result satisfies the criterion
+    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 ideal: Ideal<StaticVar, Generators> = [f, g].into_iter().collect();
+    let gb: Ideal<StaticVar, GroebnerBasis> = ideal.groebner_basis();
+    assert!(is_groebner_basis(gb.generators()));
+}
+
+#[test]
+fn test_groebner_sagemath() {
+    // I = (x³ - 2xy, x²y - 2y² + x) ⊆ k[x, y]
+    // grobner basis: {4y³, x - 2y²}
+
+    // f1 = x³ - 2xy
+    let f1: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 3u32)])),
+        (-2i32, Mono::from([("x", 1u32), ("y", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+
+    // f2 = x²y - 2y² + x
+    let f2: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32), ("y", 1u32)])),
+        (-2i32, Mono::from([("y", 2u32)])),
+        (1i32, Mono::from([("x", 1u32)])),
+    ]
+    .into_iter()
+    .collect();
+
+    let gb = Ideal::new(vec![f1.clone(), f2.clone()]).groebner_basis();
+
+    assert!(is_groebner_basis(gb.generators()));
+    assert_eq!(gb.generators().len(), 2);
+
+    // -x + 2y²
+    let neg_x_plus_2y2: Poly<StaticVar> = [
+        (-1i32, Mono::from([("x", 1u32)])),
+        (2i32, Mono::from([("y", 2u32)])),
+    ]
+    .into_iter()
+    .collect();
+
+    // -4y³
+    let neg_4y3: Poly<StaticVar> = [(-4i32, Mono::from([("y", 3u32)]))].into_iter().collect();
+
+    let expected = [neg_x_plus_2y2, neg_4y3];
+    for e in &expected {
+        assert!(
+            gb.generators()
+                .iter()
+                .any(|a| *a == *e || *a == -1i32 * e.clone())
+        );
+    }
+}