diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index f300d337f..2f69dcbd8 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -206,11 +206,11 @@ with ϕᵢⱼ | Pᵢ | Pⱼ:
 
 """
 function project!(sm::GeneralizedStabilizer, p::PauliOperator)
-    newstab, anticom_idx, res = project!(sm.stab, p)
+    newstab, anticom, res = project!(sm.stab, p)
     χ′ = expect(p, sm)
     n = nqubits(newstab)
     newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>χ′))
-    return newsm, anticom_idx, res
+    return newsm, anticom, res
 end
 
 """
@@ -225,21 +225,19 @@ function projectrand!(sm::GeneralizedStabilizer, p::PauliOperator)
 end
 
 function _proj₊(sm::GeneralizedStabilizer, p::PauliOperator)
-    newstab, anticom_idx, res = project!(sm.stab, p)
+    newstab, res = projectrand!(sm.stab, p)
     χ′ = expect(p, sm)
-    sm.stab = newstab
     n = nqubits(newstab)
     newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>χ′))
-    return newsm, anticom_idx, res
+    return newsm, res
 end
 
 function _proj₋(sm::GeneralizedStabilizer, p::PauliOperator)
-    newstab, anticom_idx, res = project!(sm.stab, -p)
+    newstab, res = projectrand!(sm.stab, -p)
     χ′ = expect(p, sm)
-    sm.stab = newstab
     n = nqubits(newstab)
     newsm = GeneralizedStabilizer(newstab, DefaultDict(0.0im, (falses(n),falses(n))=>χ′))
-    return newsm, anticom_idx, res
+    return newsm, res
 end
 
 Base.copy(sm::GeneralizedStabilizer) = GeneralizedStabilizer(copy(sm.stab),copy(sm.destabweights))
diff --git a/test/test_nonclifford_quantumoptics.jl b/test/test_nonclifford_quantumoptics.jl
index 12c3b4d82..a7ec30189 100644
--- a/test/test_nonclifford_quantumoptics.jl
+++ b/test/test_nonclifford_quantumoptics.jl
@@ -52,22 +52,53 @@ qo_tgate.data[2,2] = exp(im*pi/4)
 end
 
 @testset "project!" begin
-    checks = [(S"X",P"X"),(S"Y",P"Y"),(S"Z",P"Z"),(S"-X",P"-X"),(S"-Y",P"-Y"),(S"-Z",P"-Z")]
-    for (s, p) in checks
-        gs = GeneralizedStabilizer(s)
-        apply!(gs, p)
-        qo_state = Operator(gs)
-        project!(copy(gs), copy(p))[1]
-        qo_state_after_proj = Operator(gs)
-        qo_pauli = Operator(p)
-        qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
-        qo_proj2 = (identityoperator(qo_pauli) - qo_pauli)/2
-        qo_proj3 = (identityoperator(qo_pauli) + qo_pauli)/2
-        qo_proj4 = (identityoperator(qo_pauli) + qo_pauli)/2
-        result1  = qo_proj1*qo_state*qo_proj1'
-        result2  = qo_proj2*qo_state*qo_proj2'
-        result3  = qo_proj3*qo_state*qo_proj3'
-        result4  = qo_proj4*qo_state*qo_proj4'
-        @test qo_state_after_proj ≈ result4 || qo_state_after_proj ≈ result3 || qo_state_after_proj ≈ result2 || qo_state_after_proj ≈ result1
+    for n in 1:4
+        for repetition in 1:5
+            s = random_stabilizer(n)
+            p = random_pauli(n)
+            gs = GeneralizedStabilizer(s)
+            for i in 1:rand(1:3)
+                apply!(gs, embed(n, i, pcT))
+            end
+            qo_state = Operator(gs)
+            project!(copy(gs), copy(p))[1]
+            qo_state_after_proj = Operator(gs)
+            qo_pauli = Operator(gs)
+            qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
+            qo_proj2 = (identityoperator(qo_pauli) - qo_pauli)/2
+            qo_proj3 = (identityoperator(qo_pauli) + qo_pauli)/2
+            qo_proj4 = (identityoperator(qo_pauli) + qo_pauli)/2
+            result1  = qo_proj1*qo_state*qo_proj1'
+            result2  = qo_proj2*qo_state*qo_proj2'
+            result3  = qo_proj3*qo_state*qo_proj3'
+            result4  = qo_proj4*qo_state*qo_proj4'
+            @test qo_state_after_proj ≈ result4 || qo_state_after_proj ≈ result3 || qo_state_after_proj ≈ result2 || qo_state_after_proj ≈ result1
+        end
+    end
+end
+
+@testset "projectrand!" begin
+    for n in 1:4
+        for repetition in 1:5
+            s = random_stabilizer(n)
+            p = random_pauli(n)
+            gs = GeneralizedStabilizer(s)
+            for i in 1:rand(1:3)
+                apply!(gs, embed(n, i, pcT))
+            end
+            qo_state = Operator(gs)
+            projectrand!(copy(gs), copy(p))[1]
+            qo_state_after_proj = Operator(gs)
+            qo_pauli = Operator(gs)
+            qo_proj1 = (identityoperator(qo_pauli) - qo_pauli)/2
+            qo_proj2 = (identityoperator(qo_pauli) - qo_pauli)/2
+            qo_proj3 = (identityoperator(qo_pauli) + qo_pauli)/2
+            qo_proj4 = (identityoperator(qo_pauli) + qo_pauli)/2
+            result1  = qo_proj1*qo_state*qo_proj1'
+            result2  = qo_proj2*qo_state*qo_proj2'
+            result3  = qo_proj3*qo_state*qo_proj3'
+            result4  = qo_proj4*qo_state*qo_proj4'
+            @test qo_state_after_proj ≈ result4 || qo_state_after_proj ≈ result3 || qo_state_after_proj ≈ result2 || qo_state_after_proj ≈ result1
+        end
     end
 end