diff --git a/docs/src/hilbert.md b/docs/src/hilbert.md
new file mode 100644
index 0000000..73f3e73
--- /dev/null
+++ b/docs/src/hilbert.md
@@ -0,0 +1,40 @@
+```@meta
+CurrentModule = AlgebraicSolving
+DocTestSetup = quote
+  using AlgebraicSolving
+end
+```
+
+```@setup algebraicsolving
+using AlgebraicSolving
+```
+
+```@contents
+Pages = ["hilbert.md"]
+```
+
+# Hilbert series of an ideal
+
+## Introduction
+
+AlgebraicSolving allows to compute the Hilbert series for the ideal spanned
+by given input generators over finite fields of characteristic smaller
+$2^{31}$ and over the rationals.
+
+The underlying engine is provided by msolve.
+
+## Functionality
+
+```@docs
+    hilbert_series(I::Ideal{T}) where T <: MPolyRingElem
+```
+
+In addition, from the same input, AlgebraicSolving can also compute the dimension and degree of the ideal, as well as the Hilbert polynomial and index of regularity.
+
+```@docs
+    hilbert_dimension(I::Ideal{T}) where T <: MPolyRingElem
+
+    hilbert_degree(I::Ideal{T}) where T <: MPolyRingElem
+
+    hilbert_polynomial(I::Ideal{T}) where T <: MPolyRingElem
+```
\ No newline at end of file
diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl
index 1a2fbb6..57bdcf7 100644
--- a/src/AlgebraicSolving.jl
+++ b/src/AlgebraicSolving.jl
@@ -14,6 +14,7 @@ include("algorithms/normal-forms.jl")
 include("algorithms/solvers.jl")
 include("algorithms/dimension.jl")
 include("algorithms/decomposition.jl")
+include("algorithms/hilbert.jl")
 #= siggb =#
 include("siggb/siggb.jl")
 #= examples =#
diff --git a/src/algorithms/dimension.jl b/src/algorithms/dimension.jl
index cf9fa45..89ade19 100644
--- a/src/algorithms/dimension.jl
+++ b/src/algorithms/dimension.jl
@@ -3,7 +3,7 @@
 
 Compute the Krull dimension of a given polynomial ideal `I`.
 
-**Note**: This requires a Gröbner basis of `I`, which is computed internally if not alraedy known.
+**Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
 
 # Examples
 ```jldoctest
@@ -26,13 +26,16 @@ function dimension(I::Ideal{T}) where T <: MPolyRingElem
     R = parent(first(gb))
 
     res = Set([trues(ngens(R))])
-    lead_exps = (_drl_lead_exp).(gb)
+    lead_exps = Vector{Vector{Int}}(undef, length(gb))
+    for i in eachindex(gb)
+        lead_exps[i] = _lead_exp_ord(gb[i], :degrevlex)
+    end
     for lexp in lead_exps
         nz_exps = (!iszero).(lexp)
         nz_exps_ind = findall(nz_exps)
         next_res = Set{BitVector}()
         for mis in res
-            if all_lesseq(nz_exps, mis)
+            if _all_lesseq(nz_exps, mis)
                 @inbounds for j in nz_exps_ind
                     new_mis = copy(mis)
                     new_mis[j] = false
@@ -49,7 +52,7 @@ function dimension(I::Ideal{T}) where T <: MPolyRingElem
     return I.dim
 end
 
-function all_lesseq(a::BitVector, b::BitVector)::Bool
+function _all_lesseq(a::BitVector, b::BitVector)::Bool
     @inbounds for i in eachindex(a)
         if a[i] && !b[i]
             return false
@@ -58,12 +61,15 @@ function all_lesseq(a::BitVector, b::BitVector)::Bool
     return true
 end
 
-function _drl_exp_vector(u::Vector{Int})
-    return [sum(u), -reverse(u)...]
-end
+function _lead_exp_ord(p::MPolyRingElem, order::Symbol)
+    R = parent(p)
+    internal_ordering(R)==order && return first(exponent_vectors(p))
 
-function _drl_lead_exp(p::MPolyRingElem)
-    exps = collect(Nemo.exponent_vectors(p))
-    _, i = findmax(_drl_exp_vector.(exps))
-    return exps[i]
+    A = base_ring(R)
+    R1, _ = polynomial_ring(A, symbols(R), internal_ordering=order)
+    ctx = MPolyBuildCtx(R1)
+    for e in exponent_vectors(p)
+        push_term!(ctx, one(A), e)
+    end
+    return first(exponent_vectors(finish(ctx)))
 end
\ No newline at end of file
diff --git a/src/algorithms/hilbert.jl b/src/algorithms/hilbert.jl
new file mode 100644
index 0000000..2ff270c
--- /dev/null
+++ b/src/algorithms/hilbert.jl
@@ -0,0 +1,284 @@
+@doc Markdown.doc"""
+    hilbert_series(I::Ideal{T}) 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.
+
+**Notes**:
+* This requires a Gröbner basis of `I`, which is computed internally if not already known.
+* Significantly faster when internal_ordering is :degrevlex.
+
+# Examples
+```jldoctest
+julia> using AlgebraicSolving
+
+julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);
+
+julia> I = Ideal([x*y,x*z,y*z]);
+
+julia> hilbert_series(I)
+(-2*t - 1)//(t - 1)
+```
+"""
+function hilbert_series(I::Ideal{T}) where T <: MPolyRingElem
+
+    gb = get!(I.gb, 0) do
+        groebner_basis(I, complete_reduction = true)
+    end
+    lead_exps = Vector{Vector{Int}}(undef, length(gb))
+    for i in eachindex(gb)
+        lead_exps[i] = _lead_exp_ord(gb[i], :degrevlex)
+    end
+    return _hilbert_series_mono(lead_exps)
+end
+
+@doc Markdown.doc"""
+    hilbert_degree(I::Ideal{T}) where T <: MPolyRingElem
+
+Compute the degree of a given polynomial ideal `I` by first computing its Hilbert series.
+
+**Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
+
+# Examples
+```jldoctest
+julia> using AlgebraicSolving
+
+julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);
+
+julia> I = Ideal([x*y,x*z,y*z]);
+
+julia> hilbert_degree(I)
+3
+```
+"""
+function hilbert_degree(I::Ideal{T}) where T <: MPolyRingElem
+
+    !isnothing(I.deg) && return I.deg
+    I.deg = numerator(hilbert_series(I))(1) |> abs
+    return I.deg
+end
+
+@doc Markdown.doc"""
+    hilbert_dimension(I::Ideal{T}) where T <: MPolyRingElem
+
+Compute the Krull dimension of a given polynomial ideal `I` by first computing its Hilbert series.
+
+**Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
+
+# Examples
+```jldoctest
+julia> using AlgebraicSolving
+
+julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);
+
+julia> I = Ideal([x*y,x*z,y*z]);
+
+julia> hilbert_dimension(I)
+1
+```
+"""
+function hilbert_dimension(I::Ideal{T}) where T <: MPolyRingElem
+
+    H = hilbert_series(I)
+    I.dim = iszero(H) ? -1 : degree(denominator(H))
+    return I.dim
+end
+
+@doc Markdown.doc"""
+    hilbert_polynomial(I::Ideal{T}) 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
+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.
+
+**Note**: This requires a Gröbner basis of `I`, which is computed internally if not already known.
+
+# Examples
+```jldoctest
+julia> using AlgebraicSolving
+
+julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);
+
+julia> I = Ideal([x*y,x*z,y*z]);
+
+julia> hilbert_polynomial(I)
+(3*s + 3, 1)
+```
+"""
+function hilbert_polynomial(I::Ideal{T}) where T <: MPolyRingElem
+
+    A, s = polynomial_ring(QQ, :s)
+    H = hilbert_series(I)
+    dim = degree(denominator(H))
+    num = iseven(dim) ? numerator(H) : -numerator(H)
+    dim==0 && return num(s), 0
+
+    t = gen(parent(num))
+    La = Vector{ZZPolyRingElem}(undef, dim)
+    while dim>0
+        num, La[dim] = divrem(num, 1-t)
+        dim -= 1
+    end
+
+    Hpolyfct = d->sum(La[i](0)*binomial(i+d, i) for i in 1:length(La))
+    dim = degree(denominator(H))
+    Hpoly = interpolate(A, QQ.(0:dim+1), [QQ(Hpolyfct(d)) for d in 0:dim+1])
+    @assert(degree(Hpoly)==dim, "Degree of poly does not match the dimension")
+    # Hilbert poly, index of regularity
+    return Hpoly, degree(num)+1
+end
+
+# Computes hilbert series of a monomial ideal on input list of exponents
+function _hilbert_series_mono(exps::Vector{Vector{Int}})
+
+    h = _num_hilbert_series_mono(exps)
+    t = gen(parent(h))
+    return h//(1-t)^length(first(exps))
+end
+
+# Computes numerator hilbert series of a monomial ideal on input list of exponents
+function _num_hilbert_series_mono(exps::Vector{Vector{Int}})
+
+    A, t = polynomial_ring(ZZ, 't')
+    r = length(exps)
+    r == 0 && return one(A)
+    N = length(first(exps))
+    ## Base cases ##
+    r == 1 && return (1-t^sum(first(exps)))
+    supp = findall.(Ref(!iszero), exps)
+    pow_supp = findall(s->length(s)==1, supp)
+    # If exps is a product of simple powers
+    if length(pow_supp) == r
+        #println("Simple power")
+        return prod(1-t^(exps[i][supp[i][1]]) for i in pow_supp)
+    # Only one non-simple power P
+    elseif length(pow_supp) == r-1
+        #println("Mixed pow")
+        inpow = setdiff(eachindex(exps), pow_supp) |> first
+        # P has disjoint support with other powers
+        if all(iszero(exps[inpow][supp[i][1]]) for i in pow_supp)
+            return (1-t^sum(exps[inpow]))*prod(1-t^(exps[i][supp[i][1]]) for i in pow_supp)
+        else
+            return prod(1-t^(exps[i][supp[i][1]]) for i in pow_supp) - t^sum(exps[inpow]) *
+            prod(1-t^(exps[i][supp[i][1]]-exps[inpow][supp[i][1]]) for i in pow_supp)
+        end
+    end
+
+    # Variable index occuring the most in exps
+    counts = sum(x->x .> 0, eachcol(reduce(hcat, exps)))
+    ivarmax = argmax(counts)
+
+    ## Splitting recursive cases ##
+    # Monomials have disjoint supports
+    if counts[ivarmax] == 1
+        return prod(1-t^sum(mono) for mono in exps)
+    # Heuristic where general splitting is useful
+    elseif 8 <= r <= N
+        # Finest partition of monomial supports
+        LV, h = _monomial_support_partition(exps), one(A)
+        rem_mon = collect(1:r)
+        # If we are indeed splitting
+        if length(LV) > 1
+            for V in LV
+                JV, iJV = Vector{Vector{Int}}(), Int[]
+                for (k, i) in enumerate(rem_mon)
+                    mono = exps[i]
+                    if any(mono[j] != 0 for j in V)
+                        push!(iJV, k)
+                        push!(JV, mono)
+                    end
+                end
+                h *= _num_hilbert_series_mono(JV)
+                # Avoid re-check monomials
+                deleteat!(rem_mon, iJV)
+            end
+            return h
+        end
+    end
+
+    ## Pivot recursive case ##
+    # Exponent of ivarmax in gcd of two random generators
+    pivexp = max(1, minimum(mon[ivarmax] for mon in rand(exps, 2)))
+    h = zero(A)
+    #Compute and partition gens of (exps):pivot
+    Lquo = [Vector{Int}[] for _ in 1:pivexp+2]
+    trivialquo = false
+    for mono in exps
+        if mono[ivarmax] <= pivexp
+            monoquo = vcat(mono[1:ivarmax-1], 0, mono[ivarmax+1:end])
+            if iszero(monoquo)
+                trivialquo = true
+                break
+            end
+            push!(Lquo[mono[ivarmax]+1],  monoquo)
+        else
+            push!(Lquo[pivexp+2],
+                vcat(mono[1:ivarmax-1], mono[ivarmax]-pivexp, mono[ivarmax+1:end]))
+        end
+    end
+    if !trivialquo
+        # Interreduce generators based on partition
+        @inbounds for i in pivexp+1:-1:1
+            non_min = [ k for (k,mono) in enumerate(Lquo[i]) if
+                    any(_all_lesseq(mini, mono) for j in i+1:pivexp+1 for mini in Lquo[j])]
+            deleteat!(Lquo[i], non_min)
+        end
+        # Merge all partitions
+        h += _num_hilbert_series_mono(vcat(Lquo...))*t^pivexp
+    end
+    # Interreduce (exps) + pivot
+    filter!(e->(pivexp > e[ivarmax]), exps)
+    push!(exps,[zeros(Int64,ivarmax-1); pivexp; zeros(Int64,N-ivarmax)])
+    h += _num_hilbert_series_mono(exps)
+
+    return h
+end
+
+function _all_lesseq(a::Vector{Int}, b::Vector{Int})::Bool
+    @inbounds for i in eachindex(a)
+        if a[i] > b[i]
+            return false
+        end
+    end
+    return true
+end
+
+# Build adjacency graph: connect variables that appear together in a monomial
+function _monomial_support_partition(L::Vector{Vector{Int}})
+
+    n = length(first(L))
+    adj = [Set{Int}() for _ in 1:n]
+    active = falses(n)
+
+    for mono in L
+        support = findall(!=(0), mono)
+        foreach(i -> active[i] = true, support)
+        for i in support, j in support
+            i != j && push!(adj[i], j)
+        end
+    end
+
+    visited = falses(n)
+    components = Vector{Vector{Int}}()
+
+    function dfs(u, comp)
+        visited[u] = true
+        push!(comp, u)
+        foreach(v -> !visited[v] && dfs(v, comp), adj[u])
+    end
+
+    for v in 1:n
+        if active[v] && !visited[v]
+            comp = Int[]
+            dfs(v, comp)
+            push!(components, comp)
+        end
+    end
+
+    return components
+end
diff --git a/src/exports.jl b/src/exports.jl
index 99c51de..5f7200d 100644
--- a/src/exports.jl
+++ b/src/exports.jl
@@ -3,4 +3,5 @@ export polynomial_ring, MPolyRing, GFElem, MPolyRingElem,
        QQ, vars, nvars, ngens, ZZRingElem, QQFieldElem,
        QQMPolyRingElem, base_ring, coefficient_ring, evaluate,
        prime_field, sig_groebner_basis, cyclic, leading_coefficient,
-       equidimensional_decomposition, homogenize, dimension, FqMPolyRingElem
+       equidimensional_decomposition, homogenize, dimension, FqMPolyRingElem,
+       hilbert_series, hilbert_dimension, hilbert_degree, hilbert_polynomial
diff --git a/src/types.jl b/src/types.jl
index abd7567..f32d6d1 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -28,6 +28,7 @@ end
 mutable struct Ideal{T <: MPolyRingElem}
     gens::Vector{T}
     dim::Union{Int, Nothing}
+    deg::Union{Int, Nothing}
     gb::Dict{Int, Vector{T}}
     inter_sols::Vector{Vector{Vector{QQFieldElem}}}
     real_sols::Vector{Vector{QQFieldElem}}
@@ -39,6 +40,7 @@ mutable struct Ideal{T <: MPolyRingElem}
         I.gens = F
         I.gb   = Dict()
         I.dim  = nothing
+        I.deg  = nothing
         return I
     end
 end
diff --git a/test/algorithms/hilbert.jl b/test/algorithms/hilbert.jl
new file mode 100644
index 0000000..361e92e
--- /dev/null
+++ b/test/algorithms/hilbert.jl
@@ -0,0 +1,29 @@
+@testset "Algorithms -> Hilbert" begin
+    R, (x,y) = polynomial_ring(QQ,["x","y"])
+    I =  Ideal([x^2,x*y])
+    A, t = polynomial_ring(ZZ, :t)
+    B, s = polynomial_ring(QQ, :s)
+    HS = (t^2 - t - 1)//(t - 1)
+    HP = (s + 1, 2)
+
+    @test HS == hilbert_series(I)
+    @test HP == hilbert_polynomial(I)
+    @test isone(hilbert_dimension(I))
+    @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])
+    HS = t^2 + 2*t + 1
+    HP = (s^2 + 2*s + 1, 0)
+
+    @test HS == hilbert_series(I)
+    @test HP == hilbert_polynomial(I)
+    @test iszero(hilbert_dimension(I))
+    @test 4 == hilbert_degree(I)
+
+    I = 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
diff --git a/test/runtests.jl b/test/runtests.jl
index 8da5e27..87e9135 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ include("algorithms/groebner-bases.jl")
 include("algorithms/normal-forms.jl")
 include("algorithms/solvers.jl")
 include("algorithms/dimension.jl")
+include("algorithms/hilbert.jl")
 include("algorithms/decomposition.jl")
 include("examples/katsura.jl")
 end