99 lines
2.8 KiB
Rust
99 lines
2.8 KiB
Rust
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
|
|
}
|