diff --git a/src/poly/buchberger.rs b/src/poly/buchberger.rs new file mode 100644 index 0000000..b4fc073 --- /dev/null +++ b/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(generators: Vec>) -> Vec> { + let mut g: Vec> = 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(basis: &[Poly]) -> 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(f: &Poly, basis: &[Poly]) -> Poly { + 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 +} diff --git a/src/poly/flat.rs b/src/poly/flat.rs index fc2bd12..48516ed 100644 --- a/src/poly/flat.rs +++ b/src/poly/flat.rs @@ -11,10 +11,18 @@ pub fn lex_cmp(a: &Mono, b: &Mono) -> 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 { diff --git a/src/poly/mod.rs b/src/poly/mod.rs index 51fd3b1..94f9647 100644 --- a/src/poly/mod.rs +++ b/src/poly/mod.rs @@ -1,3 +1,4 @@ +pub mod buchberger; pub mod flat; pub mod fmt; pub mod ops; diff --git a/src/poly/ops.rs b/src/poly/ops.rs index d07a57e..1bb186b 100644 --- a/src/poly/ops.rs +++ b/src/poly/ops.rs @@ -168,9 +168,7 @@ impl Poly { /// /// Panics if `divisor` is zero. pub fn div_rem(self, divisor: &Poly) -> (u32, Poly, Poly) { - let (lt_g_mono, lt_g_coeff) = divisor - .leading_term_lex() - .expect("divisor must be nonzero"); + let (lt_g_mono, lt_g_coeff) = divisor.leading_term_lex().expect("divisor must be nonzero"); let mut p = self; let mut q = Poly::default(); diff --git a/src/poly/tests.rs b/src/poly/tests.rs index f13836d..7a9536d 100644 --- a/src/poly/tests.rs +++ b/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 = [ (1i32, Mono::from([("x", 2u32)])), (1i32, Mono::from([("y", 1u32)])), - ].into_iter().collect(); + ] + .into_iter() + .collect(); let g: Poly = [ (1i32, Mono::from([("x", 1u32), ("y", 1u32)])), (1i32, Mono::from([("z", 1u32)])), - ].into_iter().collect(); + ] + .into_iter() + .collect(); let expected: Poly = [ (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 { - Poly { mono: [(Mono { term: vec![] }, c)].into_iter().collect() } + Poly { + mono: [(Mono { term: vec![] }, c)].into_iter().collect(), + } } -fn verify_div_rem(f: Poly, g: &Poly, d: u32, q: Poly, r: Poly) { +fn verify_div_rem( + f: Poly, + g: &Poly, + d: u32, + q: Poly, + r: Poly, +) { // 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 = [ (1i32, Mono::from([("x", 3u32)])), (1i32, Mono::from([("x", 2u32), ("y", 1u32)])), - ].into_iter().collect(); + ] + .into_iter() + .collect(); let g: Poly = [(1, [("x", 2)])].into(); let expected_q: Poly = [(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 = [ (1i32, Mono::from([("x", 2u32)])), (1i32, Mono::from([("x", 1u32), ("y", 1u32)])), - ].into_iter().collect(); + ] + .into_iter() + .collect(); let g: Poly = [(1, [("x", 1)]), (1, [("y", 1)])].into(); let expected_q: Poly = [(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 = [(1, [("x", 2)])].into(); + let _basis = groebner_basis(vec![f]); + + // Ideal (x³ - x², x² - x): gcd is x² - x + let f: Poly = [ + (1i32, Mono::from([("x", 3u32)])), + (-1i32, Mono::from([("x", 2u32)])), + ] + .into_iter() + .collect(); + let g: Poly = [ + (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 = [ + (1i32, Mono::from([("x", 2u32), ("y", 1u32)])), + (-1i32, Mono::from([("x", 1u32)])), + ] + .into_iter() + .collect(); + let g: Poly = [ + (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 = [(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 = [(1, [("x", 1)])].into(); + let y: Poly = [(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 = [ + (1i32, Mono::from([("x", 2u32)])), + (1i32, Mono::from([("y", 1u32)])), + ] + .into_iter() + .collect(); + let g: Poly = [(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 = [ + (1i32, Mono::from([("x", 2u32)])), + (1i32, Mono::from([("y", 1u32)])), + ] + .into_iter() + .collect(); + let g: Poly = [(1i32, Mono::from([("x", 1u32), ("y", 1u32)]))] + .into_iter() + .collect(); + let basis = groebner_basis(vec![f, g]); + assert!(is_groebner_basis(&basis)); +}