rssn/symbolic/
grobner.rs

1//! # Gröbner Bases
2//!
3//! This module provides functions for computing Gröbner bases of polynomial ideals.
4//! It includes implementations for multivariate polynomial division, S-polynomial computation,
5//! and Buchberger's algorithm for generating a Gröbner basis. Monomial orderings are also supported.
6//!
7//! ## Overview
8//!
9//! A Gröbner basis is a particular kind of generating set of an ideal in a polynomial ring.
10//! Gröbner bases have many applications in computational algebra, including:
11//! - Solving systems of polynomial equations
12//! - Ideal membership testing
13//! - Computing polynomial remainders
14//! - Elimination theory
15//!
16//! ## Examples
17//!
18//! ### Computing a Gröbner Basis
19//! ```
20//! use rssn::symbolic::core::{Expr, Monomial, SparsePolynomial};
21//! use rssn::symbolic::grobner::{buchberger, MonomialOrder};
22//! use std::collections::BTreeMap;
23//!
24//! // Create polynomials: x^2 - y and xy - 1
25//! let mut poly1_terms = BTreeMap::new();
26//! let mut mono1 = BTreeMap::new();
27//! mono1.insert("x".to_string(), 2);
28//! poly1_terms.insert(Monomial(mono1), Expr::new_constant(1.0));
29//! let mut mono2 = BTreeMap::new();
30//! mono2.insert("y".to_string(), 1);
31//! poly1_terms.insert(Monomial(mono2), Expr::new_constant(-1.0));
32//!
33//! let poly1 = SparsePolynomial { terms: poly1_terms };
34//!
35//! // Compute Gröbner basis
36//! let basis = vec![poly1];
37//! let grobner = buchberger(&basis, MonomialOrder::Lexicographical).unwrap();
38//! ```
39
40use crate::symbolic::core::{Expr, Monomial, SparsePolynomial};
41use crate::symbolic::polynomial::{add_poly, mul_poly};
42use crate::symbolic::simplify::is_zero;
43use crate::symbolic::simplify_dag::simplify;
44use std::cmp::Ordering;
45use std::collections::BTreeMap;
46/// Defines the monomial ordering to be used in polynomial division.
47#[derive(Debug, Clone, Copy, serde::Serialize, serde::Deserialize)]
48#[repr(C)]
49pub enum MonomialOrder {
50    Lexicographical,
51    GradedLexicographical,
52    GradedReverseLexicographical,
53}
54/// Compares two monomials based on a given ordering.
55pub(crate) fn compare_monomials(m1: &Monomial, m2: &Monomial, order: MonomialOrder) -> Ordering {
56    match order {
57        MonomialOrder::Lexicographical => m1.0.iter().cmp(m2.0.iter()),
58        _ => m1.0.iter().cmp(m2.0.iter()),
59    }
60}
61/// Performs division of a multivariate polynomial `f` by a set of divisors `F`.
62///
63/// This function implements the multivariate division algorithm, which is a generalization
64/// of polynomial long division. It repeatedly subtracts multiples of the divisors from `f`
65/// until a remainder is obtained that cannot be further reduced.
66///
67/// # Arguments
68/// * `f` - The dividend polynomial.
69/// * `divisors` - A slice of `SparsePolynomial`s representing the divisors.
70/// * `order` - The `MonomialOrder` to use for division.
71///
72/// # Returns
73/// A tuple `(quotients, remainder)` where `quotients` is a `Vec<SparsePolynomial>`
74/// (one for each divisor) and `remainder` is a `SparsePolynomial`.
75pub fn poly_division_multivariate(
76    f: &SparsePolynomial,
77    divisors: &[SparsePolynomial],
78    order: MonomialOrder,
79) -> Result<(Vec<SparsePolynomial>, SparsePolynomial), String> {
80    let mut quotients = vec![
81        SparsePolynomial {
82            terms: BTreeMap::new()
83        };
84        divisors.len()
85    ];
86    let mut remainder = SparsePolynomial {
87        terms: BTreeMap::new(),
88    };
89    let mut p = f.clone();
90    while !p.terms.is_empty() {
91        let mut division_occurred = false;
92        let lead_term_p = match p.terms.keys().max_by(|a, b| compare_monomials(a, b, order)) {
93            Some(lt) => lt.clone(),
94            None => continue,
95        };
96        for (i, divisor) in divisors.iter().enumerate() {
97            if divisor.terms.is_empty() {
98                continue;
99            }
100            let lead_term_g = match divisor
101                .terms
102                .keys()
103                .max_by(|a, b| compare_monomials(a, b, order))
104            {
105                Some(lt) => lt.clone(),
106                None => unreachable!(),
107            };
108            if is_divisible(&lead_term_p, &lead_term_g) {
109                let coeff_p = match p.terms.get(&lead_term_p) {
110                    Some(c) => c,
111                    None => {
112                        return Err("Logic error: lead term not in polynomial terms".to_string());
113                    }
114                };
115                let coeff_g = match divisor.terms.get(&lead_term_g) {
116                    Some(c) => c,
117                    None => {
118                        return Err("Logic error: lead term not found in divisor terms".to_string());
119                    }
120                };
121                let coeff_ratio = simplify(&Expr::new_div(coeff_p.clone(), coeff_g.clone()));
122                let mono_ratio = subtract_monomials(&lead_term_p, &lead_term_g);
123                let mut t_terms = BTreeMap::new();
124                t_terms.insert(mono_ratio, coeff_ratio);
125                let t = SparsePolynomial { terms: t_terms };
126                quotients[i] = add_poly(&quotients[i], &t);
127                let t_g = mul_poly(&t, divisor);
128                p = subtract_poly(&p, &t_g);
129                division_occurred = true;
130                break;
131            }
132        }
133        if !division_occurred {
134            let coeff = match p.terms.remove(&lead_term_p) {
135                Some(c) => c,
136                None => {
137                    return Err("Logic error: lead term not found for removal".to_string());
138                }
139            };
140            remainder.terms.insert(lead_term_p, coeff);
141        }
142    }
143    Ok((quotients, remainder))
144}
145pub(crate) fn is_divisible(m1: &Monomial, m2: &Monomial) -> bool {
146    m2.0.iter()
147        .all(|(var, exp2)| m1.0.get(var).is_some_and(|exp1| exp1 >= exp2))
148}
149pub(crate) fn subtract_monomials(m1: &Monomial, m2: &Monomial) -> Monomial {
150    let mut result = m1.0.clone();
151    for (var, exp2) in &m2.0 {
152        let exp1 = result.entry(var.clone()).or_insert(0);
153        *exp1 -= exp2;
154    }
155    Monomial(result.into_iter().filter(|(_, exp)| *exp > 0).collect())
156}
157pub fn subtract_poly(p1: &SparsePolynomial, p2: &SparsePolynomial) -> SparsePolynomial {
158    let mut result_terms = p1.terms.clone();
159    for (mono, coeff) in &p2.terms {
160        let entry = result_terms
161            .entry(mono.clone())
162            .or_insert_with(|| Expr::Constant(0.0));
163        *entry = simplify(&Expr::new_sub(entry.clone(), coeff.clone()));
164    }
165    result_terms.retain(|_, v| !is_zero(v));
166    SparsePolynomial {
167        terms: result_terms,
168    }
169}
170/// Computes the leading term (monomial, coefficient) of a polynomial.
171pub(crate) fn leading_term(p: &SparsePolynomial, order: MonomialOrder) -> Option<(Monomial, Expr)> {
172    p.terms
173        .iter()
174        .max_by(|(m1, _), (m2, _)| compare_monomials(m1, m2, order))
175        .map(|(m, c)| (m.clone(), c.clone()))
176}
177/// Computes the least common multiple (LCM) of two monomials.
178pub(crate) fn lcm_monomial(m1: &Monomial, m2: &Monomial) -> Monomial {
179    let mut lcm_map = m1.0.clone();
180    for (var, &exp2) in &m2.0 {
181        let exp1 = lcm_map.entry(var.clone()).or_insert(0);
182        *exp1 = std::cmp::max(*exp1, exp2);
183    }
184    Monomial(lcm_map)
185}
186/// Computes the S-polynomial of two polynomials.
187/// S(f, g) = (lcm(LM(f), LM(g)) / LT(f)) * f - (lcm(LM(f), LM(g)) / LT(g)) * g
188pub(crate) fn s_polynomial(
189    p1: &SparsePolynomial,
190    p2: &SparsePolynomial,
191    order: MonomialOrder,
192) -> Option<SparsePolynomial> {
193    let (lm1, lc1) = leading_term(p1, order)?;
194    let (lm2, lc2) = leading_term(p2, order)?;
195    let lcm = lcm_monomial(&lm1, &lm2);
196    let t1_mono = subtract_monomials(&lcm, &lm1);
197    let t1_coeff = simplify(&Expr::new_div(Expr::Constant(1.0), lc1));
198    let mut t1_terms = BTreeMap::new();
199    t1_terms.insert(t1_mono, t1_coeff);
200    let t1 = SparsePolynomial { terms: t1_terms };
201    let t2_mono = subtract_monomials(&lcm, &lm2);
202    let t2_coeff = simplify(&Expr::new_div(Expr::Constant(1.0), lc2));
203    let mut t2_terms = BTreeMap::new();
204    t2_terms.insert(t2_mono, t2_coeff);
205    let t2 = SparsePolynomial { terms: t2_terms };
206    let term1 = mul_poly(&t1, p1);
207    let term2 = mul_poly(&t2, p2);
208    Some(subtract_poly(&term1, &term2))
209}
210/// Computes a Gröbner basis for a polynomial ideal using Buchberger's algorithm.
211///
212/// Buchberger's algorithm takes a set of generators for a polynomial ideal and
213/// produces a Gröbner basis for that ideal. A Gröbner basis has many desirable
214/// properties, such as simplifying polynomial division and solving systems of
215/// polynomial equations.
216///
217/// # Arguments
218/// * `basis` - A slice of `SparsePolynomial`s representing the initial generators of the ideal.
219/// * `order` - The `MonomialOrder` to use for computations.
220///
221/// # Returns
222/// A `Vec<SparsePolynomial>` representing the Gröbner basis.
223pub fn buchberger(
224    basis: &[SparsePolynomial],
225    order: MonomialOrder,
226) -> Result<Vec<SparsePolynomial>, String> {
227    if basis.is_empty() {
228        return Ok(vec![]);
229    }
230    let mut g = basis.to_vec();
231    let mut pairs: Vec<(usize, usize)> = (0..g.len())
232        .flat_map(|i| (i + 1..g.len()).map(move |j| (i, j)))
233        .collect();
234    while let Some((i, j)) = pairs.pop() {
235        if let Some(s_poly) = s_polynomial(&g[i], &g[j], order) {
236            let (_, remainder) = poly_division_multivariate(&s_poly, &g, order)?;
237            if !remainder.terms.is_empty() {
238                let new_poly_idx = g.len();
239                for k in 0..new_poly_idx {
240                    pairs.push((k, new_poly_idx));
241                }
242                g.push(remainder);
243            }
244        }
245    }
246    Ok(g)
247}
rssn/symbolic/grobner.rs

rssn/symbolic/
grobner.rs
