diff --git a/src/algorithms/hilbert.jl b/src/algorithms/hilbert.jl
index a42103e..9f86d9b 100644
--- a/src/algorithms/hilbert.jl
+++ b/src/algorithms/hilbert.jl
@@ -1,10 +1,15 @@
 @doc Markdown.doc"""
-    hilbert_series(I::Ideal{T}) where T <: MPolyRingElem
+    hilbert_series(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
 Compute the Hilbert series of a given polynomial ideal `I`.
 
-Based on: Anna M. Bigatti, Computation of Hilbert-Poincaré series,
-Journal of Pure and Applied Algebra, 1997.
+If `use_mdd == false` then the method is based on: Anna M. Bigatti,
+Computation of Hilbert-Poincaré series, Journal of Pure and Applied
+Algebra, 1997.
+
+Otherwise monomial divisibility diagrams are used, as defined in
+Pierre Lairez, Rafael Mohr, Théo Ternier, A data structure for
+monomial ideals with applications to signature Gröbner bases.
 
 **Notes**:
 * This requires a Gröbner basis of `I`, which is computed internally if not already known.
@@ -22,20 +27,31 @@ julia> hilbert_series(I)
 (-2*t - 1)//(t - 1)
 ```
 """
-function hilbert_series(I::Ideal{T}) where T <: MPolyRingElem
+function hilbert_series(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
     gb = get!(I.gb, 0) do
         groebner_basis(I, complete_reduction = true)
     end
     lead_exps = [ _lead_exp_ord(g, :degrevlex) for g in gb if !iszero(g) ]
-    return _hilbert_series_mono(lead_exps, nvars(parent(I)))
+    if use_mdd
+        n = ngens(parent(I))
+        diagram = create_diagram([monomial(SVector{n}(e)) for e in lead_exps])
+        return hilbert_series_mdd(diagram)
+    else
+        return _hilbert_series_mono(lead_exps, nvars(parent(I)))
+    end
 end
 
 @doc Markdown.doc"""
-    hilbert_degree(I::Ideal{T}) where T <: MPolyRingElem
+    hilbert_degree(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
 Compute the degree of a given polynomial ideal `I` by first computing its Hilbert series.
 
+If `use_mdd == true`, monomial divisibility diagrams are used to
+compute the Hilbert series, as defined in Pierre Lairez, Rafael Mohr,
+Théo Ternier, A data structure for monomial ideals with applications
+to signature Gröbner bases.
+
 **Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
 
 # Examples
@@ -50,18 +66,23 @@ julia> hilbert_degree(I)
 3
 ```
 """
-function hilbert_degree(I::Ideal{T}) where T <: MPolyRingElem
+function hilbert_degree(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
     !isnothing(I.deg) && return I.deg
-    I.deg = numerator(hilbert_series(I))(1) |> abs
+    I.deg = numerator(hilbert_series(I, use_mdd = use_mdd))(1) |> abs
     return I.deg
 end
 
 @doc Markdown.doc"""
-    hilbert_dimension(I::Ideal{T}) where T <: MPolyRingElem
+    hilbert_dimension(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
 Compute the Krull dimension of a given polynomial ideal `I` by first computing its Hilbert series.
 
+If `use_mdd == true`, monomial divisibility diagrams are used to
+compute the Hilbert series, as defined in Pierre Lairez, Rafael Mohr,
+Théo Ternier, A data structure for monomial ideals with applications
+to signature Gröbner bases.
+
 **Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
 
 # Examples
@@ -76,15 +97,15 @@ julia> hilbert_dimension(I)
 1
 ```
 """
-function hilbert_dimension(I::Ideal{T}) where T <: MPolyRingElem
+function hilbert_dimension(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
-    H = hilbert_series(I)
+    H = hilbert_series(I, use_mdd = use_mdd)
     I.dim = iszero(H) ? -1 : degree(denominator(H))
     return I.dim
 end
 
 @doc Markdown.doc"""
-    hilbert_polynomial(I::Ideal{T}) where T <: MPolyRingElem
+    hilbert_polynomial(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
 Compute the Hilbert polynomial and the index of regularity of a given polynomial ideal `I`
 by first computing its Hilbert series. The index of regularity is the smallest integer such that
@@ -93,6 +114,11 @@ the Hilbert function and polynomial match.
 Note that the Hilbert polynomial of I has leading term (e/d!)*t^d, where e and d are respectively
 the degree and Krull dimension of I.
 
+If `use_mdd == true`, monomial divisibility diagrams are used to
+compute the Hilbert series, as defined in Pierre Lairez, Rafael Mohr,
+Théo Ternier, A data structure for monomial ideals with applications
+to signature Gröbner bases.
+
 **Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
 
 # Examples
@@ -107,10 +133,10 @@ julia> hilbert_polynomial(I)
 (3*s + 3, 1)
 ```
 """
-function hilbert_polynomial(I::Ideal{T}) where T <: MPolyRingElem
+function hilbert_polynomial(I::Ideal{T}; use_mdd = false) where T <: MPolyRingElem
 
     A, s = polynomial_ring(QQ, :s)
-    H = hilbert_series(I)
+    H = hilbert_series(I, use_mdd = use_mdd)
     dim = degree(denominator(H))
     num = iseven(dim) ? numerator(H) : -numerator(H)
     dim==0 && return num(s), 0
diff --git a/src/siggb/helpers.jl b/src/siggb/helpers.jl
index 4ee69fc..f5c92c2 100644
--- a/src/siggb/helpers.jl
+++ b/src/siggb/helpers.jl
@@ -10,6 +10,9 @@ function new_basis(basis_size, syz_size,
     sigratios = Vector{Monomial{N}}(undef, basis_size)
     rewrite_nodes = Vector{Vector{Int}}(undef, basis_size+1)
     lm_masks = Vector{DivMask}(undef, basis_size)
+    koszul_diagram = EMPTY_DIAGRAM
+    lm_diagram = EMPTY_DIAGRAM
+    hashstate = new_hashstate()
     monomials = Vector{Vector{MonIdx}}(undef, basis_size)
     coeffs = Vector{Vector{Coeff}}(undef, basis_size)
     is_red = Vector{Bool}(undef, basis_size)
@@ -18,7 +21,8 @@ function new_basis(basis_size, syz_size,
     syz_sigs = Vector{Monomial{N}}(undef, syz_size)
     syz_masks = Vector{MaskSig}(undef, syz_size)
     basis = Basis(sigs, sigmasks, sigratios, rewrite_nodes,
-                  lm_masks, monomials, coeffs, is_red,
+                  lm_masks, lm_diagram, koszul_diagram, hashstate,
+                  monomials, coeffs, is_red,
                   mod_rep_known, mod_reps,
                   syz_sigs, syz_masks, Exp[],
                   input_length,
@@ -70,7 +74,7 @@ function make_room_new_input_el!(basis::Basis,
             basis.mod_rep_known[i] = basis.mod_rep_known[i-shift]
             basis.mod_reps[i] = basis.mod_reps[i-shift]
 
-            # adjust rewrite tree
+            # adjust rewrite diagram
             rnodes = basis.rewrite_nodes[i-shift+1]
             if rnodes[2] >= old_offset + 1
                 rnodes[2] += shift
@@ -373,7 +377,7 @@ function sort_pairset!(pairset::Pairset, from::Int, sz::Int,
 end
 
 function new_timer()
-    return Timings(0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0)
+    return Timings(0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0, 0.0, 0.0)
 end
 
 function Base.show(io::IO, timer::Timings)
@@ -382,6 +386,8 @@ function Base.show(io::IO, timer::Timings)
     @printf io "linear algebra:      %.2f\n" timer.lin_alg_time
     @printf io "select:              %.2f\n" timer.select_time
     @printf io "update:              %.2f\n" timer.update_time
+    @printf io "mdd creation:        %.2f\n" timer.time_for_mdd
+    @printf io "membership tests:    %.2f\n" timer.time_for_membership
     !iszero(timer.module_time) && @printf io "module construction: %.2f\n" timer.module_time
     !iszero(timer.comp_lc_time) && @printf io "splitting:           %.2f\n" timer.comp_lc_time
 end
diff --git a/src/siggb/monomial_diagram.jl b/src/siggb/monomial_diagram.jl
new file mode 100644
index 0000000..eb55210
--- /dev/null
+++ b/src/siggb/monomial_diagram.jl
@@ -0,0 +1,328 @@
+# Function to hash a diagram.
+function Base.hash(diagram::Diagram, h::UInt)
+    return hash(diagram.id, h)
+end
+
+# Function to test equality between two diagrams in O(1).
+function Base.:(==)(d1::Diagram, d2::Diagram)
+    return d1.id == d2.id
+end
+
+function Base.show(io::IO, diagram::Diagram)
+    print(diagram.edges)
+end
+
+# Function to print the diagram as an n-tree.
+function print_diagram(diagram::Diagram, space::Int=0)
+    for edge in diagram.edges
+        println(" " ^ space * " ", edge[1])
+        print_diagram(edge[2], space + 2)
+    end
+end
+
+# This function is to create a new hash state.
+function new_hashstate()
+    return HashState(Dict{Vector{Edge}, Diagram}(), 0, 0, 0)
+end
+
+# Function to troncate a list of edges.
+function truncate!(children::Vector{Edge}, i::Int)
+    if length(children) > i
+        resize!(children,i)
+    end
+end
+
+# Function to create a monomial divisibility diagram with hash consing
+function make_diagram(children::Vector{Edge}, hashstate::HashState)
+    i = 0
+
+    # In order to keep strict inclusions we remove the edges that are the same keeping the leftmost one.
+    last = EMPTY_DIAGRAM
+    for j in 1:length(children)
+        current_edge = children[j][2]
+        if current_edge != last
+            last = current_edge
+            children[i+1] = children[j]
+            i += 1
+        end
+    end
+
+    truncate!(children, i)
+    key = get(hashstate.hashtable, children, EMPTY_DIAGRAM)
+
+    if key != EMPTY_DIAGRAM
+        return hashstate.hashtable[children]
+    end
+
+    diagram = Diagram(hashstate.counter, children)
+    hashstate.hashtable[children] = diagram
+    hashstate.counter += 1
+    return diagram
+end
+
+# Function to insert a prefix in the diagram
+function insert_aux(diagram::Diagram, m::Monomial, i::Int, hashstate::HashState)
+    if isempty(m.exps) || i == 1
+        return diagram
+    end
+
+    i -= 1
+    a = m.exps[i]
+    new_diagram = Edge[]
+    inserted = false
+
+    for edge in diagram.edges
+        if edge[1] == a
+            push!(new_diagram, (a,insert_aux(edge[2], m, i, hashstate)))
+            inserted = true
+        else
+            push!(new_diagram, edge)
+        end
+    end
+    if !inserted
+        push!(new_diagram, (a,insert_aux(make_diagram(Edge[], hashstate),m,i,hashstate)))
+    end
+    return make_diagram(new_diagram, hashstate)
+end
+
+function insert(diagram::Diagram, m::Monomial, hashstate::HashState)
+    return insert_aux(diagram, m, length(m.exps)+1, hashstate)
+end
+
+
+# Function to insert a monomial in the diagram with memoization
+function insertion_aux(diagram::Diagram, m::Monomial, i::Int, memo::Dict{Tuple{Diagram,Int}, Diagram}, hashstate::HashState)
+    if diagram == EMPTY_DIAGRAM
+        return insert_aux(make_diagram(Edge[], hashstate), m, i, hashstate)
+    end
+    if i == 1
+         return diagram
+    end
+    if isempty(diagram.edges)
+        return diagram
+    end
+    key = get(memo, (diagram, i), EMPTY_DIAGRAM)
+    if key != EMPTY_DIAGRAM
+        return key
+    end
+
+    i -= 1
+    new_diagram = Edge[]
+    is_sub_diagram = false
+    last_sub_diagram = EMPTY_DIAGRAM
+    for edge in diagram.edges
+        if edge[1] < m.exps[i]
+            push!(new_diagram, edge)
+            last_sub_diagram = edge[2]
+        else
+            if edge[1] == m.exps[i]
+                is_sub_diagram = true
+            end
+            if edge[1] > m.exps[i] && !is_sub_diagram
+                push!(new_diagram, (m.exps[i], insertion_aux(last_sub_diagram, m, i, memo, hashstate)))
+                is_sub_diagram = true
+            end
+
+            push!(new_diagram, (edge[1], insertion_aux(edge[2], m, i, memo, hashstate)))
+        end
+    end
+
+    if !is_sub_diagram
+        push!(new_diagram, (m.exps[i], insertion_aux(last_sub_diagram, m, i, memo, hashstate)))
+    end
+
+    final_diagram = make_diagram(new_diagram, hashstate)
+    memo[(diagram, i)] = final_diagram
+    return final_diagram
+end
+
+function insertion(diagram::Diagram, m::Monomial, hashstate::HashState)
+    return insertion_aux(diagram, m, length(m.exps)+1, Dict{Tuple{Diagram,Int}, Diagram}(), hashstate)
+end
+
+
+# Function to create the monomial divisibility diagram of a list of monomials
+function create_diagram(monomial_list::Vector{<:Monomial}, hashstate::HashState = new_hashstate())
+    diagram = EMPTY_DIAGRAM
+    for m in monomial_list
+        diagram = insertion(diagram, m, hashstate)
+    end
+    return diagram
+end
+
+# Function that determines the number of nodes in the monomial divisibility diagram
+function size_of_diagram(diagram::Diagram)
+    if isempty(diagram.edges) || diagram == EMPTY_DIAGRAM
+        return 0
+    end
+    return 1 + sum(size_of_diagram(sub_diagram[2]) for sub_diagram in diagram.edges)
+end
+
+# Function that determines the number of distinct nodes in the monomial divisibility diagram
+function number_of_distinct_nodes(diagram::Diagram, mem::Dict{Diagram, Diagram} = Dict{Diagram,Diagram}())
+    if isempty(diagram.edges) || diagram == EMPTY_DIAGRAM
+        return 0
+    end
+    if haskey(mem, diagram)
+        return 0
+    end
+    number = 0
+    for sub_diagram in diagram.edges
+        if !haskey(mem, sub_diagram[2])
+            number += 1 + number_of_distinct_nodes(sub_diagram[2], mem)
+        end
+    end
+    mem[diagram] = diagram
+    return number
+end
+
+# Function that determines the largest s such that s <= exp, returns -1 otherwise
+function find_nearest_index(sub_diagram::Vector{Tuple{Exp, Diagram}}, exp::Exp)
+    j = -1
+    for k in 1:length(sub_diagram)
+        if sub_diagram[k][1] > exp
+            return j
+        end
+        j = k
+    end
+    return j
+end
+
+# Function that test if a monomial is represented by a monomial divisibility diagram
+function is_in_diagram(m::Monomial{N}, diagram::Diagram) where N
+    sub_diagram = diagram
+
+    @inbounds for j in N:-1:1
+        exp = m.exps[j]
+        i = find_nearest_index(sub_diagram.edges, exp)
+        if i == -1
+            return false
+        else
+            sub_diagram = sub_diagram.edges[i][2]
+        end
+    end
+    return true
+end
+
+# Function that converts a monomial::MPolyRingElem into an object of type Monomial{N}
+function convert_to_monomial(m::MPolyRingElem, R::MPolyRing, ::Val{N}) where N
+    exts = first(collect(exponent_vectors(R(m))))
+    sv = SVector{N, Int}(exts)
+    return monomial(sv)
+end
+
+# Function that tests if a divides b
+function divides(a::Monomial, b::Monomial)
+    return all(ai <= bi for (ai, bi) in zip(a.exps, b.exps))
+end
+
+# Function that generates a random list of monomials
+function generate_random_vectors(r::Int, d::Int, n::Int)
+    N = 0
+    list_of_monomials = Monomial[]
+    while N < n
+        rd = monomial(SVector{r}(rand(Exp(0):Exp(d-1), r)))
+        if all(!divides(mon, rd) for mon in list_of_monomials)
+            push!(list_of_monomials, rd)
+        end
+        N += 1
+    end
+
+    return list_of_monomials
+end
+
+# Function that naively test if a monomial is in a monomial ideal
+function naive_is_in_ideal(mon::Monomial, list::Vector{Monomial})
+    for i in 1:length(list)
+        if divides(list[i], mon)
+            return true
+        end
+    end
+    return false
+end
+
+# Function that returns the list of all leaves of a monomial divisibility diagram
+function get_leaves(diagram::Diagram)
+    if isempty(diagram.edges)
+        return Vector{Exp}[Exp[]]
+    end
+
+    leaves = Vector{Exp}[]
+    for (label, sub_diagram) in diagram.edges
+        for path in get_leaves(sub_diagram)
+            push!(leaves, vcat([label],path))
+        end
+    end
+    return leaves
+end
+
+function depth(diagram::Diagram)::Int
+    if isempty(diagram.edges)
+        return 0
+    end
+
+    return 1 + depth(diagram.edges[1][2])
+end
+
+function hilbert_series_mdd_aux(R, i::Int, diagram::Diagram)
+    t = gens(R)[1]
+
+    if i == 1
+        somme = R(0)
+        for j in 0:(diagram.edges[1][1]-1)
+            somme += t^j
+        end
+        return (R(1)-t) * somme
+    end
+
+    hilbert_series = R(1) - t^(diagram.edges[1][1])
+    for j in 1:length(diagram.edges)
+        if j == length(diagram.edges)
+            hilbert_series += t^(diagram.edges[j][1]) * hilbert_series_mdd_aux(R, i-1, diagram.edges[j][2])
+        else
+            hilbert_series += (t^(diagram.edges[j][1]) - t^(diagram.edges[j+1][1])) * hilbert_series_mdd_aux(R, i-1, diagram.edges[j][2])
+        end
+    end
+
+    return hilbert_series
+end
+
+function hilbert_series_mdd(diagram::Diagram)
+    n = depth(diagram)
+    R, t = polynomial_ring(ZZ, :t)
+    K = fraction_field(R)
+    return K(hilbert_series_mdd_aux(R, n, diagram)) / (K(1)-K(t))^n
+end
+
+function multi_hilbert_series_mdd_aux(R, i::Int, diagram::Diagram)
+    vars = gens(R)
+    if i == 1
+        somme = R(0)
+        for j in 0:(diagram.edges[1][1]-1)
+            somme += vars[1]^j
+        end
+        return (R(1)-vars[1]) * somme
+    end
+
+    hilbert_series = 1 - vars[i]^(diagram.edges[1][1])
+    for j in 1:length(diagram.edges)
+        if j == length(diagram.edges)
+            hilbert_series += vars[i]^(diagram.edges[j][1]) * multi_hilbert_series_mdd_aux(R, i-1, diagram.edges[j][2])
+        else
+            hilbert_series += (vars[i]^(diagram.edges[j][1]) - vars[i]^(diagram.edges[j+1][1])) * multi_hilbert_series_mdd_aux(R, i-1, diagram.edges[j][2])
+        end
+    end
+
+    return hilbert_series
+end
+
+function multi_hilbert_series_mdd(diagram::Diagram)
+    n = depth(diagram)
+    R, vars = polynomial_ring(ZZ, ["x$i" for i in 1:n])
+    K = fraction_field(R)
+    prodn = K(1)
+    for x in vars
+        prodn *= K(R(1)-x)
+    end
+    return K(multi_hilbert_series_mdd_aux(R, n, diagram)) / prodn
+end
diff --git a/src/siggb/rewriting.jl b/src/siggb/rewriting.jl
index fccb54d..a8af64e 100644
--- a/src/siggb/rewriting.jl
+++ b/src/siggb/rewriting.jl
@@ -25,27 +25,38 @@ end
                                     sigmask::DivMask,
                                     ind_order::IndOrder,
                                     tags::Tags,
-                                    check::Bool)
+                                    mod_ord::Symbol,
+                                    check::Bool,
+                                    timer::Timings)
+
 
     if !check
         return false
     end
 
-    s_ind = index(sig)
-
-    @inbounds for i in basis.basis_offset:basis.basis_load
-        basis.is_red[i] && continue
-        b_ind = index(basis.sigs[i])
-        if (ind_order.ord[b_ind] < ind_order.ord[s_ind]
-            && !are_incompat(b_ind, s_ind, ind_order))
-            if divch(basis.lm_masks[i], sigmask)
-                if divch(leading_monomial(basis, basis_ht, i), monomial(sig))
-                    return true
+    if mod_ord == :DPOT
+        s_ind = index(sig)
+
+        @inbounds for i in basis.basis_offset:basis.basis_load
+            basis.is_red[i] && continue
+            b_ind = index(basis.sigs[i])
+            if (ind_order.ord[b_ind] < ind_order.ord[s_ind]
+                && !are_incompat(b_ind, s_ind, ind_order))
+                if divch(basis.lm_masks[i], sigmask)
+                    if divch(leading_monomial(basis, basis_ht, i), monomial(sig))
+                        return true
+                    end
                 end
             end
         end
+        return false
+    else
+        m = monomial(sig)
+        basis.hashstate.numberofmembershiptests += 1
+        tmt = @elapsed begin is_in = is_in_diagram(m, basis.koszul_diagram) end
+        timer.time_for_membership += tmt
+        return is_in
     end
-    return false
 end
 
 @inline function rewriteable_basis(basis::Basis,
@@ -71,13 +82,13 @@ function find_canonical_rewriter(basis::Basis,
                                  mod_ord::Symbol)
 
     if mod_ord == :DPOT
-        return find_canonical_rewriter_tree(basis, sig, sigmask)
+        return find_canonical_rewriter_diagram(basis, sig, sigmask)
     else
         return find_canonical_rewriter_rat(basis, sig, sigmask)
     end
 end
 
-function find_canonical_rewriter_tree(basis::Basis,
+function find_canonical_rewriter_diagram(basis::Basis,
                                       sig::Sig,
                                       sigmask::DivMask)
 
@@ -130,12 +141,13 @@ function rewriteable(basis::Basis,
                      ind_order::IndOrder,
                      tags::Tags,
                      check::Bool,
-                     mod_ord::Symbol)
+                     mod_ord::Symbol,
+                     timer::Timings)
 
     s_ind = index(sig)
 
     rewriteable_syz(basis, sig, sigmask, tags, check) && return true
     rewriteable_basis(basis, idx, sig, sigmask, tags, check, mod_ord) && return true
-    rewriteable_koszul(basis, basis_ht, sig, sigmask, ind_order, tags, check) && return true
+    rewriteable_koszul(basis, basis_ht, sig, sigmask, ind_order, tags, mod_ord, check, timer) && return true
     return false
 end
diff --git a/src/siggb/siggb.jl b/src/siggb/siggb.jl
index 2e59aa7..47cc308 100644
--- a/src/siggb/siggb.jl
+++ b/src/siggb/siggb.jl
@@ -19,6 +19,7 @@ include("normalform.jl")
 include("affine_cells.jl")
 include("interfaces.jl")
 include("helpers.jl")
+include("monomial_diagram.jl")
 
 
 #---------------- user functions --------------------#
@@ -94,6 +95,9 @@ function sig_groebner_basis(sys::Vector{T}; info_level::Int=0,
                                 mod_ord)
         @info "$(arit_ops) total submul's"
         @info timer
+        @info "Size of the mdd: $(number_of_distinct_nodes(basis.lm_diagram))"
+        @info "Number of insertions: $(basis.hashstate.numberofinsertion)"
+        @info "Number of tests: $(basis.hashstate.numberofmembershiptests)"
     end
 
     # output
@@ -151,7 +155,7 @@ function siggb!(basis::Basis{N},
                                                               ind_order, tags,
                                                               mod_ord)
         timer.select_time += tim
-        tim = @elapsed symbolic_pp!(basis, matrix, basis_ht, symbol_ht,
+        tim = @elapsed symbolic_pp!(timer, basis, matrix, basis_ht, symbol_ht,
                                     ind_order, tags, sigind, compat_ind,
                                     mod_ord)
         timer.sym_pp_time += tim
@@ -163,7 +167,7 @@ function siggb!(basis::Basis{N},
             arit_ops += arit_ops_new
             timer.lin_alg_time += tim
 
-            tim = @elapsed added_unit = update_siggb!(basis, matrix, pairset, symbol_ht,
+            tim = @elapsed added_unit = update_siggb!(timer, basis, matrix, pairset, symbol_ht,
                                                       basis_ht, ind_order, tags,
                                                       tr, char, syz_queue, mod_ord)
             timer.update_time += tim
@@ -204,6 +208,10 @@ function siggb!(basis::Basis{N},
             # minimize à la F5c
             min_idx = iszero(p_idx) ? zero(SigIndex) : curr_ind
             minimize!(basis, basis_ht, min_idx, ind_order, tags)
+
+            if !iszero(p_idx)
+                basis.koszul_diagram = basis.lm_diagram
+            end
         end
     end
     return false, arit_ops, nz_conds
@@ -339,7 +347,7 @@ function siggb_for_split!(basis::Basis{N},
         tim = @elapsed deg, _, _ = select_normal!(pairset, basis, matrix,
                                                   basis_ht, symbol_ht, ind_order, tags)
         timer.select_time += tim
-        tim = @elapsed symbolic_pp!(basis, matrix, basis_ht, symbol_ht,
+        tim = @elapsed symbolic_pp!(timer, basis, matrix, basis_ht, symbol_ht,
                                     ind_order, tags)
         timer.sym_pp_time += tim
 
@@ -348,7 +356,7 @@ function siggb_for_split!(basis::Basis{N},
         tim = @elapsed echelonize!(matrix, tags, ind_order, char, shift, tr)
         timer.lin_alg_time += tim
 
-        time = @elapsed update_siggb!(basis, matrix, pairset,
+        time = @elapsed update_siggb!(timer, basis, matrix, pairset,
                                       symbol_ht, basis_ht,
                                       ind_order, tags,
                                       tr, char, syz_queue)
diff --git a/src/siggb/symbolic_pp.jl b/src/siggb/symbolic_pp.jl
index 5271206..a92a2f7 100644
--- a/src/siggb/symbolic_pp.jl
+++ b/src/siggb/symbolic_pp.jl
@@ -187,7 +187,8 @@ function select_normal!(pairset::Pairset{N},
     return deg, compat_ind, sigind
 end
 
-function symbolic_pp!(basis::Basis{N},
+function symbolic_pp!(timer::Timings,
+                      basis::Basis{N},
                       matrix::MacaulayMatrix,
                       ht::MonomialHashtable,
                       symbol_ht::MonomialHashtable,
@@ -231,6 +232,12 @@ function symbolic_pp!(basis::Basis{N},
 
         exp = symbol_ht.exponents[i]
         divm = symbol_ht.hashdata[i].divmask
+
+        basis.hashstate.numberofmembershiptests += 1
+        if !is_in_diagram(exp, basis.lm_diagram)
+            i += one(MonIdx)
+            continue
+        end
         
         j = basis.basis_offset 
         @label target
@@ -293,7 +300,7 @@ function symbolic_pp!(basis::Basis{N},
             if rewriteable(basis, ht, j, mul_cand_sig,
                            cand_sig_mask,
                            ind_order, tags, true,
-                           mod_ord)
+                           mod_ord, timer)
                 j += 1
                 @goto target
             end
diff --git a/src/siggb/typedefs.jl b/src/siggb/typedefs.jl
index 48377b8..6193245 100644
--- a/src/siggb/typedefs.jl
+++ b/src/siggb/typedefs.jl
@@ -61,13 +61,35 @@ const MaskSig = Tuple{SigIndex, DivMask}
 const Polynomial = Tuple{Vector{Coeff}, Vector{MonIdx}}
 const ModCache{N} = Dict{Tuple{Sig{N}, SigIndex}, Polynomial}
 
+# For monomial diagram
+
+# Define a monomial structure for representing monomial ideals
+struct Diagram
+    id::Int
+    edges::Vector{Tuple{Exp, Diagram}}
+end
+
+# Define an Edge as a tuple of an exponent and a Diagram
+const Edge = Tuple{Exp, Diagram}
+
+# Define a mutable struct to hold the hash state
+mutable struct HashState
+    hashtable::Dict{Vector{Edge}, Diagram}
+    counter::Int64
+    numberofinsertion::Int64
+    numberofmembershiptests::Int64
+end
+
+# Define a constant for an empty diagram
+const EMPTY_DIAGRAM = Diagram(-1, Tuple{Int,Diagram}[])
+
 mutable struct Basis{N}
     sigs::Vector{Sig{N}}
     sigmasks::Vector{MaskSig}
 
     sigratios::Vector{Monomial{N}}
 
-    # tree structure:
+    # monomial divisibility diagram:
     #   - length of data for i is rewrite_nodes[i][1]
     #   - parent of i is rewrite_nodes[i][2]
     #   - children of i are rewrite_nodes[i][3:end]
@@ -76,6 +98,11 @@ mutable struct Basis{N}
 
     lm_masks::Vector{DivMask}
 
+    # monomial divisibility diagram
+    lm_diagram::Diagram
+    koszul_diagram::Diagram
+    hashstate::HashState
+
     monomials::Vector{Vector{MonIdx}}
     coefficients::Vector{Vector{Coeff}}
 
@@ -200,6 +227,7 @@ mutable struct LocClosedSet{T<:MPolyRingElem}
     gbs::Vector{Vector{T}}
 end
 
+
 # for benchmarking
 mutable struct Timings
     sym_pp_time::Float32
@@ -209,4 +237,6 @@ mutable struct Timings
     arit_ops::Int64
     module_time::Float32
     comp_lc_time::Float32
+    time_for_mdd::Float32
+    time_for_membership::Float32
 end
diff --git a/src/siggb/update.jl b/src/siggb/update.jl
index fc602f4..490b964 100644
--- a/src/siggb/update.jl
+++ b/src/siggb/update.jl
@@ -1,7 +1,8 @@
 # updating the pairset and basis
 
 # add new reduced rows to basis/syzygies
-function update_siggb!(basis::Basis,
+function update_siggb!(timer::Timings,
+                       basis::Basis,
                        matrix::MacaulayMatrix,
                        pairset::Pairset{N},
                        symbol_ht::MonomialHashtable,
@@ -45,7 +46,7 @@ function update_siggb!(basis::Basis,
             added_unit = add_basis_elem!(basis, pairset, basis_ht, symbol_ht,
                                          row, coeffs,
                                          new_sig, new_sig_mask, parent_ind,
-                                         tr, ind_order, tags, mod_ord)
+                                         tr, ind_order, tags, mod_ord, timer)
         end
     end
 
@@ -71,7 +72,8 @@ function add_basis_elem!(basis::Basis{N},
                          tr::Tracer,
                          ind_order::IndOrder,
                          tags::Tags,
-                         mod_ord::Symbol) where N
+                         mod_ord::Symbol,
+                         timer::Timings) where N
 
 
     # make sure we have enough space
@@ -80,6 +82,10 @@ function add_basis_elem!(basis::Basis{N},
     # add to basis hashtable
     insert_in_basis_hash_table_pivots!(row, basis_ht, symbol_ht)
     lm = basis_ht.exponents[first(row)]
+    # add to lm diagramn
+    timeinsertion = @elapsed begin basis.lm_diagram = insertion(basis.lm_diagram, lm, basis.hashstate) end
+    timer.time_for_mdd += timeinsertion
+    basis.hashstate.numberofinsertion += 1
     @debug "new lm $(lm)"
     lm_mask = divmask(lm, basis_ht.divmap, basis_ht.ndivbits)
     s = new_sig
@@ -105,18 +111,18 @@ function add_basis_elem!(basis::Basis{N},
 
     basis.is_red[l] = false
 
-    tree_data = basis.rewrite_nodes[parent_ind+1]
+    diagram_data = basis.rewrite_nodes[parent_ind+1]
     insind = 3 
-    @inbounds for j in insind:insind+tree_data[1]
-        child_ind = tree_data[j]
+    @inbounds for j in insind:insind+diagram_data[1]
+        child_ind = diagram_data[j]
         rat = basis.sigratios[child_ind-1]
         if lt_drl(new_sig_ratio, rat)
             break
         end
         insind += 1
     end
-    insert!(tree_data, insind, l+1)
-    tree_data[1] += 1
+    insert!(diagram_data, insind, l+1)
+    diagram_data[1] += 1
 
     # if an existing sig further reduced we dont need the old element
     # if basis.sigs[parent_ind] == new_sig && parent_ind >= basis.basis_offset
@@ -130,7 +136,7 @@ function add_basis_elem!(basis::Basis{N},
     store_basis_elem!(tr, new_sig, l, basis.basis_size)
     
     # build new pairs
-    update_pairset!(pairset, basis, basis_ht, l, ind_order, tags, mod_ord)
+    update_pairset!(timer, pairset, basis, basis_ht, l, ind_order, tags, mod_ord)
 
     return false
 end
@@ -250,7 +256,8 @@ end
 # construct all pairs with basis element at new_basis_idx
 # and perform corresponding rewrite checks
 # assumes DPOT
-function update_pairset!(pairset::Pairset{N},
+function update_pairset!(timer::Timings,
+                         pairset::Pairset{N},
                          basis::Basis,
                          basis_ht::MonomialHashtable,
                          new_basis_idx::Int,
@@ -335,11 +342,11 @@ function update_pairset!(pairset::Pairset{N},
                         basis_pair_sig_mask, tags,
                         mod_ord == :DPOT || cmp_ind(new_sig_idx, basis_sig_idx, ind_order)) && continue
         rewriteable_koszul(basis, basis_ht, new_pair_sig,
-                           new_pair_sig_mask, ind_order, tags,
-                           mod_ord == :DPOT || cmp_ind(basis_sig_idx, new_sig_idx, ind_order)) && continue
+                           new_pair_sig_mask, ind_order, tags, mod_ord,
+                           mod_ord == :DPOT || cmp_ind(basis_sig_idx, new_sig_idx, ind_order), timer) && continue
         rewriteable_koszul(basis, basis_ht, basis_pair_sig,
-                           basis_pair_sig_mask, ind_order, tags,
-                           mod_ord == :DPOT || cmp_ind(new_sig_idx, basis_sig_idx, ind_order)) && continue
+                           basis_pair_sig_mask, ind_order, tags, mod_ord,
+                           mod_ord == :DPOT || cmp_ind(new_sig_idx, basis_sig_idx, ind_order), timer) && continue
 
         top_sig, top_sig_mask, top_index,
         bot_sig, bot_sig_mask, bot_index = begin
diff --git a/test/algorithms/decomposition.jl b/test/algorithms/decomposition.jl
index e8b7981..fd81884 100644
--- a/test/algorithms/decomposition.jl
+++ b/test/algorithms/decomposition.jl
@@ -4,7 +4,7 @@
          3038*x^2*y-3686*x^2*z+2288*x*y^2-16544*x*y*z-3344*x*z^2+315*y^3+4111*y^2*z+157*y*z^2-663*z^3-27166*x*y+4168*x*z-2769*y^2+13942*y*z+1527*z^2+25806*y-690*z,
          650*x^2*y+3350*x^2*z+195*x*y^2-8450*x*y*z+845*x*z^2-260*y^3-3410*y^2*z-390*y*z^2+13390*x*y-1630*x*z-276*y^2+8372*y*z-15088*y,
          -225*x^2*y+1685*x^2*z-705*x*y^2-450*x*y*z+345*x*z^2-420*y^3+1325*y^2*z+645*y*z^2+8056*x*y-705*x*z+918*y^2+1527*y*z-5658*y]
-    I = Ideal(F)
+    I = AlgebraicSolving.Ideal(F)
     dec = equidimensional_decomposition(I)
     @test AlgebraicSolving._check_decomp(I, dec)
 
@@ -13,7 +13,7 @@
          10816588*x1^3+38077162*x1^2*x2+25140479*x1^2*x3-124056256*x1^2*x4-193635828*x1*x2^2-6790460*x1*x2*x3+250307870*x1*x2*x4-10319588*x1*x3^2+5839912*x1*x3*x4+112915639*x1*x4^2+162291312*x2^3-126154506*x2^2*x3-52085058*x2^2*x4+133237254*x2*x3^2-122040686*x2*x3*x4-132982128*x2*x4^2+16298428*x3^3+46623254*x3^2*x4-26720875*x3*x4^2-16047619*x4^3-9400295*x1^2-82452832*x1*x2+22255617*x1*x3+115474646*x1*x4+157662837*x2^2-2637098*x2*x3-183748338*x2*x4+119183502*x3^2-131012023*x3*x4-88321124*x4^2-17533608*x1+40968512*x2+10206682*x3-58141098*x4+12230767,
          -19696652*x1^3+25140479*x1^2*x2-38705922*x1^2*x3+30086207*x1^2*x4-3395230*x1*x2^2-20639176*x1*x2*x3+5839912*x1*x2*x4+36342096*x1*x3^2+64629228*x1*x3*x4-5646418*x1*x4^2-42051502*x2^3+133237254*x2^2*x3-61020343*x2^2*x4+48895284*x2*x3^2+93246508*x2*x3*x4-26720875*x2*x4^2+60295600*x3^3+28934928*x3^2*x4-19580172*x3*x4^2-9372130*x4^3-1289811*x1^2+22255617*x1*x2+109167336*x1*x3+150465505*x1*x4-1318549*x2^2+238367004*x2*x3-131012023*x2*x4+124816848*x3^2+116104228*x3*x4-70772721*x4^2+133696391*x1+10206682*x2+145409424*x3-111938286*x4-23665371,
          12014836*x1^3-124056256*x1^2*x2+30086207*x1^2*x3-31723130*x1^2*x4+125153935*x1*x2^2+5839912*x1*x2*x3+225831278*x1*x2*x4+32314614*x1*x3^2-11292836*x1*x3*x4-11471604*x1*x4^2-17361686*x2^3-61020343*x2^2*x3-132982128*x2^2*x4+46623254*x2*x3^2-53441750*x2*x3*x4-48142857*x2*x4^2+9644976*x3^3-19580172*x3^2*x4-28116390*x3*x4^2+18646000*x4^3+64756487*x1^2+115474646*x1*x2+150465505*x1*x3-38448834*x1*x4-91874169*x2^2-131012023*x2*x3-176642248*x2*x4+58052114*x3^2-141545442*x3*x4+9369657*x4^2-71771275*x1-58141098*x2-111938286*x3-11893932*x4+7665457];
-    I = Ideal(F)
+    I = AlgebraicSolving.Ideal(F)
     dec = equidimensional_decomposition(I)
     @test AlgebraicSolving._check_decomp(I, dec)
 end
diff --git a/test/algorithms/dimension.jl b/test/algorithms/dimension.jl
index dc836a4..7acc3a7 100644
--- a/test/algorithms/dimension.jl
+++ b/test/algorithms/dimension.jl
@@ -1,19 +1,19 @@
 @testset "Algorithms -> Dimension" begin
     R, (x,y) = polynomial_ring(QQ,["x","y"])
-    I =  Ideal([x^2,x*y])
+    I =  AlgebraicSolving.Ideal([x^2,x*y])
     @test isone(dimension(I))
     @test isone(I.dim)
 
     R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"])
-    I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
+    I = AlgebraicSolving.Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     @test iszero(dimension(I))
     @test iszero(I.dim)
 
-    I = Ideal([R(0)])
+    I = AlgebraicSolving.Ideal([R(0)])
     @test dimension(I) == ngens(R)
     @test I.dim == ngens(R)
 
-    I = Ideal([R(1)])
+    I = AlgebraicSolving.Ideal([R(1)])
     @test dimension(I) == -1
     @test I.dim == -1
 end
diff --git a/test/algorithms/groebner-bases.jl b/test/algorithms/groebner-bases.jl
index a436ab5..6e4c4db 100644
--- a/test/algorithms/groebner-bases.jl
+++ b/test/algorithms/groebner-bases.jl
@@ -1,6 +1,6 @@
 @testset "Algorithms -> Gröbner bases" begin
     R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"], internal_ordering=:degrevlex)
-    I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
+    I = AlgebraicSolving.Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     G = groebner_basis(I)
     H = MPolyRingElem[
          x + 2*y + 2*z + 100
@@ -17,7 +17,7 @@
         ]
     @test G == H
 
-    I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
+    I = AlgebraicSolving.Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     G = groebner_basis(I, eliminate=2, intersect=false)
     H = MPolyRingElem[
         z^4 + 38*z^3 + 95*z^2 + 95*z
@@ -31,19 +31,19 @@
     @test L == H
     @test I.gb[2] == H
 
-    I = Ideal([R(0)])
+    I = AlgebraicSolving.Ideal([R(0)])
     G = groebner_basis(I)
     @test G == [R(0)]
 
     R, (x1, x2) = polynomial_ring(QQ, ["x1", "x2"])
-    I = Ideal([3*x1^2 + ZZRingElem(2)^100, 2*x1*x2 + 5*x1 + ZZRingElem(2)^100])
+    I = AlgebraicSolving.Ideal([3*x1^2 + ZZRingElem(2)^100, 2*x1*x2 + 5*x1 + ZZRingElem(2)^100])
     G = groebner_basis(I)
     H = MPolyRingElem[
         3*x1 - 2*x2 - 5
         4*x2^2 + 20*x2 + 3802951800684688204490109616153
         ]
     @test G == H
-    J = Ideal([3*x1^2 + ZZRingElem(2)^100, 2*x1*x2 + 5*x1 + ZZRingElem(2)^100])
+    J = AlgebraicSolving.Ideal([3*x1^2 + ZZRingElem(2)^100, 2*x1*x2 + 5*x1 + ZZRingElem(2)^100])
     G = groebner_basis(J, normalize=true)
     H = MPolyRingElem[
         x1 - 2//3*x2 - 5//3
@@ -51,14 +51,14 @@
         ]
     @test G == H
     R, (x,y,z) = polynomial_ring(QQ,["x","y","z"], internal_ordering=:degrevlex)
-    I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
+    I = AlgebraicSolving.Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     G = eliminate(I, 2)
     H = MPolyRingElem[
         84*z^4 - 40*z^3 + z^2 + z
     ]
     @test G == H
     R, (a,b,c,d,x,y,z,w) = polynomial_ring(QQ, ["a", "b", "c", "d", "x", "y", "z", "w"])
-    I = Ideal([x - a*c, y - a*c*d, z - a*c^2 - b, w - a*c^2*d - b*d]);
+    I = AlgebraicSolving.Ideal([x - a*c, y - a*c*d, z - a*c^2 - b, w - a*c^2*d - b*d]);
     G = eliminate(I, 4)
     H = MPolyRingElem[
         -x*w + y*z
diff --git a/test/algorithms/hilbert.jl b/test/algorithms/hilbert.jl
index 80fedb5..790885f 100644
--- a/test/algorithms/hilbert.jl
+++ b/test/algorithms/hilbert.jl
@@ -1,6 +1,6 @@
 @testset "Algorithms -> Hilbert" begin
     R, (x,y) = polynomial_ring(QQ,["x","y"])
-    I =  Ideal([x^2,x*y])
+    I =  AlgebraicSolving.Ideal([x^2,x*y])
     A, t = polynomial_ring(ZZ, :t)
     B, s = polynomial_ring(QQ, :s)
     HS = (t^2 - t - 1)//(t - 1)
@@ -12,7 +12,7 @@
     @test isone(hilbert_degree(I))
 
     R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"])
-    I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
+    I = AlgebraicSolving.Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     HS = t^2 + 2*t + 1
     HP = (s^2 + 2*s + 1, 0)
 
@@ -21,7 +21,7 @@
     @test iszero(hilbert_dimension(I))
     @test 4 == hilbert_degree(I)
 
-    I = Ideal([R(0)])
+    I = AlgebraicSolving.Ideal([R(0)])
     HS = 1//(1-t)^(nvars(R))
     HP = (1//6*s^3 + s^2 + 11//6*s + 1, 0)
 
@@ -30,9 +30,9 @@
     @test nvars(R) == hilbert_dimension(I)
     @test isone(hilbert_degree(I))
 
-    I = Ideal([R(1)])
+    I = AlgebraicSolving.Ideal([R(1)])
     @test iszero(hilbert_series(I))
     @test all(iszero, hilbert_polynomial(I))
     @test hilbert_dimension(I) == -1
     @test iszero(hilbert_degree(I))
-end
+end
\ No newline at end of file
diff --git a/test/algorithms/monomial_diagram.jl b/test/algorithms/monomial_diagram.jl
new file mode 100644
index 0000000..3d00934
--- /dev/null
+++ b/test/algorithms/monomial_diagram.jl
@@ -0,0 +1,38 @@
+using Nemo
+using Random
+using StaticArrays
+
+@testset "Algorithms -> Monomial Diagram" begin
+    R, (x1,x2,x3,x4,x5,x6) = polynomial_ring(GF(101), ["x1", "x2", "x3", "x4", "x5", "x6"], internal_ordering=:degrevlex)
+    list_of_polynomials1 = [x1+2*x2+2*x3+2*x4+2*x5+2*x6-1, 
+    x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2-x1,
+    2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6-x2,
+    x2^2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6-x3,
+    2*x2*x3+2*x1*x4+2*x2*x5+2*x3*x6-x4,
+    x3^2+2*x2*x4+2*x1*x5+2*x2*x6-x5]
+    I = AlgebraicSolving.Ideal(list_of_polynomials1)
+    G = AlgebraicSolving.groebner_basis(I)
+    list_of_monomials1 = [AlgebraicSolving.convert_to_monomial(Nemo.leading_monomial(g), R, Val(6)) for g in G]
+    hashstate = AlgebraicSolving.new_hashstate()
+    diagram1 = AlgebraicSolving.create_diagram(list_of_monomials1, hashstate)
+    @test AlgebraicSolving.size_of_diagram(diagram1) == 155 
+    @test AlgebraicSolving.number_of_distinct_nodes(diagram1) == 21
+
+    hashstate = AlgebraicSolving.new_hashstate()
+    lst = AlgebraicSolving.generate_random_vectors(10, 5, 10)
+    
+    d1 = AlgebraicSolving.create_diagram(lst, hashstate)
+    for k in 1:30
+        lst2 = shuffle(lst)
+        d2 = AlgebraicSolving.create_diagram(lst2, hashstate)
+        @test d1 == d2
+    end
+
+    for k in 1:30
+        mon = AlgebraicSolving.generate_random_vectors(10, 5, 1)[1]
+        @test AlgebraicSolving.naive_is_in_ideal(mon, lst) == AlgebraicSolving.is_in_diagram(mon, d1)
+    end
+
+    leaves = AlgebraicSolving.get_leaves(d1)
+    @test is_empty([leaf for leaf in leaves if !AlgebraicSolving.naive_is_in_ideal(AlgebraicSolving.monomial(AlgebraicSolving.SVector{10}(reverse(leaf))), lst)])
+end
\ No newline at end of file
diff --git a/test/algorithms/normal-forms.jl b/test/algorithms/normal-forms.jl
index 6ad4ce7..11659e4 100644
--- a/test/algorithms/normal-forms.jl
+++ b/test/algorithms/normal-forms.jl
@@ -1,18 +1,18 @@
 @testset "Algorithms -> Normal forms" begin
     R, (x,y) = polynomial_ring(GF(101),["x","y"])
-    I = Ideal([y*x+17-y, x+13*y])
+    I = AlgebraicSolving.Ideal([y*x+17-y, x+13*y])
     f = x^4+y
-    @test normal_form(f, I) == 100*y+77
+    @test AlgebraicSolving.normal_form(f, I) == 100*y+77
     F = [x+13*y, x*y-4]
-    @test normal_form(F, I) == [0, y + 80]
+    @test AlgebraicSolving.normal_form(F, I) == [0, y + 80]
     R, (x,y) = polynomial_ring(GF(536870923),["x","y"])
-    I = Ideal([y*x+17-y, x+13*y])
+    I = AlgebraicSolving.Ideal([y*x+17-y, x+13*y])
     F = [x+13*y, x*y-4]
-    @test normal_form(F, I) == [0, y + 536870902]
+    @test AlgebraicSolving.normal_form(F, I) == [0, y + 536870902]
     R, (x,y) = polynomial_ring(GF(65521),["x","y"])
-    I = Ideal([y*x+17-y, x+13*y])
+    I = AlgebraicSolving.Ideal([y*x+17-y, x+13*y])
     F = [x+13*y, x*y+16]
-    @test normal_form(F, I) == [0, y + 65520]
-    I = Ideal([R(0)])
-    @test normal_form(F, I) == F
+    @test AlgebraicSolving.normal_form(F, I) == [0, y + 65520]
+    I = AlgebraicSolving.Ideal([R(0)])
+    @test AlgebraicSolving.normal_form(F, I) == F
 end
diff --git a/test/algorithms/param-curves.jl b/test/algorithms/param-curves.jl
index dd03592..47d9894 100644
--- a/test/algorithms/param-curves.jl
+++ b/test/algorithms/param-curves.jl
@@ -1,6 +1,6 @@
 @testset "Algorithms -> Curve Parametrization" begin
     R, (x1,x2,x3) = polynomial_ring(QQ, ["x1","x2","x3"])
-    I = Ideal([x1+2*x2+2*x3-1, x1^2+2*x2^2+2*x3^2-x1])
+    I = AlgebraicSolving.Ideal([x1+2*x2+2*x3-1, x1^2+2*x2^2+2*x3^2-x1])
 
     C, (x,y) = polynomial_ring(QQ, ["x","y"])
     elim  = x^2 + 4//3*x*y - 1//3*x + y^2 - 1//3*y
@@ -20,15 +20,15 @@
 
     R, (x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"])
 
-    I = Ideal([x1^2-x2, x1*x3-x4, x2*x4-12, x4^3-x3^2])
+    I = AlgebraicSolving.Ideal([x1^2-x2, x1*x3-x4, x2*x4-12, x4^3-x3^2])
     rational_curve_parametrization(I)
     @test I.rat_param.vars == Symbol[]
     @test I.rat_param.elim == -one(C)
     @test I.dim == -1
 
-    I = Ideal([x1^2-x2, x1*x3, x2-12])
+    I = AlgebraicSolving.Ideal([x1^2-x2, x1*x3, x2-12])
 	@test_throws AssertionError rational_curve_parametrization(I)
-    @test_throws ErrorException rational_curve_parametrization(Ideal([x1^2-x2, x1*x3, x2-12]), check_gen=false)
+    @test_throws ErrorException rational_curve_parametrization(AlgebraicSolving.Ideal([x1^2-x2, x1*x3, x2-12]), check_gen=false)
 
     elim    = 5041//2500*x^2 - 71//25*x*y - 3834//625*x + y^2 + 108//25*y - 6492//625
     denom   = -71//25*x + 2*y + 108//25
@@ -49,21 +49,21 @@
     @test param.param[4] == p4
     @test I.dim == 1
 
-    I = Ideal([x1^2-x2, x1*x3])
+    I = AlgebraicSolving.Ideal([x1^2-x2, x1*x3])
 	@test_throws AssertionError rational_curve_parametrization(I)
     @test_throws AssertionError rational_curve_parametrization(I, use_lfs=true)
     @test_throws AssertionError rational_curve_parametrization(I, check_gen=false)
 
-    I = Ideal([R(0)])
+    I = AlgebraicSolving.Ideal([R(0)])
     @test_throws AssertionError rational_curve_parametrization(I)
 
-    I = Ideal([R(1)])
+    I = AlgebraicSolving.Ideal([R(1)])
     @test rational_curve_parametrization(I).vars == Symbol[]
 
     # Non-generic variables
     A,(x,t,u) = polynomial_ring(QQ, [:x,:u,:t])
-    I = Ideal([t^2+u, u-x])
+    I = AlgebraicSolving.Ideal([t^2+u, u-x])
 	@test_throws AssertionError rational_curve_parametrization(I)
-    I = Ideal([u^2+t, t-x])
+    I = AlgebraicSolving.Ideal([u^2+t, t-x])
     @test_throws AssertionError rational_curve_parametrization(I)
 end
diff --git a/test/algorithms/solvers.jl b/test/algorithms/solvers.jl
index 7bc1ca0..faebbaf 100644
--- a/test/algorithms/solvers.jl
+++ b/test/algorithms/solvers.jl
@@ -1,6 +1,6 @@
 @testset "Algorithms -> Solvers" begin
     R, (x1,x2,x3,x4) = polynomial_ring(QQ,["x1","x2","x3","x4"], internal_ordering=:degrevlex)
-    I = Ideal([x1 + 2*x2 + 2*x3 + 2*x4 - 1,
+    I = AlgebraicSolving.Ideal([x1 + 2*x2 + 2*x3 + 2*x4 - 1,
          x1^2 + 2*x2^2 + 2*x3^2 + 2*x4^2 - x1,
          2*x1*x2 + 2*x2*x3 + 2*x3*x4 - x2,
          x2^2 + 2*x1*x3 + 2*x2*x4 - x3])
@@ -53,12 +53,12 @@
     @test I.rat_param.param[2] == p2
     @test I.rat_param.param[3] == p3
 
-    I = Ideal([x1^2-x2, x1*x3-x4, x2*x4-12, x4^3-x3^2])
+    I = AlgebraicSolving.Ideal([x1^2-x2, x1*x3-x4, x2*x4-12, x4^3-x3^2])
     real_solutions(I)
     @test I.rat_param.vars == Symbol[]
     @test I.dim == -1
 
-    I = Ideal([x1^2-x2, x1*x3, x2-12])
+    I = AlgebraicSolving.Ideal([x1^2-x2, x1*x3, x2-12])
 	@test_throws ErrorException real_solutions(I)
     @test_throws ErrorException real_solutions(I, interval=true)
 	@test_throws ErrorException rational_solutions(I)
@@ -66,25 +66,25 @@
 
     R, (x, y) = polynomial_ring(QQ,["x","y"])
     # issue 54
-    I = Ideal([R(0)])
+    I = AlgebraicSolving.Ideal([R(0)])
 	@test_throws ErrorException real_solutions(I)
-    I = Ideal([x-1,y+2,R(0)])
+    I = AlgebraicSolving.Ideal([x-1,y+2,R(0)])
     @test issetequal(real_solutions(I), Vector{QQFieldElem}[[1, -2]])
 
     # check variable permutation
-    I = Ideal([x^2-1, y])
+    I = AlgebraicSolving.Ideal([x^2-1, y])
     sols = Vector{QQFieldElem}[
         [-1, 0],
         [1, 0]
             ]
     @test issetequal(sols, real_solutions(I))
 
-    I = Ideal([x^2-1, y^2])
+    I = AlgebraicSolving.Ideal([x^2-1, y^2])
     @test issetequal(sols, real_solutions(I))
 
     # issue 5741 from Oscar
     R, (x,y,z) = polynomial_ring(QQ, ["x","y","z"])
-    I = Ideal([(x-1)*x, y-2, z-3])
+    I = AlgebraicSolving.Ideal([(x-1)*x, y-2, z-3])
     rat_sols = Vector{QQFieldElem}[[1, 2, 3], [0, 2, 3]]
     @test issetequal(rat_sols, rational_solutions(I))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 97e0071..e6c77da 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,7 @@ using Test
 @testset verbose = true "AlgebraicSolving Tests" begin
 include("interfaces/nemo.jl")
 include("algorithms/groebner-bases.jl")
+include("algorithms/monomial_diagram.jl")
 include("algorithms/normal-forms.jl")
 include("algorithms/solvers.jl")
 include("algorithms/dimension.jl")