diff --git a/CHANGELOG.md b/CHANGELOG.md
index a150277a..bf99894d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,8 +5,16 @@
 
 # News
 
-## v0.9.8 - 2024-07-24
-
+## v0.9.8 - 2024-08-03
+
+- New group-theoretical tools:
+    - `groupify` to generate full stabilizer group from generating set
+    - `minimal_generating_set` function to find the minimal generating set of a set of operators
+    - `pauligroup` to generate the full Pauli group of a certain number of qubits
+    - `normalizer` to generate all Paulis that commute with a set of Paulis
+    - `centralizer` to find a subset of a set of Paulis such that each element in the subset commutes with each element in the set
+    - `contractor` to find a subset of Paulis in a tableau that have an identity operator on a certain qubit
+    - `delete_columns` to remove the operators corresponding to a certain qubit from all Paulis in a Stabilizer
 - `PauliError` can now encode biased noise during Pauli frame simulation, i.e. one can simulate only X errors, or only Y errors, or only Z errors, or some weighted combination of these.
 
 ## v0.9.7 - 2024-07-23
diff --git a/Project.toml b/Project.toml
index 196ddf83..550c8385 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org> and QuantumSavory community members"]
-version = "0.9.7"
+version = "0.9.8"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 44fb2331..65771059 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -74,6 +74,8 @@ export
     single_z, single_x, single_y,
     # Graphs
     graphstate, graphstate!, graph_gatesequence, graph_gate,
+    # Group theory tools
+    groupify, minimal_generating_set, pauligroup, normalizer, centralizer, contractor, delete_columns,
     # Clipped Gauge
     canonicalize_clip!, bigram, entanglement_entropy,
     # mctrajectories
@@ -1202,6 +1204,7 @@ include("sumtypes.jl")
 include("precompiles.jl")
 include("ecc/ECC.jl")
 include("nonclifford.jl")
+include("grouptableaux.jl")
 include("plotting_extensions.jl")
 #
 include("gpu_adapters.jl")
diff --git a/src/grouptableaux.jl b/src/grouptableaux.jl
new file mode 100644
index 00000000..b0f76e1c
--- /dev/null
+++ b/src/grouptableaux.jl
@@ -0,0 +1,211 @@
+using Graphs
+using LinearAlgebra
+
+"""
+Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)).
+
+The returned object is exponentially long.
+
+```jldoctest
+julia> groupify(S"XZ ZX")
++ __
++ XZ
++ ZX
++ YY
+```
+"""
+function groupify(s::Stabilizer)
+    # Create a `Tableau` of 2ⁿ n-qubit identity Pauli operators(where n is the size of 
+    # `Stabilizer` s), then multiply each one by a different subset of the elements in s to 
+    # create all 2ⁿ unique elements in the group generated by s, then return the `Tableau`.
+    n = length(s)::Int
+    group = zero(Tableau, 2^n, nqubits(s))
+    for i in 0:2^n-1
+        for (digit_order, j) in enumerate(digits(i, base=2, pad=n))
+            if j == 1
+                group[i+1] *= s[digit_order]
+            end
+        end
+    end
+    return group 
+end
+
+
+"""
+For a not-necessarily-minimal generating set,
+return the minimal generating set.
+
+The input has to have only real phases.
+
+```jldoctest
+julia> minimal_generating_set(S"__ XZ ZX YY")
++ XZ
++ ZX
+```
+"""
+function minimal_generating_set(s::Stabilizer)
+    # Canonicalize `Stabilizer` s, then return a `Stabilizer` with all non-identity Pauli operators 
+    # in the result. If s consists of only identity operators, return the negative 
+    # identity operator if one is contained in s, and the positive identity operator otherwise.
+    s, _, r = canonicalize!(copy(s), ranks=true)
+    if r == 0
+        gs = zero(Stabilizer, 1, nqubits(s))
+        if 0x02 in phases(s)
+            gs[1] = -1 * gs[1]
+        end
+        return gs
+    else
+        return s[1:r, :]
+    end
+end
+
+"""
+Return the full Pauli group of a given length. Phases are ignored by default, 
+but can be included by setting `phases=true`.
+
+```jldoctest
+julia> pauligroup(1)
++ _
++ X
++ Z
++ Y
+
+julia> pauligroup(1, phases=true)
++ _
++ X
++ Z
++ Y
+- _
+- X
+- Z
+- Y
++i_
++iX
++iZ
++iY
+-i_
+-iX
+-iZ
+-iY
+```
+"""
+function pauligroup(n::Int; phases=false)
+    if phases
+        s = zero(Tableau, 4^(n + 1), n)
+        paulis = ((false, false), (true, false), (false, true), (true, true))
+        for (i, P) in enumerate(Iterators.product(Iterators.repeated(paulis, n)...))
+            for (j, p) in enumerate(P)
+                s[i, j] = p
+            end
+        end
+        for i in 1:4^n
+            s[i+4^n] = -1 * s[i]
+        end
+        for i in 4^n+1:2*4^n
+            s[i+4^n] = -1im * s[i]
+        end
+        for i in 2*4^n+1:3*4^n
+            s[i+4^n] = -1 * s[i]
+        end
+    end
+    if !phases
+        s = zero(Tableau, 4^n, n)
+        paulis = ((false, false), (true, false), (false, true), (true, true))
+        for (i, P) in enumerate(Iterators.product(Iterators.repeated(paulis, n)...))
+            for (j, p) in enumerate(P)
+                s[i, j] = p
+            end
+        end
+    end
+    return s
+end
+
+"""
+Return all Pauli operators with the same number of qubits as the given `Tableau` `t`
+that commute with all operators in `t`.
+
+```jldoctest
+julia> normalizer(T"X")
++ _
++ X
+```
+"""
+function normalizer(t::Tableau)
+    # For each `PauliOperator` p in the with same number of qubits as the `Stabilizer` s, iterate through s and check each 
+    # operator's commutivity with p. If they all commute, add p a vector of `PauliOperators`. Return the vector 
+    # converted to `Tableau`.
+    n = nqubits(t)
+    pgroup = pauligroup(n, phases=false)
+    ptype = typeof(t[1])
+    normalizer = ptype[]
+    for p in pgroup
+        commutes = true
+        for q in t
+            if comm(p, q) == 0x01
+                commutes = false
+            end
+        end
+        if commutes
+            push!(normalizer, p)
+        end
+    end
+
+    return Tableau(normalizer)
+end
+
+"""
+For a given set of Paulis (in the form of a `Tableau`), return the subset of Paulis that commute with all Paulis in set.
+
+```jldoctest
+julia> centralizer(T"XX ZZ _Z")
++ ZZ
+```
+"""
+function centralizer(t::Tableau) 
+    center = typeof(t[1])[]
+    for P in t
+        commutes = 0
+        for Q in t
+            if comm(P, Q) == 0x01
+                commutes = 1
+                break
+            end
+        end
+        if commutes == 0
+            push!(center, P)
+        end
+    end
+    if length(center) == 0 
+        return Tableau(zeros(Bool, 1,1))
+    end
+    c = Tableau(center)
+    return c
+end
+
+"""
+Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to 
+the given subset, without the entries corresponding to subset.
+
+```jldoctest
+julia> contractor(S"_X X_", [1])
++ X
+```
+"""
+function contractor(s::Stabilizer, subset)
+    result = typeof(s[1])[]
+    for p in s
+        contractable = true
+        for i in subset
+            if p[i] != (false, false) 
+                contractable = false 
+                break 
+            end
+        end
+        if contractable push!(result, p[setdiff(1:length(p), subset)]) end
+    end
+    if length(result) > 0
+        return Tableau(result)
+    else
+        return Tableau(zeros(Bool, 1,1))
+    end
+end
\ No newline at end of file
diff --git a/src/pauli_operator.jl b/src/pauli_operator.jl
index 640092b5..f607ab07 100644
--- a/src/pauli_operator.jl
+++ b/src/pauli_operator.jl
@@ -122,6 +122,11 @@ Base.hash(p::PauliOperator, h::UInt) = hash(p.phase,hash(p.nqubits,hash(p.xz, h)
 
 Base.copy(p::PauliOperator) = PauliOperator(copy(p.phase),p.nqubits,copy(p.xz))
 
+function Base.deleteat!(p::PauliOperator, subset) 
+    p =p[setdiff(1:length(p), subset)]
+    return p
+end
+
 _nchunks(i::Int,T::Type{<:Unsigned}) = 2*( (i-1) ÷ (8*sizeof(T)) + 1 )
 Base.zero(::Type{PauliOperator{Tₚ, Tᵥ}}, q) where {Tₚ,T<:Unsigned,Tᵥ<:AbstractVector{T}} = PauliOperator(zeros(UInt8), q, zeros(T, _nchunks(q,T)))
 Base.zero(::Type{PauliOperator}, q) = zero(PauliOperator{Array{UInt8, 0}, Vector{UInt}}, q)
diff --git a/src/project_trace_reset.jl b/src/project_trace_reset.jl
index 8ea01f38..19e49521 100644
--- a/src/project_trace_reset.jl
+++ b/src/project_trace_reset.jl
@@ -627,6 +627,8 @@ end
 $TYPEDSIGNATURES
 
 Trace out a qubit.
+
+See also: [`delete_columns`](@ref)
 """ # TODO all of these should raise an error if length(qubits)>rank
 function traceout!(s::Stabilizer, qubits; phases=true, rank=false)
     _,i = canonicalize_rref!(s,qubits;phases=phases)
@@ -866,3 +868,22 @@ function traceoutremove!(s::MixedDestabilizer, qubit)
     traceout!(s,[qubit]) # TODO this can be optimized thanks to the information already known from projfunc
     s = _remove_rowcol!(s, nqubits(s), qubit)
 end
+
+
+"""
+Return the given stabilizer without all the qubits in the given iterable.
+
+The resulting tableaux is not guaranteed to be valid (to retain its commutation relationships).
+
+```jldoctest
+julia> delete_columns(S"XYZ YZX ZXY", [1,3])
++ Y
++ Z
++ X
+```
+
+See also: [`traceout!`](@ref)
+"""
+function delete_columns(𝒮::Stabilizer, subset)
+    return 𝒮[:, setdiff(1:nqubits(𝒮), subset)]
+end
\ No newline at end of file
diff --git a/test/test_group_tableaux.jl b/test/test_group_tableaux.jl
new file mode 100644
index 00000000..afec47e9
--- /dev/null
+++ b/test/test_group_tableaux.jl
@@ -0,0 +1,126 @@
+@testitem "Classical" begin
+    using Test
+
+    using Random
+    using QuantumClifford
+
+    # Including sizes that would test off-by-one errors in the bit encoding.
+    test_sizes = [1, 2, 3, 4, 5, 7, 8, 9, 15, 16, 17] 
+    # Zero function(in groupify) slows down around 2^30(n=30),eventually breaks
+    small_test_sizes = [1, 2, 3, 4, 5, 7] # Pauligroup slows around n = 8
+
+    @testset "group_tableaux" begin
+        #Test groupify
+        for n in [1, test_sizes...]
+            s = random_stabilizer(n)
+            s_test = copy(s)
+            group = groupify(s)
+            @test length(group) == 2^n
+            unchanged = true
+            for stabilizer in group
+                apply!(s, stabilizer)
+                if !(s == s_test)
+                    unchanged = false
+                end
+                @test unchanged == true
+            end
+        end
+        #Test minimal_generating_set
+        for n in [1, small_test_sizes...]
+            s = random_stabilizer(n)
+            group = groupify(s)
+            gen_set = minimal_generating_set(Stabilizer(group))
+            new_group = groupify(gen_set)
+            canonicalize!(Stabilizer(group))
+            canonicalize!(Stabilizer(new_group))
+            @test group == new_group
+            s = zero(Stabilizer, rand(1:(2*n)), n)
+            for i in 1:length(s)
+                s[i] = random_pauli(n)
+            end
+            gen_set = minimal_generating_set(s)
+            new_group = groupify(s)
+            for operator in s
+                @test operator in new_group
+            end
+        end
+        #Test pauligroup
+        for n in [1, small_test_sizes...] 
+            @test length(QuantumClifford.pauligroup(n, phases=false)) == 4^n
+            @test length(QuantumClifford.pauligroup(n, phases=true)) == 4^(n+1)
+        end
+        #Test normalizer
+        for n in [1, small_test_sizes...] # pauligroup is very slow at n=14
+            t = zero(QuantumClifford.Tableau, rand(1:(2*n)), n)
+            for i in eachindex(t) t[i] = random_pauli(n) end
+            normalized = normalizer(t)
+            paulis = QuantumClifford.pauligroup(n, phases=false)
+            for n_pauli in normalized
+                for t_pauli in t
+                    @test comm(n_pauli, t_pauli) == 0x0
+                end
+            end
+            for pauli in paulis
+                commutes = true
+                for t_pauli in t
+                    if comm(t_pauli, pauli) == 0x01
+                        commutes = false
+                    end
+                end
+            @test (!commutes) || (pauli in normalized)
+            end
+        end
+        #Test centralizer
+        for n in [1, test_sizes...]
+            t = zero(QuantumClifford.Tableau, rand(1:(2*n)), n)
+            for i in eachindex(t) t[i] = random_pauli(n) end
+            c = centralizer(t)
+            for c_pauli in c
+                for t_pauli in t
+                    @test comm(c_pauli, t_pauli) == 0x0
+                end
+            end
+            for pauli in t
+                commutes = true
+                for t_pauli in t
+                    if comm(t_pauli, pauli)==0x01 commutes = false end
+                end
+                @test !commutes || pauli in c
+            end
+
+        end
+        #Test contractor
+        for n in [1, test_sizes...]
+            s = random_stabilizer(n)
+            subset = []
+            for i in 1:nqubits(s) #create a random subset
+                if rand(1:2) == 1 push!(subset, i) end
+            end
+            c = contractor(s, subset)
+            count = 0
+            for stabilizer in s 
+                contractable = true
+                for i in subset
+                    if stabilizer[i] != (false, false) contractable = false end
+                end
+                if contractable count+=1 end
+            end
+            if length(c[1]) > 0
+            @test count == length(c)
+                for contracted in c
+                    p = zero(PauliOperator, nqubits(s))
+                    index = 0
+                    for i in 1:nqubits(s)
+                        if !(i in subset) 
+                            index+=1 
+                            p[i] = contracted[index] 
+                        end
+                    end
+                    @test p in s || -1* p in s || 1im * p in s || -1im * p in s
+                end
+            else
+                @test count ==0
+            end
+        end
+    end
+end
\ No newline at end of file