diff --git a/src/poly/buchberger.rs b/src/poly/buchberger.rs
index b4fc073..f5817f7 100644
--- a/src/poly/buchberger.rs
+++ b/src/poly/buchberger.rs
@@ -24,6 +24,30 @@ pub fn groebner_basis<V: Var>(generators: Vec<Poly<V>>) -> Vec<Poly<V>> {
         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
 }
 
@@ -47,7 +71,7 @@ pub fn is_groebner_basis<V: Var>(basis: &[Poly<V>]) -> bool {
 /// by the leading monomial of any element in `basis`.
 ///
 /// Uses the pseudo-division remainder and repeats until stable.
-fn reduce<V: Var>(f: &Poly<V>, basis: &[Poly<V>]) -> Poly<V> {
+pub(crate) fn reduce<V: Var>(f: &Poly<V>, basis: &[Poly<V>]) -> Poly<V> {
     let mut p = f.clone();
 
     'outer: loop {
diff --git a/src/poly/fmt.rs b/src/poly/fmt.rs
index 29967af..ab59df8 100644
--- a/src/poly/fmt.rs
+++ b/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!(
diff --git a/src/poly/ideal.rs b/src/poly/ideal.rs
new file mode 100644
index 0000000..718148e
--- /dev/null
+++ b/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()
+    }
+}
diff --git a/src/poly/mod.rs b/src/poly/mod.rs
index 94f9647..4cc7c33 100644
--- a/src/poly/mod.rs
+++ b/src/poly/mod.rs
@@ -1,6 +1,7 @@
 pub mod buchberger;
 pub mod flat;
 pub mod fmt;
+pub mod ideal;
 pub mod ops;
 pub mod var;
 
diff --git a/src/poly/tests.rs b/src/poly/tests.rs
index 7a9536d..d29c155 100644
--- a/src/poly/tests.rs
+++ b/src/poly/tests.rs
@@ -1,5 +1,6 @@
 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]
@@ -377,3 +378,87 @@ fn test_is_groebner_basis() {
     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())
+        );
+    }
+}