diff --git a/src/entanglement.jl b/src/entanglement.jl
index 1803e0eb5..0c63989de 100644
--- a/src/entanglement.jl
+++ b/src/entanglement.jl
@@ -135,15 +135,35 @@ canonicalize_clip!(ps::Base.AbstractVecOrTuple{PauliOperator}, args...; kwargs..
 """
 $TYPEDSIGNATURES
 
-Get the bigram of a tableau.
+The Bigram `B` of stabilizer endpoints represents the "span" of each stabilizer within a set of Pauli operators `𝒢 = {g₁,…,gₙ}`.
 
-It is the list of endpoints of a tableau in the clipped gauge.
+For each stabilizer `g`, the left endpoint `𝓁(g)` is defined as the minimum site `x` where `g` acts non-trivially, while the
+right endpoint `𝓇(g)` is the maximum site where `g` acts non-trivially. 
+
+The site `x` represent the position within the system, taking values from `{1,2,…,n}` where `n` is the number of qubits.
+
+The bigram set `B(𝒢)` encodes these endpoints as pairs:
+
+`B(𝒢) ≡ {(𝓁(g₁),𝓇(g₁)),…,(𝓁(gₙ),𝓇(gₙ))}`
 
 If `clip=true` (the default) the tableau is converted to the clipped gauge in-place before calculating the bigram.
 Otherwise, the clip gauge conversion is skipped (for cases where the input is already known to be in the correct gauge).
 
 Introduced in [nahum2017quantum](@cite), with a more detailed explanation of the algorithm in [li2019measurement](@cite) and [gullans2021quantum](@cite).
 
+```jldoctest
+julia> s = ghz(3)
++ XXX
++ ZZ_
++ _ZZ
+
+julia> bigram(s)
+3×2 Matrix{Int64}:
+ 1  3
+ 1  2
+ 2  3
+```
+
 See also: [`canonicalize_clip!`](@ref)
 """
 function bigram(state::AbstractStabilizer; clip::Bool=true)::Matrix{Int} # JET-XXX The ::Matrix{Int} should not be necessary, but they help with inference
@@ -174,16 +194,41 @@ the most performant one depending on the particular case.
 
 Currently implemented are the `:clip` (clipped gauge), `:graph` (graph state), and `:rref` (Gaussian elimination) algorithms.
 Benchmark your particular case to choose the best one.
+
+See Appendix C of [nahum2017quantum](@cite).
 """
 function entanglement_entropy end
 
 
 """
+$TYPEDSIGNATURES
+
 Get bipartite entanglement entropy of a contiguous subsystem by passing through the clipped gauge.
 
 If `clip=false` is set the canonicalization step is skipped, useful if the input state is already in the clipped gauge.
 
-See also: [`bigram`](@ref), [`canonicalize_clip!`](@ref)
+```jldoctest
+julia> using Graphs # hide
+
+julia> s = ghz(3)
++ XXX
++ ZZ_
++ _ZZ
+
+julia> entanglement_entropy(s, 1:3, Val(:clip))
+0
+
+julia> s = Stabilizer(Graph(ghz(4)))
++ XZZZ
++ ZX__
++ Z_X_
++ Z__X
+
+julia> entanglement_entropy(s, [1,4], Val(:graph))
+1
+```
+
+See also: [`bigram`](@ref), [`canonicalize_clip!`](@ref).
 """
 function entanglement_entropy(state::AbstractStabilizer, subsystem_range::UnitRange, algorithm::Val{:clip}; clip::Bool=true)
     # JET-XXX The ::Matrix{Int} should not be necessary, but they help with inference
@@ -196,6 +241,8 @@ end
 
 
 """
+$TYPEDSIGNATURES
+
 Get bipartite entanglement entropy by first converting the state to a graph and computing the rank of the adjacency matrix.
 
 Based on "Entanglement in graph states and its applications".
@@ -210,6 +257,8 @@ end
 
 
 """
+$TYPEDSIGNATURES
+
 Get bipartite entanglement entropy by converting to RREF form (i.e., partial trace form).
 
 The state will be partially canonicalized in an RREF form.
@@ -231,3 +280,59 @@ function entanglement_entropy(state::AbstractStabilizer, subsystem::AbstractVect
 end
 
 entanglement_entropy(state::MixedDestabilizer, subsystem::AbstractVector, a::Val{:rref}) = entanglement_entropy(state, subsystem, a; pure=nqubits(state)==rank(state))
+
+"""
+$TYPEDSIGNATURES
+
+The mutual information between subsystems `𝒶` and `𝒷` in a stabilizer state is given by `I(𝒶, 𝒷) = S𝒶 + S𝒷 - S𝒶𝒷`.
+
+```jldoctest
+julia> using QuantumClifford
+
+julia> using Graphs; using QuantumClifford: mutual_information # hide
+
+julia> mutual_information(ghz(3), 1:2, 3:4, Val(:clip))
+2
+
+julia> s = Stabilizer(Graph(ghz(4)))
++ XZZZ
++ ZX__
++ Z_X_
++ Z__X
+
+julia> mutual_information(s, [1,2], [3, 4], Val(:graph))
+2
+```
+
+See Eq. E6 of [li2019measurement](@cite). See also: [`entanglement_entropy`](@ref)
+"""
+function mutual_information(state::AbstractStabilizer, A::AbstractVector{<:Integer}, B::AbstractVector{<:Integer}, alg::Val{:clip})
+    if !isempty(intersect(A, B))
+        throw(ArgumentError("Ranges A and B must not overlap."))
+    end
+    union_AB = union(A, B)
+    min_AB = minimum(union_AB)
+    max_AB = maximum(union_AB)
+    if length(union_AB) != max_AB - min_AB + 1
+        throw(ArgumentError("For the :clip algorithm, the union of A and B must form a contiguous range."))
+    end
+    contiguous_union = min_AB:max_AB
+    S_A = entanglement_entropy(state, A, alg)
+    S_B = entanglement_entropy(state, B, alg)
+    S_AB = entanglement_entropy(state, contiguous_union, alg)
+    return S_A + S_B - S_AB
+end
+
+function mutual_information(state::AbstractStabilizer, A::AbstractVector{<:Integer}, B::AbstractVector{<:Integer}, alg::Val{T}) where T
+    if !isempty(intersect(A, B))
+        throw(ArgumentError("Ranges A and B must not overlap."))
+    end
+    S_A = entanglement_entropy(state, A, alg)
+    S_B = entanglement_entropy(state, B, alg)
+    S_AB = entanglement_entropy(state, union(A, B), alg)
+    return S_A + S_B - S_AB
+end
+
+# Default method: use the :rref algorithm
+mutual_information(state::AbstractStabilizer, A::AbstractVector{<:Integer}, B::AbstractVector{<:Integer}) = 
+    mutual_information(state, A, B, Val{:rref}())
diff --git a/test/test_entanglement.jl b/test/test_entanglement.jl
index fa50fb369..88d77d968 100644
--- a/test/test_entanglement.jl
+++ b/test/test_entanglement.jl
@@ -3,7 +3,8 @@
 
     test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
 
-    using QuantumClifford: stab_looks_good, destab_looks_good, mixed_stab_looks_good, mixed_destab_looks_good
+    using QuantumClifford
+    using QuantumClifford: stab_looks_good, destab_looks_good, mixed_stab_looks_good, mixed_destab_looks_good, mutual_information, entanglement_entropy
 
     @testset "Clipped gauge of stabilizer states" begin
         for n in test_sizes
@@ -50,4 +51,60 @@
         @test entanglement_entropy(copy(s), subsystem, Val(:graph))==2
         @test entanglement_entropy(copy(s), subsystem, Val(:rref))==2
     end
+
+    @testset "Mutual information for Clifford circuits with entanglement_entropy check" begin
+        using QuantumOpticsBase
+        import QuantumOpticsBase: entanglement_entropy
+
+        for n in [4, 5, 6, 7] # exclude larger test sizes to avoid out of memory error
+            s = random_stabilizer(n)
+            a_start = rand(1:max(1, n-2))
+            a_end = rand(a_start:min(n-1, a_start+2))
+            subsystem_rangeA = a_start:a_end
+            b_start = min(n, a_end + rand(1:2))
+            b_end = rand(b_start:n)
+            subsystem_rangeB = b_start:b_end
+
+            if !isempty(intersect(subsystem_rangeA, subsystem_rangeB))
+                @test_throws ArgumentError mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:clip))
+                @test_throws ArgumentError mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:rref))
+                @test_throws ArgumentError mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:graph))
+            else
+                # Now we need to check if the union is contiguous for :clip
+                union_AB = union(subsystem_rangeA, subsystem_rangeB)
+                min_AB = minimum(union_AB)
+                max_AB = maximum(union_AB)
+                is_contiguous = (length(union_AB) == (max_AB - min_AB + 1))
+                if is_contiguous
+                    mi_clip  = mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:clip))
+                    mi_rref  = mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:rref))
+                    mi_graph = mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:graph))
+                    @test mi_clip == mi_rref == mi_graph
+                    @test mi_clip ≥ 0
+                    ψ = Ket(s)
+                    ρ = dm(ψ)
+                    S_A = QuantumOpticsBase.entanglement_entropy(ρ, subsystem_rangeA, entropy_vn)
+                    S_B = QuantumOpticsBase.entanglement_entropy(ρ, subsystem_rangeB, entropy_vn)
+                    # If A ∪ B covers the full system (1:n), set S_AB = 0 to avoid an invalid full-system trace in entanglement_entropy
+                    S_AB = union_AB == (1:n) ? 0.0 : QuantumOpticsBase.entanglement_entropy(ρ, union_AB, entropy_vn)
+                    # For a pure state: I(A:B) = [S(A) + S(B) - S(A∪B)] / 2, and convert nats → bits by dividing by log(2).
+                    mi_indep = (S_A + S_B - S_AB) / (2 * log(2))
+                    @test isapprox(mi_clip, mi_indep; atol=1e-6)
+                else
+                    @test_throws ArgumentError mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:clip))
+                    mi_rref  = mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:rref))
+                    mi_graph = mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:graph))
+                    @test mi_rref == mi_graph # We can still test consistency between rref and graph
+                    @test mi_rref ≥ 0
+                end
+            end
+        end
+        @testset "Explicit non-contiguous ranges with :clip" begin
+            n = 6
+            s = random_stabilizer(n)
+            subsystem_rangeA = 1:2
+            subsystem_rangeB = 5:6
+            @test_throws ArgumentError mutual_information(copy(s), subsystem_rangeA, subsystem_rangeB, Val(:clip))
+        end
+    end
 end