improve var creation for SumProdCircuit

README.md

@@ -0,0 +1,134 @@
+# circuit-cas
+
+A computer algebra system in Rust for algebraic circuit verification over polynomial rings with integer coefficients.
+
+## Overview
+
+`circuit-cas` provides two main layers:
+
+- **`poly`** — multivariate polynomials over ℤ, Gröbner bases, and ideals
+- **`circuit`** — arithmetic circuits (DAGs of add/mul nodes) quotiented by an ideal
+
+The intended use case is verifying arithmetic circuits modulo an ideal: two circuits computing the same polynomial in the quotient ring are equivalent.
+
+## Modules
+
+### `poly`
+
+#### Variables — `poly::var`
+
+Variables are any type implementing the `Var` trait. The built-in `StaticVar` carries a static string name and optional integer indices, with a `var!` macro for construction:
+
+```rust
+use circuit_cas::var;
+
+let x   = var!("x");          // x
+let x0  = var!("x", 0);       // x₀
+let x01 = var!("x", 0, 1);    // x₀,₁
+```
+
+#### Polynomials — `poly::flat`
+
+`Poly<V>` is a multivariate polynomial over ℤ backed by a `HashMap<Mono<V>, i32>`. Arithmetic operators (`+`, `-`, `*`) and scalar multiplication are provided. Monomial ordering uses pure lexicographic order (lex) with variable names sorted alphabetically.
+
+```rust
+// x² - 2xy  (using array-of-pairs syntax)
+let f: Poly<StaticVar> = [
+    (1i32,  Mono::from([("x", 2u32)])),
+    (-2i32, Mono::from([("x", 1u32), ("y", 1u32)])),
+].into_iter().collect();
+```
+
+Key operations:
+
+| Method | Description |
+|--------|-------------|
+| `Poly::is_zero` | whether the polynomial is the zero polynomial |
+| `Poly::leading_term_lex` | leading `(monomial, coefficient)` under lex |
+| `Poly::s_poly(&other)` | S-polynomial of two polynomials |
+| `Poly::div_rem(divisor)` | pseudo-division: returns `(d, q, r)` with `lc(g)^d·f = q·g + r` |
+
+#### Ideals — `poly::ideal`
+
+`Ideal<V, S>` uses the typestate pattern to track at the type level whether its generators form a Gröbner basis:
+
+```rust
+use circuit_cas::poly::ideal::{Ideal, Generators, GroebnerBasis};
+
+// Arbitrary generators
+let ideal: Ideal<StaticVar, Generators> = Ideal::new(vec![f1, f2]);
+
+// Compute the reduced Gröbner basis — consumes the original ideal
+let gb: Ideal<StaticVar, GroebnerBasis> = ideal.groebner_basis();
+
+// Ideal membership test (only available on GroebnerBasis state)
+assert!(gb.contains(&some_poly));
+```
+
+`Ideal` implements `Display` as `<g1, g2, ...>` and `FromIterator<Poly<V>>` for construction from an iterator.
+
+#### Buchberger's algorithm — `poly::buchberger`
+
+```rust
+use circuit_cas::poly::buchberger::{groebner_basis, is_groebner_basis};
+
+let basis = groebner_basis(vec![f1, f2]);
+assert!(is_groebner_basis(&basis));
+```
+
+`groebner_basis` returns the **reduced** Gröbner basis: after running Buchberger's algorithm it minimizes (removes redundant generators) and interreduces (fully reduces each generator modulo the others).
+
+### `circuit`
+
+#### DAG — `circuit::dag`
+
+`Circuit<V>` is a DAG of `Node<V>` values, with structural sharing via interning. Nodes are `Leaf(V)`, `Scale(id, i32)`, `Sum(id, id)`, or `Prod(id, id)`.
+
+#### Quotient ring — `circuit::quotient`
+
+`Quotient<V>` pairs an arithmetic circuit with a `Ideal<V, GroebnerBasis>`, representing a circuit in the quotient ring `ℤ[vars] / I`.
+
+```rust
+use circuit_cas::circuit::quotient::Quotient;
+
+// Build from a list of generators — Gröbner basis is computed automatically
+let quotient: Quotient<StaticVar> = vec![f1, f2, f3].into_iter().collect();
+
+// Or from a pre-computed Gröbner basis
+let quotient = Quotient::from(gb);
+
+println!("{quotient}"); // C/<g1, g2>
+```
+
+## Example
+
+```rust
+use circuit_cas::{var, circuit::quotient::Quotient, poly::var::StaticVar};
+
+fn main() {
+    let x  = var!("x");
+    let xb = var!("x\u{0304}");  // x̄
+
+    // Idempotency relations: x² = x, x̄² = x̄, x·x̄ = x
+    let idem = vec![
+        1 * (&x  ^ 2) - 1 * (&x  ^ 1),
+        1 * (&xb ^ 2) - 1 * (&xb ^ 1),
+        1 * ((&x ^ 1) * (&xb ^ 1)) - 1 * (&x ^ 1),
+    ];
+
+    let quotient: Quotient<StaticVar> = idem.into_iter().collect();
+    println!("{quotient}");
+}
+```
+
+Run with:
+
+```
+cargo run --example quotient
+```
+
+## Development
+
+```
+cargo test     # run all tests
+```

examples/dag.rs

@@ -1,26 +1,25 @@
 use std::cell::RefCell;
 use std::rc::Rc;
-
+use circuit_cas::circuit::probabilistic::ProbCircuit;
-use circuit_cas::circuit::dag::{Circuit, CircuitExt};
+use circuit_cas::circuit::dag::CircuitExt;
-use circuit_cas::var;
 
 fn main() {
-    let circuit = Rc::new(RefCell::new(Circuit::new()));
+    let circuit: Rc<RefCell<ProbCircuit>> = ProbCircuit::new();
 
-    // Build (x + y) * (x + z)
+    // vars accept anything that implements Into<StaticVar>: &'static str, (&str, u32), (&str, u32, u32)
-    let x = circuit.leaf(var!("x"));
+    let x   = circuit.var("x");
-    let y = circuit.leaf(var!("y"));
+    let y   = circuit.var(("y", 1));       // indexed variable y_1
-    let z = circuit.leaf(var!("z"));
+    let z   = circuit.var(("z", 0, 1));    // doubly-indexed variable z_{0,1}
 
-    let x_plus_y = circuit.leaf(var!("x")) + y;
+    let x_plus_y = circuit.var("x") + y;
-    let x_plus_z = circuit.leaf(var!("x")) + z;
+    let x_plus_z = circuit.var("x") + z;
     let expr = x_plus_y * x_plus_z;
 
     // Deduplication: both x leaves share the same NodeId
-    let x2 = circuit.leaf(var!("x"));
+    let x2 = circuit.var("x");
     assert_eq!(x.id, x2.id);
 
-    println!("(x + y) * (x + z) root node id: {:?}", expr.id);
+    println!("(x + y_1) * (x + z_{{0,1}}) root node id: {:?}", expr.id);
     println!("x node id: {:?}", x.id);
     println!("x deduplicated node id: {:?}", x2.id);
 }

examples/quotient.rs

@@ -0,0 +1,31 @@
+use std::cell::RefCell;
+use std::rc::Rc;
+use circuit_cas::circuit::dag::CircuitExt;
+use circuit_cas::circuit::quotient::QuotientCircuit;
+use circuit_cas::circuit::traits::Circuit;
+use circuit_cas::poly::ideal::{Generators, Ideal};
+use circuit_cas::var;
+
+fn main() {
+    let x = var!("x");
+    let nx = var!("x\u{0304}");
+
+    let idem: Ideal<_, Generators> = vec![
+        1 * (&x ^ 2) - 1 * (&x ^ 1),
+        1 * (&nx ^ 2) - 1 * (&nx ^ 1),
+        1 * ((&x ^ 1) * (&nx ^ 1)) - 1 * (&x ^ 1),
+    ].into();
+
+    let quotient: Rc<RefCell<QuotientCircuit>> = idem.into();
+
+    // Build x * x̄ + x in the DAG
+    // var accepts anything that implements Into<StaticVar>: &'static str, (&str, u32), (&str, u32, u32)
+    let xn   = quotient.var("x");
+    let nxn  = quotient.var("x\u{0304}");
+    let prod = xn * nxn;
+    let xn2  = quotient.var("x");
+    let expr = prod + xn2;
+
+    println!("dag size: {}", quotient.borrow().len());
+    println!("expr node id: {:?}", expr.id);
+}

src/circuit/dag.rs

@@ -1,45 +1,29 @@
 use slotmap::{SlotMap, new_key_type};
 use std::cell::RefCell;
 use std::collections::HashMap;
-use std::ops::{Add, Mul};
 use std::rc::Rc;
 
-use crate::poly::var::Var;
+use super::traits::{Circuit, Node};
 
 new_key_type! { pub struct NodeId; }
 
-#[derive(Clone, PartialEq, Eq, Hash)]
+#[derive(Clone, Debug)]
-pub enum Node<V: Var> {
+pub struct Dag<N: Node> {
-    Leaf(V),
+    nodes: SlotMap<NodeId, N>,
-    Sum(NodeId, NodeId),
+    intern: HashMap<N, NodeId>,
-    Prod(NodeId, NodeId),
 }
 
-impl<V: Var> Node<V> {
+impl<N: Node> Default for Dag<N> {
-    pub fn children(&self) -> impl Iterator<Item = NodeId> {
+    fn default() -> Self {
-        match self {
+        Dag {
-            Node::Leaf(_) => [None, None],
+            nodes: Default::default(),
-            Node::Sum(l, r) | Node::Prod(l, r) => [Some(*l), Some(*r)],
+            intern: Default::default(),
         }
-        .into_iter()
-        .flatten()
     }
 }
 
-pub struct Circuit<V: Var> {
+impl<N: Node> Dag<N> {
-    nodes: SlotMap<NodeId, Node<V>>,
+    pub(super) fn node(&mut self, n: N) -> NodeId {
-    intern: HashMap<Node<V>, NodeId>,
-}
-
-impl<V: Var> Circuit<V> {
-    pub fn new() -> Self {
-        Circuit {
-            nodes: SlotMap::with_key(),
-            intern: HashMap::new(),
-        }
-    }
-
-    pub fn node(&mut self, n: Node<V>) -> NodeId {
         if let Some(&id) = self.intern.get(&n) {
             return id;
         }
@@ -48,90 +32,33 @@ impl<V: Var> Circuit<V> {
         id
     }
 
-    pub fn get(&self, id: NodeId) -> Option<&Node<V>> {
+    pub(super) fn get(&self, id: NodeId) -> Option<&N> {
         self.nodes.get(id)
     }
 
-    pub fn add(&mut self, left: NodeId, right: NodeId) -> NodeId {
+    pub(super) fn len(&self) -> usize {
-        let (l, r) = if left <= right {
-            (left, right)
-        } else {
-            (right, left)
-        };
-        self.node(Node::Sum(l, r))
-    }
-
-    pub fn mul(&mut self, left: NodeId, right: NodeId) -> NodeId {
-        let (l, r) = if left <= right {
-            (left, right)
-        } else {
-            (right, left)
-        };
-        self.node(Node::Prod(l, r))
-    }
-
-    pub fn len(&self) -> usize {
         self.nodes.len()
     }
 
-    pub fn children(&self, id: NodeId) -> impl Iterator<Item = NodeId> + '_ {
+    pub(super) fn children(&self, id: NodeId) -> impl Iterator<Item = NodeId> + '_ {
         self.nodes.get(id).into_iter().flat_map(Node::children)
     }
 
-    pub fn remove(&mut self, id: NodeId) {
+    pub(super) fn remove(&mut self, id: NodeId) {
         if let Some(node) = self.nodes.remove(id) {
             self.intern.remove(&node);
         }
     }
 }
 
-pub struct CircuitNode<V: Var> {
+pub struct RefNode<C: Circuit> {
     pub id: NodeId,
-    circuit: Rc<RefCell<Circuit<V>>>,
+    pub(super) circuit: Rc<RefCell<C>>,
-}
-pub trait CircuitExt<V: Var> {
-    fn leaf(&self, v: V) -> CircuitNode<V>;
-    fn get_node(&self, id: NodeId) -> Option<CircuitNode<V>>;
 }
 
-impl<V: Var> CircuitExt<V> for Rc<RefCell<Circuit<V>>> {
-    fn leaf(&self, v: V) -> CircuitNode<V> {
-        let id = self.borrow_mut().node(Node::Leaf(v));
-        CircuitNode {
-            id,
-            circuit: self.clone(),
-        }
-    }
 
-    fn get_node(&self, id: NodeId) -> Option<CircuitNode<V>> {
+pub trait CircuitExt {
-        self.borrow().get(id)?;
+    type C: Circuit;
-        Some(CircuitNode {
+    type Var;
-            id,
+    fn var<T: Into<Self::Var>>(&self, v: T) -> RefNode<Self::C>;
-            circuit: self.clone(),
-        })
-    }
-}
-
-impl<V: Var> Add for CircuitNode<V> {
-    type Output = Self;
-
-    fn add(self, rhs: Self) -> Self {
-        let id = self.circuit.borrow_mut().add(self.id, rhs.id);
-        CircuitNode {
-            id,
-            circuit: self.circuit,
-        }
-    }
-}
-
-impl<V: Var> Mul for CircuitNode<V> {
-    type Output = Self;
-
-    fn mul(self, rhs: Self) -> Self {
-        let id = self.circuit.borrow_mut().mul(self.id, rhs.id);
-        CircuitNode {
-            id,
-            circuit: self.circuit,
-        }
-    }
 }

src/circuit/mod.rs

@@ -1,6 +1,7 @@
+pub mod traits;
 pub mod dag;
+pub mod probabilistic;
 pub mod quotient;
 
 #[cfg(test)]
 mod tests;
-

src/circuit/probabilistic.rs

@@ -0,0 +1,68 @@
+use std::cell::RefCell;
+use std::ops::{Deref, DerefMut};
+use std::rc::Rc;
+
+use crate::poly::var::{StaticVar, Var};
+use super::dag::{Dag, NodeId};
+use super::traits::{Node, SumProdCircuit};
+
+
+#[derive(Clone, Debug, PartialEq, Eq, Hash)]
+pub enum PNode<V: Var> {
+    Var(V),
+    Sum(NodeId, NodeId),
+    Prod(NodeId, NodeId),
+}
+
+impl<V: Var> Node for PNode<V> {
+    fn children(&self) -> impl Iterator<Item = NodeId> {
+        match self {
+            Self::Var(_) => [None, None],
+            Self::Sum(l, r) | Self::Prod(l, r) => [Some(*l), Some(*r)],
+        }
+        .into_iter()
+        .flatten()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct ProbCircuit<V: Var = StaticVar> {
+    dag: Dag<PNode<V>>,
+}
+
+impl<V: Var> Default for ProbCircuit<V> {
+    fn default() -> Self { ProbCircuit { dag: Dag::default() } }
+}
+
+impl<V: Var> ProbCircuit<V> {
+    pub fn new() -> Rc<RefCell<Self>> {
+        Rc::new(RefCell::new(Self::default()))
+    }
+}
+
+impl<V: Var> SumProdCircuit for ProbCircuit<V> {
+    type Var = V;
+
+    fn var<T: Into<V>>(&mut self, v: T) -> NodeId { self.node(PNode::Var(v.into())) }
+
+    fn add(&mut self, l: NodeId, r: NodeId) -> NodeId {
+        let (l, r) = if l <= r { (l, r) } else { (r, l) };
+        self.node(PNode::Sum(l, r))
+    }
+
+    fn mul(&mut self, l: NodeId, r: NodeId) -> NodeId {
+        let (l, r) = if l <= r { (l, r) } else { (r, l) };
+        self.node(PNode::Prod(l, r))
+    }
+}
+
+impl<V: Var> Deref for ProbCircuit<V> {
+    type Target = Dag<PNode<V>>;
+    fn deref(&self) -> &Self::Target { &self.dag }
+}
+
+impl<V: Var> DerefMut for ProbCircuit<V> {
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.dag }
+}
+
+

src/circuit/quotient.rs

@@ -1 +1,73 @@
+use std::cell::RefCell;
+use std::ops::{Deref, DerefMut};
+use std::rc::Rc;
+
+use crate::poly::var::{StaticVar, Var};
+use crate::poly::ideal::{Generators, GroebnerBasis, Ideal};
+use super::dag::{Dag, NodeId};
+use super::traits::{Node, SumProdCircuit};
+
+
+#[derive(Clone, Debug, PartialEq, Eq, Hash)]
+pub enum QNode<V: Var> {
+    Var(V),
+    Sum(NodeId, NodeId),
+    Prod(NodeId, NodeId),
+    DivStep(NodeId, NodeId),
+}
+
+impl<V: Var> Node for QNode<V> {
+    fn children(&self) -> impl Iterator<Item = NodeId> {
+        match self {
+            Self::Var(_) => [None, None],
+            Self::Sum(l, r) | Self::Prod(l, r) | Self::DivStep(l, r) => [Some(*l), Some(*r)],
+        }
+        .into_iter()
+        .flatten()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct QuotientCircuit<V: Var = StaticVar> {
+    basis: Ideal<V, GroebnerBasis>,
+    dag: Dag<QNode<V>>,
+}
+
+impl<V: Var> From<Ideal<V, Generators>> for Rc<RefCell<QuotientCircuit<V>>> {
+    fn from(ideal: Ideal<V, Generators>) -> Self {
+        ideal.groebner_basis().into()
+    }
+}
+
+impl<V: Var> From<Ideal<V, GroebnerBasis>> for Rc<RefCell<QuotientCircuit<V>>> {
+    fn from(basis: Ideal<V, GroebnerBasis>) -> Self {
+        Rc::new(RefCell::new(QuotientCircuit { basis, dag: Default::default() }))
+    }
+}
+
+impl<V: Var> Deref for QuotientCircuit<V> {
+    type Target = Dag<QNode<V>>;
+    fn deref(&self) -> &Self::Target { &self.dag }
+}
+
+impl<V: Var> DerefMut for QuotientCircuit<V> {
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.dag }
+}
+
+impl<V: Var> SumProdCircuit for QuotientCircuit<V> {
+    type Var = V;
+
+    fn var<T: Into<V>>(&mut self, v: T) -> NodeId { self.node(QNode::Var(v.into())) }
+
+    fn add(&mut self, l: NodeId, r: NodeId) -> NodeId {
+        let (l, r) = if l <= r { (l, r) } else { (r, l) };
+        self.node(QNode::Sum(l, r))
+    }
+
+    fn mul(&mut self, l: NodeId, r: NodeId) -> NodeId {
+        let (l, r) = if l <= r { (l, r) } else { (r, l) };
+        self.node(QNode::Prod(l, r))
+    }
+}
+
 

src/circuit/tests.rs

@@ -1,51 +1,48 @@
-
 use std::cell::RefCell;
 use std::rc::Rc;
-
+use super::probabilistic::ProbCircuit;
-use super::dag::{Circuit, CircuitExt};
+use super::dag::CircuitExt;
-use crate::poly::var::StaticVar;
-
 #[test]
 fn test_deduplication() {
-    let circuit = Rc::new(RefCell::new(Circuit::new()));
+    let circuit: Rc<RefCell<ProbCircuit>> = ProbCircuit::new();
 
     // Same leaf constructed twice returns the same NodeId
-    let x1 = circuit.leaf(StaticVar::from("x"));
+    let x1 = circuit.var("x");
-    let x2 = circuit.leaf(StaticVar::from("x"));
+    let x2 = circuit.var("x");
     assert_eq!(x1.id, x2.id);
     assert_eq!(circuit.borrow().len(), 1);
 
     // Same sum constructed twice returns the same NodeId
-    let _y = circuit.leaf(StaticVar::from("y"));
+    let _y = circuit.var("y");
-    let sum1 = circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y"));
+    let sum1 = circuit.var("x") + circuit.var("y");
-    let sum2 = circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y"));
+    let sum2 = circuit.var("x") + circuit.var("y");
     assert_eq!(sum1.id, sum2.id);
     assert_eq!(circuit.borrow().len(), 3); // x, y, x+y
 
     // Shared subexpression: (x + y) * (x + y) reuses the x+y node
-    let xy = circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y"));
+    let xy = circuit.var("x") + circuit.var("y");
-    let xy2 = circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y"));
+    let xy2 = circuit.var("x") + circuit.var("y");
     let _sq = xy * xy2;
     assert_eq!(circuit.borrow().len(), 4); // x, y, x+y, (x+y)*(x+y)
 
     // Commutativity: x+y and y+x are the same node
-    let xy = circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y"));
+    let xy = circuit.var("x") + circuit.var("y");
-    let yx = circuit.leaf(StaticVar::from("y")) + circuit.leaf(StaticVar::from("x"));
+    let yx = circuit.var("y") + circuit.var("x");
     assert_eq!(xy.id, yx.id);
 
     // Associativity: (x+y)+z and x+(y+z) are distinct nodes
-    let _z = circuit.leaf(StaticVar::from("z"));
+    let _z = circuit.var("z");
-    let xy_z = (circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y")))
+    let xy_z = (circuit.var("x") + circuit.var("y"))
-        + circuit.leaf(StaticVar::from("z"));
+        + circuit.var("z");
-    let x_yz = circuit.leaf(StaticVar::from("x"))
+    let x_yz = circuit.var("x")
-        + (circuit.leaf(StaticVar::from("y")) + circuit.leaf(StaticVar::from("z")));
+        + (circuit.var("y") + circuit.var("z"));
     assert_ne!(xy_z.id, x_yz.id);
 
     // Deep shared structure: (x+y)*z appears twice in ((x+y)*z) + ((x+y)*z)
-    let xyz1 = (circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y")))
+    let xyz1 = (circuit.var("x") + circuit.var("y"))
-        * circuit.leaf(StaticVar::from("z"));
+        * circuit.var("z");
-    let xyz2 = (circuit.leaf(StaticVar::from("x")) + circuit.leaf(StaticVar::from("y")))
+    let xyz2 = (circuit.var("x") + circuit.var("y"))
-        * circuit.leaf(StaticVar::from("z"));
+        * circuit.var("z");
     assert_eq!(xyz1.id, xyz2.id);
     let _sum = xyz1 + xyz2;
     // x, y, z, x+y(==y+x), (x+y)*z, (x+y)+z, y+z, x+(y+z), (x+y)*z+(x+y)*z, sq

src/circuit/traits.rs

@@ -0,0 +1,63 @@
+use std::cell::RefCell;
+use std::hash::Hash;
+use std::ops::{Add, DerefMut, Mul};
+use std::rc::Rc;
+
+use super::dag::{CircuitExt, Dag, NodeId, RefNode};
+
+pub trait Node: Clone + PartialEq + Eq + Hash {
+    fn children(&self) -> impl Iterator<Item = NodeId>;
+}
+
+pub trait Circuit: Clone + DerefMut<Target = Dag<Self::Node>> {
+    type Node: Node;
+
+    fn node(&mut self, n: Self::Node) -> NodeId { (**self).node(n) }
+    fn remove(&mut self, id: NodeId) { (**self).remove(id) }
+    fn get(&self, id: NodeId) -> Option<&Self::Node> { (**self).get(id) }
+    fn len(&self) -> usize { (**self).len() }
+    fn children(&self, id: NodeId) -> impl Iterator<Item = NodeId> + '_ { (**self).children(id) }
+}
+
+impl<T, N> Circuit for T
+where
+    T: Clone + DerefMut<Target = Dag<N>>,
+    N: Node,
+{
+    type Node = N;
+}
+
+pub trait SumProdCircuit: Circuit {
+    type Var;
+    fn var<T: Into<Self::Var>>(&mut self, v: T) -> NodeId;
+    fn add(&mut self, l: NodeId, r: NodeId) -> NodeId;
+    fn mul(&mut self, l: NodeId, r: NodeId) -> NodeId;
+}
+
+impl<C: SumProdCircuit> Add for RefNode<C> {
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self {
+        let id = self.circuit.borrow_mut().add(self.id, rhs.id);
+        RefNode { id, circuit: self.circuit }
+    }
+}
+
+impl<C: SumProdCircuit> Mul for RefNode<C> {
+    type Output = Self;
+
+    fn mul(self, rhs: Self) -> Self {
+        let id = self.circuit.borrow_mut().mul(self.id, rhs.id);
+        RefNode { id, circuit: self.circuit }
+    }
+}
+
+impl<C: SumProdCircuit> CircuitExt for Rc<RefCell<C>> {
+    type C = C;
+    type Var = C::Var;
+
+    fn var<T: Into<C::Var>>(&self, v: T) -> RefNode<C> {
+        let id = self.borrow_mut().var(v);
+        RefNode { id, circuit: self.clone() }
+    }
+}

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

@@ -1,7 +1,54 @@
+use std::cmp::Ordering;
 use std::collections::HashMap;
 
 use super::var::Var;
 
+pub fn lex_cmp<V: Var>(a: &Mono<V>, b: &Mono<V>) -> Ordering {
+    let mut a_it = a.vars.iter().peekable();
+    let mut b_it = b.vars.iter().peekable();
+
+    loop {
+        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
+                };
+            }
+            (None, Some((_, b_exp))) => {
+                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 {
+                    if *a_exp > 0 {
+                        return Ordering::Greater;
+                    }
+                    a_it.next();
+                } else if a_var > b_var {
+                    if *b_exp > 0 {
+                        return Ordering::Less;
+                    }
+                    b_it.next();
+                } else {
+                    match a_exp.cmp(b_exp) {
+                        Ordering::Equal => {
+                            a_it.next();
+                            b_it.next();
+                        }
+                        ord => return ord,
+                    }
+                }
+            }
+        }
+    }
+}
+
 #[derive(Clone, Debug, PartialEq)]
 pub struct Poly<V: Var> {
     pub(crate) mono: HashMap<Mono<V>, i32>,
@@ -15,6 +62,20 @@ impl<V: Var> Default for Poly<V> {
     }
 }
 
+impl<V: Var> Poly<V> {
+    pub fn is_zero(&self) -> bool {
+        self.mono.is_empty()
+    }
+
+    /// Returns (leading monomial, leading coefficient) under lex order, or None if zero.
+    pub fn leading_term_lex(&self) -> Option<(Mono<V>, i32)> {
+        self.mono
+            .iter()
+            .max_by(|(m1, _), (m2, _)| lex_cmp(m1, m2))
+            .map(|(m, &c)| (m.clone(), c))
+    }
+}
+
 impl<V: Var, T: IntoIterator> From<T> for Poly<V>
 where
     Poly<V>: FromIterator<<T as IntoIterator>::Item>,
@@ -34,20 +95,20 @@ impl<V: Var, U: Into<Mono<V>>> FromIterator<(i32, U)> for Poly<V> {
 
 #[derive(Clone, Debug, Default, PartialEq, Eq, Hash, PartialOrd, Ord)]
 pub struct Mono<V: Var> {
-    pub term: Vec<(V, u32)>,
+    pub vars: Vec<(V, u32)>,
 }
 
 impl<V: Var> Mono<V> {
     pub fn contains(&self, other: &Mono<V>) -> bool {
-        let mut self_it = self.term.iter().peekable();
+        let mut self_it = self.vars.iter().peekable();
-        let mut other_it = other.term.iter().peekable();
+        let mut other_it = other.vars.iter().peekable();
 
-        while let Some((o_term, o_exp)) = other_it.peek() {
+        while let Some((o_var, o_exp)) = other_it.peek() {
-            if let Some((s_term, s_exp)) = self_it.peek() {
+            if let Some((s_var, s_exp)) = self_it.peek() {
-                if s_term < o_term {
+                if s_var < o_var {
                     self_it.next();
                     continue;
-                } else if s_term > o_term {
+                } else if s_var > o_var {
                     return false;
                 } else if o_exp <= s_exp {
                     self_it.next();
@@ -62,6 +123,64 @@ 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.vars.iter().peekable();
+        let mut other_it = other.vars.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 { vars: result }
+    }
+
+    /// Divides self by other. Assumes `self.contains(other)`.
+    pub fn div(self, other: &Mono<V>) -> Mono<V> {
+        let mut self_it = self.vars.into_iter().peekable();
+        let mut other_it = other.vars.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()),
+                (None, Some(_)) => unreachable!("divisor not contained in dividend"),
+                (Some((s_var, _)), Some((o_var, _))) => {
+                    if s_var < o_var {
+                        result.push(self_it.next().unwrap());
+                    } else if s_var > o_var {
+                        unreachable!("divisor not contained in dividend");
+                    } else {
+                        let (var, s_exp) = self_it.next().unwrap();
+                        let (_, o_exp) = other_it.next().unwrap();
+                        if s_exp > *o_exp {
+                            result.push((var, s_exp - o_exp));
+                        }
+                    }
+                }
+            }
+        }
+
+        Mono { vars: result }
+    }
 }
 
 impl<V: Var, T: IntoIterator> From<T> for Mono<V>
@@ -75,19 +194,19 @@ where
 
 impl<V: Var, U: Into<V>> FromIterator<(U, u32)> for Mono<V> {
     fn from_iter<T: IntoIterator<Item = (U, u32)>>(iter: T) -> Self {
-        let mut term = iter
+        let mut vars = iter
             .into_iter()
             .map(|(t, pow)| (t.into(), pow))
             .collect::<Vec<(V, u32)>>();
-        term.sort();
+        vars.sort();
 
         // Check duplicate variables
         assert!(
-            (term[..])
+            (vars[..])
                 .windows(2)
                 .all(|window| window[0].0 != window[1].0)
         );
 
-        Mono { term }
+        Mono { vars }
     }
 }

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 {
@@ -33,19 +34,37 @@ impl<V: Var> Display for Poly<V> {
                 "{}",
                 self.mono
                     .iter()
-                    .map(|(m, c)| format!("{}{}", c, m))
+                    .map(|(m, c)| match c {
+                        1 => format!("{}", m),
+                        -1 => format!("-{}", m),
+                        _ => format!("{}{}", c, m),
+                    })
                     .join(" + ")
             ),
         }
     }
 }
 
+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!(
             fmt,
             "{}",
-            self.term
+            self.vars
                 .iter()
                 .map(|(t, p)| match p {
                     1 => format!("{t}"),

src/poly/ideal.rs

@@ -0,0 +1,69 @@
+use std::marker::PhantomData;
+
+use super::buchberger;
+use super::flat::Poly;
+use super::var::Var;
+
+/// Marker: the ideal's generators are arbitrary polynomials.
+#[derive(Clone, Debug)]
+pub struct Generators;
+
+/// Marker: the ideal's generators form a Gröbner basis.
+#[derive(Clone, Debug)]
+pub struct GroebnerBasis;
+
+#[derive(Clone, Debug)]
+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

@@ -24,31 +24,31 @@ impl<V: Var> Mul for Mono<V> {
     type Output = Self;
 
     fn mul(self, other: Mono<V>) -> Self::Output {
-        let mut a_term = self.term.into_iter().peekable();
+        let mut a_vars = self.vars.into_iter().peekable();
-        let mut b_term = other.term.into_iter().peekable();
+        let mut b_vars = other.vars.into_iter().peekable();
 
         let mut result: Vec<(V, u32)> = Default::default();
 
         loop {
-            match (a_term.peek(), b_term.peek()) {
+            match (a_vars.peek(), b_vars.peek()) {
                 (Some((a_var, _)), Some((b_var, _))) => {
                     if a_var < b_var {
-                        result.push(a_term.next().unwrap());
+                        result.push(a_vars.next().unwrap());
                     } else if a_var > b_var {
-                        result.push(b_term.next().unwrap());
+                        result.push(b_vars.next().unwrap());
                     } else {
-                        let (var, a_exp) = a_term.next().unwrap();
+                        let (var, a_exp) = a_vars.next().unwrap();
-                        let (_, b_exp) = b_term.next().unwrap();
+                        let (_, b_exp) = b_vars.next().unwrap();
                         result.push((var, a_exp + b_exp));
                     }
                 }
                 (Some(a), None) => {
                     result.push(a.clone());
-                    a_term.next();
+                    a_vars.next();
                 }
                 (None, Some(b)) => {
                     result.push(b.clone());
-                    b_term.next();
+                    b_vars.next();
                 }
                 (None, None) => {
                     break;
@@ -56,7 +56,7 @@ impl<V: Var> Mul for Mono<V> {
             }
         }
 
-        Mono { term: result }
+        Mono { vars: result }
     }
 }
 
@@ -72,6 +72,20 @@ impl<V: Var> Mul<Mono<V>> for i32 {
     }
 }
 
+impl<V: Var> Mul<Poly<V>> for i32 {
+    type Output = Poly<V>;
+
+    fn mul(self, mut poly: Poly<V>) -> Self::Output {
+        if self == 0 {
+            return Poly::default();
+        }
+        for coeff in poly.mono.values_mut() {
+            *coeff *= self;
+        }
+        poly
+    }
+}
+
 impl<V: Var> Mul for Poly<V> {
     type Output = Poly<V>;
 
@@ -113,3 +127,79 @@ impl<V: Var> Sub for Poly<V> {
         self
     }
 }
+
+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.
+    ///
+    /// Returns `(d, q, r)` satisfying `lc(divisor)^d * self = q * divisor + r`,
+    /// where no term of `r` has its monomial divisible by `LM(divisor)` under lex order.
+    ///
+    /// 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.leading_term_lex().expect("divisor must be nonzero");
+
+        let mut p = self;
+        let mut q = Poly::default();
+        let mut r = Poly::default();
+        let mut d = 0u32;
+
+        while !p.is_zero() {
+            let (lt_p_mono, lt_p_coeff) = p.leading_term_lex().unwrap();
+
+            if lt_p_mono.contains(&lt_g_mono) {
+                let t_mono = lt_p_mono.div(&lt_g_mono);
+                if lt_p_coeff % lt_g_coeff == 0 {
+                    // Exact division: no need to multiply through by lc(g)
+                    let t_poly: Poly<V> = (lt_p_coeff / lt_g_coeff) * t_mono;
+                    p = p - t_poly.clone() * divisor.clone();
+                    q = q + t_poly;
+                } else {
+                    // Pseudo-division: multiply through by lc(g) to stay in ℤ
+                    let t_poly: Poly<V> = lt_p_coeff * t_mono;
+                    p = lt_g_coeff * p - t_poly.clone() * divisor.clone();
+                    q = lt_g_coeff * q + t_poly;
+                    r = lt_g_coeff * r;
+                    d += 1;
+                }
+            } else {
+                let lt_poly: Poly<V> = lt_p_coeff * lt_p_mono;
+                p = p - lt_poly.clone();
+                r = r + lt_poly;
+            }
+        }
+
+        (d, q, r)
+    }
+}

src/poly/tests.rs

@@ -1,4 +1,6 @@
-use super::flat::{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]
@@ -55,7 +57,7 @@ fn test_mono_mul() {
 
     // Multiply by constant monomial (empty term vec = 1)
     let a: Mono<StaticVar> = [("x", 4)].into();
-    let one: Mono<StaticVar> = Mono { term: vec![] };
+    let one: Mono<StaticVar> = Mono { vars: vec![] };
     assert_eq!(a.clone() * one, a);
 }
 
@@ -99,3 +101,364 @@ fn test_poly_sub() {
     let b: Poly<StaticVar> = [(4, [("x", 2)]), (1, [("y", 1)])].into();
     assert_eq!(a - b, Poly::default());
 }
+
+#[test]
+fn test_lex_cmp() {
+    use std::cmp::Ordering;
+
+    let x2: Mono<StaticVar> = [("x", 2)].into();
+    let xy: Mono<StaticVar> = [("x", 1), ("y", 1)].into();
+    let y2: Mono<StaticVar> = [("y", 2)].into();
+    let x: Mono<StaticVar> = [("x", 1)].into();
+    let one: Mono<StaticVar> = Mono { vars: vec![] };
+
+    // x² > xy (x exponent 2 vs 1)
+    assert_eq!(lex_cmp(&x2, &xy), Ordering::Greater);
+    // x > y² (x has higher priority)
+    assert_eq!(lex_cmp(&x, &y2), Ordering::Greater);
+    // xy > y² (x present in xy but not y²)
+    assert_eq!(lex_cmp(&xy, &y2), Ordering::Greater);
+    // 1 < x
+    assert_eq!(lex_cmp(&one, &x), Ordering::Less);
+    // reflexive
+    assert_eq!(lex_cmp(&x2, &x2), Ordering::Equal);
+}
+
+#[test]
+fn test_mono_div() {
+    // x² / x = x
+    let a: Mono<StaticVar> = [("x", 2)].into();
+    let b: Mono<StaticVar> = [("x", 1)].into();
+    assert_eq!(a.div(&b), Mono::from([("x", 1)]));
+
+    // x²y / xy = x
+    let a: Mono<StaticVar> = [("x", 2), ("y", 1)].into();
+    let b: Mono<StaticVar> = [("x", 1), ("y", 1)].into();
+    assert_eq!(a.div(&b), Mono::from([("x", 1)]));
+
+    // x²y / y = x²
+    let a: Mono<StaticVar> = [("x", 2), ("y", 1)].into();
+    let b: Mono<StaticVar> = [("y", 1)].into();
+    assert_eq!(a.div(&b), Mono::from([("x", 2)]));
+
+    // x / x = 1
+    let a: Mono<StaticVar> = [("x", 1)].into();
+    let b: Mono<StaticVar> = [("x", 1)].into();
+    assert_eq!(a.div(&b), Mono { vars: vec![] });
+}
+
+#[test]
+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);
+}
+
+#[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 { vars: 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;
+    let rhs = q * g.clone() + r;
+    assert_eq!(lhs, rhs);
+}
+
+#[test]
+fn test_div_rem() {
+    // x³ / x² = x, r=0, d=0
+    let f: Poly<StaticVar> = [(1, [("x", 3)])].into();
+    let g: Poly<StaticVar> = [(1, [("x", 2)])].into();
+    let expected_q: Poly<StaticVar> = [(1, [("x", 1)])].into();
+    let (d, q, r) = f.clone().div_rem(&g);
+    assert_eq!(q, expected_q);
+    assert!(r.is_zero());
+    verify_div_rem(f, &g, d, q, r);
+
+    // (x³ + x²y) / x² = x + y, r=0
+    // f = x²(x + y), g = x²
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 3u32)])),
+        (1i32, Mono::from([("x", 2u32), ("y", 1u32)])),
+    ]
+    .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);
+    assert_eq!(q, expected_q);
+    assert!(r.is_zero());
+    verify_div_rem(f, &g, d, q, r);
+
+    // (x³ + y) / x² = x, r=y
+    // LT(x³) divisible by x², LT(y) is not
+    let f: Poly<StaticVar> = [(1, [("x", 3)]), (1, [("y", 1)])].into();
+    let g: Poly<StaticVar> = [(1, [("x", 2)])].into();
+    let expected_q: Poly<StaticVar> = [(1, [("x", 1)])].into();
+    let expected_r: Poly<StaticVar> = [(1, [("y", 1)])].into();
+    let (d, q, r) = f.clone().div_rem(&g);
+    assert_eq!(q, expected_q);
+    assert_eq!(r, expected_r);
+    verify_div_rem(f, &g, d, q, r);
+
+    // 3x² / (2x): 2 ∤ 3, needs pseudo-division
+    // 2¹ · 3x² = 3x · 2x  =>  d=1, q=3x, r=0
+    let f: Poly<StaticVar> = [(3, [("x", 2)])].into();
+    let g: Poly<StaticVar> = [(2, [("x", 1)])].into();
+    let expected_q: Poly<StaticVar> = [(3, [("x", 1)])].into();
+    let (d, q, r) = f.clone().div_rem(&g);
+    assert_eq!(d, 1);
+    assert_eq!(q, expected_q);
+    assert!(r.is_zero());
+    verify_div_rem(f, &g, d, q, r);
+
+    // 6xy / 2 = 3xy, r=0, d=0
+    let f: Poly<StaticVar> = [(6, [("x", 1), ("y", 1)])].into();
+    let g = make_const_poly(2);
+    let expected_q: Poly<StaticVar> = [(3, [("x", 1), ("y", 1)])].into();
+    let (d, q, r) = f.clone().div_rem(&g);
+    assert_eq!(d, 0);
+    assert_eq!(q, expected_q);
+    assert!(r.is_zero());
+    verify_div_rem(f, &g, d, q, r);
+
+    // (x² + xy) / (x + y) = x, r=0
+    // f = x(x + y), g = x + y
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32)])),
+        (1i32, Mono::from([("x", 1u32), ("y", 1u32)])),
+    ]
+    .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);
+    assert_eq!(q, expected_q);
+    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())
+        );
+    }
+}