From 6a8dbb426da67c05c3dc9f74bb67941307d938c3 Mon Sep 17 00:00:00 2001
From: Tommy Hofmann <thofma@gmail.com>
Date: Tue, 22 Oct 2024 04:30:19 +0200
Subject: [PATCH 01/15] ci: tweak ci.yml (#401)

---
 .github/workflows/ci.yml | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index e00b68835..1f7eba665 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -4,6 +4,14 @@ on:
     branches: [master, main]
     tags: ["*"]
   pull_request:
+
+concurrency:
+  # group by workflow and ref; the last slightly strange component ensures that for pull
+  # requests, we limit to 1 concurrent job, but for the master branch we don't
+  group: ${{ github.workflow }}-${{ github.ref }}-${{ github.ref != 'refs/heads/master' || github.run_number }}
+  # Cancel intermediate builds, but only if it is a pull request build.
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
 env:
   PYTHON: ~
 jobs:

From 73ccb34553070bbf5f31c589aa744352929ff1a0 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Thu, 24 Oct 2024 19:57:10 -0400
Subject: [PATCH 02/15] remove unnecessary imports

---
 src/grouptableaux.jl        | 33 +++++++++++++++------------------
 test/test_group_tableaux.jl | 17 ++++++++---------
 2 files changed, 23 insertions(+), 27 deletions(-)

diff --git a/src/grouptableaux.jl b/src/grouptableaux.jl
index b0f76e1c2..208f1ddfb 100644
--- a/src/grouptableaux.jl
+++ b/src/grouptableaux.jl
@@ -1,6 +1,3 @@
-using Graphs
-using LinearAlgebra
-
 """
 Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)).
 
@@ -15,8 +12,8 @@ julia> groupify(S"XZ ZX")
 ```
 """
 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 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))
@@ -27,7 +24,7 @@ function groupify(s::Stabilizer)
             end
         end
     end
-    return group 
+    return group
 end
 
 
@@ -44,8 +41,8 @@ julia> minimal_generating_set(S"__ XZ ZX YY")
 ```
 """
 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 
+    # 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
@@ -60,7 +57,7 @@ function minimal_generating_set(s::Stabilizer)
 end
 
 """
-Return the full Pauli group of a given length. Phases are ignored by default, 
+Return the full Pauli group of a given length. Phases are ignored by default,
 but can be included by setting `phases=true`.
 
 ```jldoctest
@@ -131,8 +128,8 @@ julia> normalizer(T"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 
+    # 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)
@@ -161,7 +158,7 @@ julia> centralizer(T"XX ZZ _Z")
 + ZZ
 ```
 """
-function centralizer(t::Tableau) 
+function centralizer(t::Tableau)
     center = typeof(t[1])[]
     for P in t
         commutes = 0
@@ -175,7 +172,7 @@ function centralizer(t::Tableau)
             push!(center, P)
         end
     end
-    if length(center) == 0 
+    if length(center) == 0
         return Tableau(zeros(Bool, 1,1))
     end
     c = Tableau(center)
@@ -183,7 +180,7 @@ function centralizer(t::Tableau)
 end
 
 """
-Return the subset of Paulis in a Stabilizer that have identity operators on all qubits corresponding to 
+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
@@ -196,9 +193,9 @@ function contractor(s::Stabilizer, subset)
     for p in s
         contractable = true
         for i in subset
-            if p[i] != (false, false) 
-                contractable = false 
-                break 
+            if p[i] != (false, false)
+                contractable = false
+                break
             end
         end
         if contractable push!(result, p[setdiff(1:length(p), subset)]) end
@@ -208,4 +205,4 @@ function contractor(s::Stabilizer, subset)
     else
         return Tableau(zeros(Bool, 1,1))
     end
-end
\ No newline at end of file
+end
diff --git a/test/test_group_tableaux.jl b/test/test_group_tableaux.jl
index afec47e90..d208afdff 100644
--- a/test/test_group_tableaux.jl
+++ b/test/test_group_tableaux.jl
@@ -1,11 +1,10 @@
-@testitem "Classical" begin
+@testitem "group theory routines" 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] 
+    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
 
@@ -45,7 +44,7 @@
             end
         end
         #Test pauligroup
-        for n in [1, small_test_sizes...] 
+        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
@@ -98,7 +97,7 @@
             end
             c = contractor(s, subset)
             count = 0
-            for stabilizer in s 
+            for stabilizer in s
                 contractable = true
                 for i in subset
                     if stabilizer[i] != (false, false) contractable = false end
@@ -111,9 +110,9 @@
                     p = zero(PauliOperator, nqubits(s))
                     index = 0
                     for i in 1:nqubits(s)
-                        if !(i in subset) 
-                            index+=1 
-                            p[i] = contracted[index] 
+                        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
@@ -123,4 +122,4 @@
             end
         end
     end
-end
\ No newline at end of file
+end

From 617843c81a4362d1ef02997311dd771a196657c2 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Thu, 24 Oct 2024 20:58:55 -0400
Subject: [PATCH 03/15] more tests for sZCrY

---
 test/test_symcontrolled.jl | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/test/test_symcontrolled.jl b/test/test_symcontrolled.jl
index 256f62be3..4fbd86a87 100644
--- a/test/test_symcontrolled.jl
+++ b/test/test_symcontrolled.jl
@@ -125,4 +125,21 @@
             end
         end
     end
+
+    @testset "Ket-based definition for Y and rY" begin
+        for control in (:Z,)
+            target = :Y
+            for r in (true, false)
+                s = Stabilizer(QuantumClifford._T_str(string(control)))
+                k1 = Ket(s)
+                s.tab.phases[1] = 0x2
+                k2 = Ket(s)
+                i = Operator(tId1)
+                o = Operator(CliffordOperator(eval(Symbol(:s,target,))(1),1))
+                gate = projector(k1)⊗i + (r ? -1 : -im) * projector(k2)⊗o # XXX there is a -1 here because global phase is arbitrary and thus our Operator(tY) is NOT im*Operator(X)*Operator(Z)
+                implemented_gate = Operator(CliffordOperator(eval(Symbol(:s,control,:C,(r ? :rY : :Y)))(1,2),2))
+                @test gate≈implemented_gate
+            end
+        end
+    end
 end

From 0849b1c9f2d320ce3306961e15642da7dca83a6c Mon Sep 17 00:00:00 2001
From: Feroz Ahmad <feroz.ahmad.email@gmail.com>
Date: Fri, 25 Oct 2024 07:45:38 +0500
Subject: [PATCH 04/15] adding all missing types of two qubit SWAP gates and
 sSQRTZZ/sInvSQRTZZ gate (#336)

---------

Co-authored-by: Stefan Krastanov <github.acc@krastanov.org>
Co-authored-by: Stefan Krastanov <stefan@krastanov.org>
---
 CHANGELOG.md              |  5 ++++
 src/QuantumClifford.jl    |  3 ++-
 src/symbolic_cliffords.jl | 49 ++++++++++++++++++++++++++++-----------
 test/test_symcliff.jl     | 13 +++++++++++
 4 files changed, 55 insertions(+), 15 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e260d5c02..06856e122 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,11 @@
 
 # News
 
+## v0.9.13 - dev
+
+- Implementing additional named two-qubit gates: `sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP,
+    sSQRTZZ, sInvSQRTZZ`
+
 ## v0.9.12 - 2024-10-18
 
 - Minor compat fixes for julia 1.11 in the handling of `hgp`
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 32ca34a06..6837483cb 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -51,7 +51,8 @@ export
     sHadamardXY, sHadamardYZ, sSQRTX, sInvSQRTX, sSQRTY, sInvSQRTY, sCXYZ, sCZYX,
     sCNOT, sCPHASE, sSWAP,
     sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZCX, sZCY, sZCZ,
-    sZCrY, sInvZCrY,
+    sZCrY, sInvZCrY, sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP,
+    sSQRTZZ, sInvSQRTZZ,
     # Misc Ops
     SparseGate,
     sMX, sMY, sMZ, PauliMeasurement, Reset, sMRX, sMRY, sMRZ,
diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl
index 601696855..98f0e5d61 100644
--- a/src/symbolic_cliffords.jl
+++ b/src/symbolic_cliffords.jl
@@ -308,6 +308,16 @@ macro qubitop2(name, kernel)
 end
 #                 x1   z1      x2      z2
 @qubitop2 SWAP   (x2 , z2    , x1    , z1    , false)
+
+@qubitop2 SWAPCX    (x2    , z2⊻z1 , x2⊻x1 , z1    , ~iszero((x1 & z1 & x2 & z2) | (~x1 & z1 & x2 & ~z2)))
+@qubitop2 InvSWAPCX (x2⊻x1 , z2    , x1    , z2⊻z1 , ~iszero((x1 & z1 & x2 & z2) | ( x1 &~z1 &~x2 &  z2)))
+
+@qubitop2 ISWAP    (x2       , x1⊻z2⊻x2 , x1       , x1⊻x2⊻z1 , ~iszero((x1 & z1 & ~x2) | (~x1 & x2 & z2)))
+@qubitop2 InvISWAP (x2       , x1⊻z2⊻x2 , x1       , x1⊻x2⊻z1 , ~iszero((x1 &~z1 & ~x2) | (~x1 & x2 &~z2)))
+
+@qubitop2 CZSWAP (x2    , z2⊻x1 , x1    , x2⊻z1 , ~iszero((x1 & ~z1 & x2 & z2) | (x1 & z1 & x2 & ~z2)))
+@qubitop2 CXSWAP (x2⊻x1 , z2    , x1    , z2⊻z1 , ~iszero((x1 & ~z1 &~x2 & z2) | (x1 & z1 & x2 &  z2)))
+
 @qubitop2 CNOT   (x1 , z1⊻z2 , x2⊻x1 , z2    , ~iszero( (x1 & z1 & x2 & z2)  | (x1 & z2 &~(z1|x2)) ))
 @qubitop2 CPHASE (x1 , z1⊻x2 , x2    , z2⊻x1 , ~iszero( (x1 & z1 & x2 &~z2)  | (x1 &~z1 & x2 & z2) ))
 
@@ -326,6 +336,9 @@ end
 @qubitop2 ZCrY    (x1, x1⊻z1⊻x2⊻z2, x1⊻x2, x1⊻z2, ~iszero((x1 &~z1 & x2) | (x1 & ~z1 & ~z2) | (x1 & x2 & ~z2)))
 @qubitop2 InvZCrY (x1, x1⊻z1⊻x2⊻z2, x1⊻x2, x1⊻z2, ~iszero((x1 & z1 &~x2) | (x1 &  z1 &  z2) | (x1 &~x2 &  z2)))
 
+@qubitop2 SQRTZZ    (x1       , x1⊻x2⊻z1 , x2       , x1⊻z2⊻x2 , ~iszero((x1 & z1 & ~x2) | (~x1 & x2 & z2)))
+@qubitop2 InvSQRTZZ (x1       , x1⊻x2⊻z1 , x2       , x1⊻z2⊻x2 , ~iszero((x1 &~z1 & ~x2) | (~x1 & x2 &~z2)))
+
 #=
 To get the boolean formulas for the phase, it is easiest to first write down the truth table for the phase:
 for i in 0:15
@@ -370,20 +383,28 @@ function Base.show(io::IO, op::AbstractTwoQubitOperator)
     end
 end
 
-LinearAlgebra.inv(op::sSWAP) = sSWAP(op.q1, op.q2)
-LinearAlgebra.inv(op::sCNOT) = sCNOT(op.q1, op.q2)
-LinearAlgebra.inv(op::sCPHASE) = sCPHASE(op.q1, op.q2)
-LinearAlgebra.inv(op::sZCX) = sZCX(op.q1, op.q2)
-LinearAlgebra.inv(op::sZCY) = sZCY(op.q1, op.q2)
-LinearAlgebra.inv(op::sZCZ) = sZCZ(op.q1, op.q2)
-LinearAlgebra.inv(op::sXCX) = sXCX(op.q1, op.q2)
-LinearAlgebra.inv(op::sXCY) = sXCY(op.q1, op.q2)
-LinearAlgebra.inv(op::sXCZ) = sXCZ(op.q1, op.q2)
-LinearAlgebra.inv(op::sYCX) = sYCX(op.q1, op.q2)
-LinearAlgebra.inv(op::sYCY) = sYCY(op.q1, op.q2)
-LinearAlgebra.inv(op::sYCZ) = sYCZ(op.q1, op.q2)
-LinearAlgebra.inv(op::sZCrY) = sInvZCrY(op.q1, op.q2)
-LinearAlgebra.inv(op::sInvZCrY) = sZCrY(op.q1, op.q2)
+LinearAlgebra.inv(op::sSWAP)      = sSWAP(op.q1, op.q2)
+LinearAlgebra.inv(op::sCNOT)      = sCNOT(op.q1, op.q2)
+LinearAlgebra.inv(op::sCPHASE)    = sCPHASE(op.q1, op.q2)
+LinearAlgebra.inv(op::sZCX)       = sZCX(op.q1, op.q2)
+LinearAlgebra.inv(op::sZCY)       = sZCY(op.q1, op.q2)
+LinearAlgebra.inv(op::sZCZ)       = sZCZ(op.q1, op.q2)
+LinearAlgebra.inv(op::sXCX)       = sXCX(op.q1, op.q2)
+LinearAlgebra.inv(op::sXCY)       = sXCY(op.q1, op.q2)
+LinearAlgebra.inv(op::sXCZ)       = sXCZ(op.q1, op.q2)
+LinearAlgebra.inv(op::sYCX)       = sYCX(op.q1, op.q2)
+LinearAlgebra.inv(op::sYCY)       = sYCY(op.q1, op.q2)
+LinearAlgebra.inv(op::sYCZ)       = sYCZ(op.q1, op.q2)
+LinearAlgebra.inv(op::sZCrY)      = sInvZCrY(op.q1, op.q2)
+LinearAlgebra.inv(op::sInvZCrY)   = sZCrY(op.q1, op.q2)
+LinearAlgebra.inv(op::sSWAPCX)    = sInvSWAPCX(op.q1, op.q2)
+LinearAlgebra.inv(op::sInvSWAPCX) = sSWAPCX(op.q1, op.q2)
+LinearAlgebra.inv(op::sCZSWAP)    = sCZSWAP(op.q1, op.q2)
+LinearAlgebra.inv(op::sCXSWAP)    = sSWAPCX(op.q1, op.q2)
+LinearAlgebra.inv(op::sISWAP)     = sInvISWAP(op.q1, op.q2)
+LinearAlgebra.inv(op::sInvISWAP)  = sISWAP(op.q1, op.q2)
+LinearAlgebra.inv(op::sSQRTZZ)    = sInvSQRTZZ(op.q1, op.q2)
+LinearAlgebra.inv(op::sInvSQRTZZ) = sSQRTZZ(op.q1, op.q2)
 
 ##############################
 # Functions that perform direct application of common operators without needing an operator instance
diff --git a/test/test_symcliff.jl b/test/test_symcliff.jl
index 1d9f2dbe1..40ed505ac 100644
--- a/test/test_symcliff.jl
+++ b/test/test_symcliff.jl
@@ -99,6 +99,19 @@
             @test CliffordOperator(inv(sXCZ(n₁, n₂)), n₁) == inv(CliffordOperator(sCNOT(n₂, n₁), n₁))
             @test CliffordOperator(inv(sZCrY(n₁, n₂)), n₁) == inv(CliffordOperator(sZCrY(n₁, n₂), n₁))
             @test CliffordOperator(inv(sInvZCrY(n₁, n₂)), n₁) == inv(CliffordOperator(sInvZCrY(n₁, n₂), n₁))
+            @test CliffordOperator(inv(sCXSWAP(n₁, n₂)), n₁) == inv(CliffordOperator(sInvSWAPCX(n₁, n₂), n₁))
         end
     end
+
+    @testset "Consistency check with STIM conventions" begin
+        # see https://github.com/quantumlib/Stim/blob/main/doc/gates.md
+        @test CliffordOperator(sSWAPCX)    == C"IX XX ZZ ZI"
+        @test CliffordOperator(sInvSWAPCX) == C"XX XI IZ ZZ"
+        @test CliffordOperator(sCZSWAP)    == C"ZX XZ IZ ZI"
+        @test CliffordOperator(sCXSWAP)    == C"XX XI IZ ZZ"
+        @test CliffordOperator(sISWAP)     == C"ZY YZ IZ ZI"
+        @test CliffordOperator(sInvISWAP)  == C"-ZY -YZ IZ ZI"
+        @test CliffordOperator(sSQRTZZ)    == C"YZ ZY ZI IZ"
+        @test CliffordOperator(sInvSQRTZZ) == C"-YZ -ZY ZI IZ"
+    end
 end

From f3eb7cd7fd352fcc6f4994b7e42b7feca8112b7b Mon Sep 17 00:00:00 2001
From: IsaacP1234 <111547953+IsaacP1234@users.noreply.github.com>
Date: Thu, 24 Oct 2024 22:46:50 -0400
Subject: [PATCH 05/15] Additional group theory functions for error correction
 (#351)

---------

Co-authored-by: Kenneth Goodenough <kdgoodenough@gmail.com>
Co-authored-by: Stefan Krastanov <github.acc@krastanov.org>
Co-authored-by: Stefan Krastanov <stefan@krastanov.org>
---
 CHANGELOG.md                |   6 +
 docs/src/references.bib     |  25 +++
 src/QuantumClifford.jl      |   1 +
 src/grouptableaux.jl        | 296 ++++++++++++++++++++++++++++++++++--
 src/mul_leftright.jl        |   6 +
 src/project_trace_reset.jl  |   1 +
 test/test_group_tableaux.jl |  77 +++++++---
 7 files changed, 380 insertions(+), 32 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 06856e122..2fc286752 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,12 @@
 
 ## v0.9.13 - dev
 
+- New error correction group theory tools:
+    - `canonicalize_noncomm` function to find a generating set with minimal anticommutivity
+    - `SubsystemCodeTableau` data structure to represent the output of `canonicalize_noncomm`
+    - `commutify` function to find a commutative version of a non-commutative set of Paulis with minimal changes
+    - `matroid_parent` to, for set of Paulis that doesn't represent a state, find a version
+    that does.
 - Implementing additional named two-qubit gates: `sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP,
     sSQRTZZ, sInvSQRTZZ`
 
diff --git a/docs/src/references.bib b/docs/src/references.bib
index 29500a034..281ace9c1 100644
--- a/docs/src/references.bib
+++ b/docs/src/references.bib
@@ -403,6 +403,31 @@ @inproceedings{brown2013short
   doi = {10.1109/ISIT.2013.6620245}
 }
 
+@article{RevModPhys.87.307,
+  title = {Quantum error correction for quantum memories},
+  author = {Terhal, Barbara M.},
+  journal = {Rev. Mod. Phys.},
+  volume = {87},
+  issue = {2},
+  pages = {307--346},
+  numpages = {40},
+  year = {2015},
+  month = {Apr},
+  publisher = {American Physical Society},
+  doi = {10.1103/RevModPhys.87.307},
+  url = {https://link.aps.org/doi/10.1103/RevModPhys.87.307}
+}
+
+@misc{goodenough2024bipartiteentanglementnoisystabilizer,
+  title={Bipartite entanglement of noisy stabilizer states through the lens of stabilizer codes}, 
+  author={Kenneth Goodenough and Aqil Sajjad and Eneet Kaur and Saikat Guha and Don Towsley},
+  year={2024},
+  eprint={2406.02427},
+  archivePrefix={arXiv},
+  primaryClass={quant-ph},
+  url={https://arxiv.org/abs/2406.02427}, 
+}
+
 @article{panteleev2021degenerate,
   title = {Degenerate {{Quantum LDPC Codes With Good Finite Length Performance}}},
   author = {Panteleev, Pavel and Kalachev, Gleb},
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 6837483cb..9d892c958 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -78,6 +78,7 @@ export
     graphstate, graphstate!, graph_gatesequence, graph_gate,
     # Group theory tools
     groupify, minimal_generating_set, pauligroup, normalizer, centralizer, contractor, delete_columns,
+    canonicalize_noncomm, commutify, matroid_parent, SubsystemCodeTableau,
     # Clipped Gauge
     canonicalize_clip!, bigram, entanglement_entropy,
     # mctrajectories
diff --git a/src/grouptableaux.jl b/src/grouptableaux.jl
index 208f1ddfb..6f76802cf 100644
--- a/src/grouptableaux.jl
+++ b/src/grouptableaux.jl
@@ -1,3 +1,110 @@
+"""
+A tableau representation of the non-commutative canonical form of a set of Paulis,
+which is used in [`commutify`](@ref).
+
+They are organized in the same form as [`MixedDestabilizer`](@ref)
+with a stabilizer, destabilizer, logical X, and logical Z components.
+"""
+mutable struct SubsystemCodeTableau{T <: Tableau} <: AbstractStabilizer
+    tab::T
+    index::Int
+    r::Int
+    m::Int
+    k::Int
+end
+
+function SubsystemCodeTableau(t::Tableau)
+    index = 1
+    for i in range(1, stop=length(t), step=2)
+        if i + 1 > length(t)
+            break
+        end
+        if comm(t[i], t[i+1]) == 0x01
+            index = i+2 # index to split t into non-commuting pairs and commuting operators
+        end
+    end
+    s = Stabilizer(t[index:length(t)])
+    ind = 1
+    if length(s)>nqubits(s)#if stabilizer is overdetermined, Destabilizer constructor throws error
+        m = length(s)
+        tab = zero(Tableau, length(t)+length(s), nqubits(t))
+        for i in s
+            tab[ind] = zero(PauliOperator, nqubits(s))
+            ind+=1
+        end
+    else
+        d = Destabilizer(s)
+        m = length(d)
+        tab = zero(Tableau, length(t)+length(d), nqubits(t))
+        for p in destabilizerview(d)
+            tab[ind] = p
+            ind+=1
+        end
+    end
+    for i in range(1, stop=index-1, step=2)
+        tab[ind] = t[i]
+        ind+=1
+    end
+    for p in s
+        tab[ind] = p
+        ind+=1
+    end
+    for i in range(2, stop=index, step=2)
+        tab[ind] = t[i]
+        ind+=1
+    end
+    return SubsystemCodeTableau(tab, index, length(s), m, Int((index-1)/2))
+end
+
+
+Base.copy(t::SubsystemCodeTableau) = SubsystemCodeTableau(copy(t.tab), t.index, t.r, t.m, t.k)
+
+function Base.show(io::IO, t::SubsystemCodeTableau)
+    r = t.r+t.m+2*t.k 
+    q = nqubits(t)
+    if get(io, :compact, false) | haskey(io, :typeinfo)
+        print(io, "MixedDestablizer $r×$q")
+    elseif get(io, :limit, false)
+        h,w = displaysize(io)
+        println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(min(w-9,size(t.tab,2)-4),0))
+        _show(io, destabilizerview(t).tab, w, h÷4)
+        if r != q
+            println(io)
+            println(io, "𝒳" * "━"^max(min(w-5,size(t.tab,2)),0))
+            _show(io, logicalxview(t).tab, w, h÷4)
+        end
+        println(io)
+        println(io, "𝒮𝓉𝒶𝒷" * "━"^max(min(w-7,size(t.tab,2)-2),0))
+        _show(io, stabilizerview(t).tab, w, h÷4)
+        if r != q
+            println(io)
+            println(io, "𝒵" * "━"^max(min(w-5,size(t.tab,2)),0))
+            _show(io, logicalzview(t).tab, w, h÷4)
+        end
+    else
+        println(io, "𝒟ℯ𝓈𝓉𝒶𝒷" * "━"^max(size(t.tab,2)-4,0))
+        _show(io, destabilizerview(t).tab, missing, missing)
+        if r != q
+            println(io)
+            println(io, "𝒳ₗ" * "━"^max(size(t.tab,2),0))
+            _show(io, logicalxview(t).tab, missing, missing)
+        end
+        println(io)
+        println(io, "𝒮𝓉𝒶𝒷" * "━"^max(size(t.tab,2)-2,0))
+        _show(io, stabilizerview(t).tab, missing, missing)
+        if r != q
+            println(io)
+            println(io, "𝒵ₗ" * "━"^max(size(t.tab,2)))
+            _show(io, logicalzview(t).tab, missing, missing)
+        end
+    end
+end
+
+@inline stabilizerview(s::SubsystemCodeTableau) = Stabilizer(@view tab(s)[s.m+s.k+1:s.m+s.k+s.r])
+@inline destabilizerview(s::SubsystemCodeTableau) = Stabilizer(@view tab(s)[1:s.r])
+@inline logicalxview(s::SubsystemCodeTableau) = Stabilizer(tab(s)[s.m+1:s.m+s.k])
+@inline logicalzview(s::SubsystemCodeTableau) = Stabilizer(tab(s)[length(s.tab)-s.k+1:length(s.tab)])
+
 """
 Return the full stabilizer group represented by the input generating set (a [`Stabilizer`](@ref)).
 
@@ -11,16 +118,20 @@ julia> groupify(S"XZ ZX")
 + YY
 ```
 """
-function groupify(s::Stabilizer)
+function groupify(s::Stabilizer; phases=false)
     # 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))
+    if phases == true
+        throw(ArgumentError("in groupify phases=true functionality not yet implemented"))
+    end
+    gen_set = minimal_generating_set(s)
+    n = length(gen_set)::Int
+    group = zero(Tableau, 2^n, nqubits(gen_set))
     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]
+                mul_left!(group, i+1, tab(gen_set), digit_order, phases=Val(false))
             end
         end
     end
@@ -48,7 +159,7 @@ function minimal_generating_set(s::Stabilizer)
     if r == 0
         gs = zero(Stabilizer, 1, nqubits(s))
         if 0x02 in phases(s)
-            gs[1] = -1 * gs[1]
+            phases(gs)[1] = 0x2
         end
         return gs
     else
@@ -56,6 +167,153 @@ function minimal_generating_set(s::Stabilizer)
     end
 end
 
+"""
+For a not-necessarily commutative set of Paulis, return a generating set of the form
+⟨A₁, A₂, ... Aₖ, Aₖ₊₁, ... Aₘ, B₁, B₂, ... Bₖ⟩ where pairs Aₖ, Bₖ anticommute and all other pairings commute. Based on [RevModPhys.87.307](@cite)
+
+Returns the generating set as a data structure of type [`SubsystemCodeTableau`](@ref). The
+`logicalxview` function returns the ⟨A₁, A₂,... Aₖ⟩, and the `logicalzview`
+function returns ⟨B₁, B₂, ... Bₖ⟩. `stabilizerview` returns ⟨Aₖ₊₁, ... Aₘ⟩
+as a Stabilizer, and `destabilizerview` returns the Destabilizer of that Stabilizer.
+
+Phases are zeroed-out in this canonicalization.
+
+```jldoctest
+julia> canonicalize_noncomm(T"XX XZ XY")
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z_
+𝒳━━
++ XX
+𝒮𝓉𝒶𝒷
++ X_
+𝒵━━
++ XZ
+```
+"""
+function canonicalize_noncomm(t::Tableau)
+    for i in eachindex(t) phases(t)[i] = 0x00 end
+    loc = zero(Tableau, 2*length(t), nqubits(t))
+    x_index = 0
+    z_index = length(loc)
+    for i in eachindex(t)
+        for j in eachindex(t)
+            if comm(t[i], t[j]) == 0x01
+                for k in eachindex(t)
+                    if k !=i && k != j
+                        if comm(t[k], t[i]) == 0x01
+                            mul_left!(t, k, j, phases=Val(false))
+                        end
+                        if comm(t[k], t[j]) == 0x01
+                            mul_left!(t, k, i, phases=Val(false))
+                        end
+                    end
+                end
+                if !(t[i] in loc)
+                    x_index+=1
+                    loc[x_index]= t[i]
+                end
+                if !(t[j] in loc)
+                    loc[z_index]= t[j]
+                    z_index-=1
+                end
+            end
+        end
+    end
+    ind = 2*x_index+1 #format tableau for SubsystemCodeTableau
+    k = x_index
+    m = length(t)-2*k
+    for i in range(k, stop =1, step=-1)
+        loc[i+m]= loc[i]
+    end
+    i = m+1
+    j = m+k
+    while i < j #ensure anticommutative pairs are lined up correctly
+        loc[i], loc[j], = loc[j], loc[i]
+        i+=1
+        j-=1
+    end
+    s_index = m+k+1
+    for i in eachindex(t)
+        if !(t[i] in loc)
+            loc[s_index]= t[i]
+            s_index+=1
+        end
+    end
+    r = s_index-m-k-1
+    if r <= nqubits(t)
+        d = destabilizerview(Destabilizer(Stabilizer(loc[m+k+1:s_index-1])))
+        for i in eachindex(d)
+            loc[i] =d[i]
+        end
+    else for i in 1:m zero!(loc, i) end end
+    return SubsystemCodeTableau(loc, ind, r, m, k)
+end
+
+"""
+For a not-necessarily commutative set of Paulis S,
+computed S', the [non-commutative canonical form](@ref canonicalize_noncomm) of of S.
+For each pair Aₖ, Bₖ of anticommutative Paulis in S', add a qubit to each Pauli in the set:
+X to Aₖ, Z to Bₖ, and I to each other operator to produce S'', a fully commutative set. Return
+S'' as well as a list of the indices of the added qubits.
+
+The returned object is a Stabilizer that is also useful for the [`matroid_parent`](@ref) function.
+
+```jldoctest
+julia> commutify(T"XX XZ XY")[1]
++ XXX
++ X__
++ XZZ
+
+julia> commutify(T"XX XZ XY")[2]
+3:3
+```
+"""
+function commutify(t)
+    loc = canonicalize_noncomm(t)
+    commutative = zero(Stabilizer, 2*loc.k+loc.r, nqubits(loc)+loc.k)
+    puttableau!(commutative, logicalxview(loc), 0,0)
+    puttableau!(commutative, stabilizerview(loc), loc.k,0)
+    puttableau!(commutative, logicalzview(loc), loc.k+loc.r,0)
+    x= T"X"
+    z=T"Z"
+    for i in 0:(loc.k-1)
+        puttableau!(commutative, x, i, nqubits(loc)+i)
+        puttableau!(commutative, z, i+loc.r+loc.k, nqubits(loc)+i)
+    end
+    to_delete = nqubits(loc)+1:nqubits(loc)+loc.k
+    return commutative, to_delete
+end
+
+"""
+For a given set S of Paulis that does not necessarily represent a state,
+return a set of Paulis S' that represents a state.
+S' is a superset of [commutified](@ref commutify) S.
+Additionally returns two arrays representing deletions needed to produce S.
+Based on [goodenough2024bipartiteentanglementnoisystabilizer](@cite)
+
+By deleting the qubits in the first output array from S', taking the [`normalizer`](@ref) of S', then
+deleting the qubits in the second returned array from the [`normalizer`](@ref) of S', S is reproduced.
+
+```jldoctest
+julia> matroid_parent(T"XX")[1]
++ X_X
++ XX_
++ ZZZ
+
+julia> matroid_parent(T"XX")[2]
+3:3
+
+julia> matroid_parent(T"XX")[3]
+3:2
+```
+"""
+function matroid_parent(t::Tableau)
+    com, d1= commutify(t)
+    norm = normalizer(com.tab)
+    state, d2 = commutify(norm)
+    return state, d2, d1
+end
+
 """
 Return the full Pauli group of a given length. Phases are ignored by default,
 but can be included by setting `phases=true`.
@@ -127,15 +385,21 @@ julia> normalizer(T"X")
 + X
 ```
 """
-function normalizer(t::Tableau)
+function normalizer(t::Tableau; phases=false)
     # 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
+    ptype =typeof(P"I")
+    norm = ptype[]
+
+    p = zero(PauliOperator, n)
+    paulis = ((false, false), (true, false), (false, true), (true, true))
+    for i in Iterators.product(Iterators.repeated(paulis, n)...)
+        zero!(p)
+        for (j, k) in enumerate(i)
+            p[j] = k
+        end
         commutes = true
         for q in t
             if comm(p, q) == 0x01
@@ -143,11 +407,15 @@ function normalizer(t::Tableau)
             end
         end
         if commutes
-            push!(normalizer, p)
+            push!(norm, copy(p))
+        end
+        if phases
+            for phase in [-1, 1im, -1im]
+                push!(norm, phase *p)
+            end
         end
     end
-
-    return Tableau(normalizer)
+    return Tableau(norm)
 end
 
 """
@@ -181,7 +449,7 @@ 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.
+the given subset, without the entries corresponding to subset. Based on [goodenough2024bipartiteentanglementnoisystabilizer](@cite)
 
 ```jldoctest
 julia> contractor(S"_X X_", [1])
diff --git a/src/mul_leftright.jl b/src/mul_leftright.jl
index 9803ed2b5..08f905664 100644
--- a/src/mul_leftright.jl
+++ b/src/mul_leftright.jl
@@ -157,6 +157,12 @@ end
 # On Tableaux
 ##############################
 
+@inline function mul_left!(s::Tableau, m, t::Tableau, i; phases::Val{B}=Val(true)) where B
+    extra_phase = mul_left!((@view s.xzs[:,m]), (@view t.xzs[:,i]); phases=phases)
+    B && (s.phases[m] = (extra_phase+s.phases[m]+s.phases[i])&0x3)
+    s
+end
+
 @inline function mul_left!(s::Tableau, m, i; phases::Val{B}=Val(true)) where B
     extra_phase = mul_left!((@view s.xzs[:,m]), (@view s.xzs[:,i]); phases=phases)
     B && (s.phases[m] = (extra_phase+s.phases[m]+s.phases[i])&0x3)
diff --git a/src/project_trace_reset.jl b/src/project_trace_reset.jl
index 4b2e6e2e5..25aed09ce 100644
--- a/src/project_trace_reset.jl
+++ b/src/project_trace_reset.jl
@@ -885,5 +885,6 @@ julia> delete_columns(S"XYZ YZX ZXY", [1,3])
 See also: [`traceout!`](@ref)
 """
 function delete_columns(𝒮::Stabilizer, subset)
+    if length(𝒮) == 0 return 𝒮 end
     return 𝒮[:, setdiff(1:nqubits(𝒮), subset)]
 end
diff --git a/test/test_group_tableaux.jl b/test/test_group_tableaux.jl
index d208afdff..fb3240441 100644
--- a/test/test_group_tableaux.jl
+++ b/test/test_group_tableaux.jl
@@ -1,4 +1,4 @@
-@testitem "group theory routines" begin
+@testitem "group theory tools" begin
     using Test
     using Random
     using QuantumClifford
@@ -8,41 +8,82 @@
     # 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
+    @testset "group theory tools" 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
+                @test s == s_test
             end
         end
         #Test minimal_generating_set
-        for n in [1, small_test_sizes...]
+        for n in [1, test_sizes...]
             s = random_stabilizer(n)
             group = groupify(s)
             gen_set = minimal_generating_set(Stabilizer(group))
+            @test length(group) == 2^(length(gen_set))
             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)
+            s1, _, r = canonicalize!(Stabilizer(group), ranks = true)
+            s2, _, r = canonicalize!(Stabilizer(new_group), ranks=true)
+            @test group[1:r, :] == new_group[1:r, :]
+        end
+        #Test canonicalize_noncomm
+        for n in [1, small_test_sizes...]
+            t = zero(QuantumClifford.Tableau, rand(1:(2*n)), n)
+            for i in eachindex(t) t[i] = random_pauli(n) end
+            loc = canonicalize_noncomm(t)
+            for i in 1:loc.k
+                for j in 1:loc.k
+                    if i == j
+                        @test comm(logicalxview(loc)[i], logicalzview(loc)[j]) == 0x01
+                    else @test comm(logicalxview(loc)[i], logicalzview(loc)[j]) == 0x00 end
+                end
             end
-            gen_set = minimal_generating_set(s)
-            new_group = groupify(s)
-            for operator in s
-                @test operator in new_group
+            for i in stabilizerview(loc)
+                for j in stabilizerview(loc) @test comm(i, j) == 0x00 end
+                for j in logicalzview(loc) @test comm(i, j) == 0x00 end
+                for j in logicalxview(loc) @test comm(i, j) == 0x00 end
             end
         end
+        #Test commutify
+        for n in [1, small_test_sizes...]
+            t = zero(QuantumClifford.Tableau, rand(1:(2*n)), n)
+            for i in eachindex(t) t[i] = random_pauli(n) end
+            c, d = commutify(t)
+            for i in c
+                for j in c
+                    @test comm(i, j) == 0x00
+                end
+            end
+            for i in d
+                for p in c
+                    @test p[i] != (true, true)
+                end
+            end
+            loc1= delete_columns(c, d)
+            loc2 = canonicalize_noncomm(t).tab
+            for i in eachindex(delete_columns(c, d))
+            end
+        end
+        #Test matroid_parent
+        for n in [1,2,3,4,5]
+            t = zero(QuantumClifford.Tableau, 2*n, n)
+            for i in eachindex(t) t[i] = random_pauli(n) end
+            e, d2, d1 = matroid_parent(t)
+            s = Stabilizer(groupify(e))
+            for i in e for j in e @test comm(i, j)==0x00 end end
+            @test 2^(nqubits(s)) == length(s) #assumes commutativise works
+            #find original tableau from matroid_parentded state, ignoring phases
+            inverted = delete_columns(Stabilizer(normalizer(delete_columns(Stabilizer(e), d2).tab)), d1)
+            original = Stabilizer(groupify(Stabilizer(t)))
+            canonicalize!(inverted)
+            canonicalize!(original)
+            @test inverted[1:length(inverted)].tab.xzs == original[1:length(inverted)].tab.xzs
+        end
         #Test pauligroup
         for n in [1, small_test_sizes...]
             @test length(QuantumClifford.pauligroup(n, phases=false)) == 4^n

From 47b8d4631c3f64df084c8a5c76edb37ce0eadc9b Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Wed, 30 Oct 2024 18:17:19 -0400
Subject: [PATCH 06/15] use explicit imports in ECC (#404)

---
 CHANGELOG.md   |  7 +++----
 Project.toml   |  2 +-
 src/ecc/ECC.jl | 15 +++++++++++----
 3 files changed, 15 insertions(+), 9 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2fc286752..9c5077104 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,16 +5,15 @@
 
 # News
 
-## v0.9.13 - dev
+## v0.9.13 - 2024-10-30
 
-- New error correction group theory tools:
+- New error-correction group theory tools:
     - `canonicalize_noncomm` function to find a generating set with minimal anticommutivity
     - `SubsystemCodeTableau` data structure to represent the output of `canonicalize_noncomm`
     - `commutify` function to find a commutative version of a non-commutative set of Paulis with minimal changes
     - `matroid_parent` to, for set of Paulis that doesn't represent a state, find a version
     that does.
-- Implementing additional named two-qubit gates: `sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP,
-    sSQRTZZ, sInvSQRTZZ`
+- Implementing additional named two-qubit gates: `sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP, sSQRTZZ, sInvSQRTZZ`
 
 ## v0.9.12 - 2024-10-18
 
diff --git a/Project.toml b/Project.toml
index 4d300bc66..83e608beb 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.12"
+version = "0.9.13"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl
index cdda7742e..ff4444694 100644
--- a/src/ecc/ECC.jl
+++ b/src/ecc/ECC.jl
@@ -1,9 +1,16 @@
 module ECC
 
-using LinearAlgebra
-using LinearAlgebra: I
-using QuantumClifford
-using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer
+using LinearAlgebra: LinearAlgebra, I, rank, tr
+using QuantumClifford: QuantumClifford, AbstractOperation, AbstractStabilizer,
+    AbstractTwoQubitOperator, Stabilizer, PauliOperator,
+    random_brickwork_clifford_circuit, random_all_to_all_clifford_circuit,
+    canonicalize!, canonicalize_gott!,
+    logicalxview, logicalzview, stabilizerview, destabilizerview, tab, phases,
+    sCNOT, sSWAP, sHadamard, sPhase, sInvPhase,
+    sZCX, sZCY, sZCZ, sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZ, sX, sY, sMRZ, sMRX,
+    single_x, single_y, single_z, random_pauli!, PauliError,
+    apply!, comm, comm!, stab_to_gf2, embed, @S_str, affectedqubits, affectedbits,
+    pftrajectories, pfmeasurements, mctrajectories
 import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
 using DocStringExtensions
 using Combinatorics: combinations

From f69e49ccfb91d7fd85bda16a846c34d22c500399 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Quinten=20Prei=C3=9F?=
 <90014830+J-C-Q@users.noreply.github.com>
Date: Thu, 31 Oct 2024 16:30:30 +0100
Subject: [PATCH 07/15] Update noisycircuits_mc.md (#410)

---
 docs/src/noisycircuits_mc.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/noisycircuits_mc.md b/docs/src/noisycircuits_mc.md
index 37e900df4..d24585573 100644
--- a/docs/src/noisycircuits_mc.md
+++ b/docs/src/noisycircuits_mc.md
@@ -15,7 +15,7 @@ Import with `using QuantumClifford.Experimental.NoisyCircuits`.
 
 This module enables the simulation of noisy Clifford circuits through a Monte Carlo method where the same circuit is evaluated multiple times with random errors interspersed through it as prescribed by a given error model.
 
-Below is an example of a purification circuit. We first prepare the circuit we desire to use, including a noise model. `Quantikz.jl` was is used to visualize the circuit.
+Below is an example of a purification circuit. We first prepare the circuit we desire to use, including a noise model. `Quantikz.jl` is used to visualize the circuit.
 
 ```@example 1
 using QuantumClifford # hide
@@ -55,8 +55,8 @@ If you want to create a custom gate type (e.g. calling it `Operation`), you need
 The `Symbol` is the status of the operation. Predefined statuses are kept in the `registered_statuses` list, but you can add more.
 Be sure to expand this list if you want the trajectory simulators using your custom statuses to output all trajectories.
 
-There is also [`applynoise!`](@ref) which is convenient wait to create a noise model that can then be plugged into the [`NoisyGate`](@ref) struct,
+There is also [`applynoise!`](@ref) which is a convenient way to create a noise model that can then be plugged into the [`NoisyGate`](@ref) struct,
 letting you reuse the predefined perfect gates and measurements.
 However, you can also just make up your own noise operator simply by implementing [`applywstatus!`](@ref) for it.
 
-You can also consult the [list of implemented operators](@ref noisycircuits_ops).
\ No newline at end of file
+You can also consult the [list of implemented operators](@ref noisycircuits_ops).

From abfcd3cc7c4b7207d6964e3d607b46b946755138 Mon Sep 17 00:00:00 2001
From: "github-actions[bot]"
 <41898282+github-actions[bot]@users.noreply.github.com>
Date: Fri, 1 Nov 2024 23:02:25 -0400
Subject: [PATCH 08/15] CompatHelper: bump compat for ILog2 to 2, (keep
 existing compat) (#380)

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
Co-authored-by: Stefan Krastanov <github.acc@krastanov.org>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 83e608beb..81df346bd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -50,7 +50,7 @@ DocStringExtensions = "0.9"
 Graphs = "1.9"
 Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
 HostCPUFeatures = "0.1.6"
-ILog2 = "0.2.3"
+ILog2 = "0.2.3, 1, 2"
 InteractiveUtils = "1.9"
 LDPCDecoders = "0.3.1"
 LinearAlgebra = "1.9"

From 4cc664eb6cea74703a402ecc79845f1fcb6078a1 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Sat, 2 Nov 2024 10:23:08 -0400
Subject: [PATCH 09/15] fix to affectedbit for symbolic measurements (#414)

---
 CHANGELOG.md          |  4 ++++
 Project.toml          |  2 +-
 src/affectedqubits.jl |  3 +--
 test/test_symmeas.jl  | 35 +++++++++++++++++++++++++++++++++++
 4 files changed, 41 insertions(+), 3 deletions(-)
 create mode 100644 test/test_symmeas.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9c5077104..e324370b0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,10 @@
 
 # News
 
+## v0.9.14 - 2024-11-02
+
+- **(fix)** `affectedqubits()` on `sMX`, `sMY`, and `sMR*`
+
 ## v0.9.13 - 2024-10-30
 
 - New error-correction group theory tools:
diff --git a/Project.toml b/Project.toml
index 81df346bd..e49f22503 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.13"
+version = "0.9.14"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/src/affectedqubits.jl b/src/affectedqubits.jl
index da6314a65..3609710e3 100644
--- a/src/affectedqubits.jl
+++ b/src/affectedqubits.jl
@@ -15,6 +15,5 @@ affectedqubits(c::CliffordOperator) = 1:nqubits(c)
 affectedqubits(c::ClassicalXOR) = ()
 
 affectedbits(o) = ()
-affectedbits(m::sMRZ) = (m.bit,)
-affectedbits(m::sMZ) = (m.bit,)
+affectedbits(m::Union{sMRZ,sMZ,sMRX,sMX,sMRY,sMY}) = m.bit==0 ? () : (m.bit,)
 affectedbits(c::ClassicalXOR) = (c.bits..., c.store)
diff --git a/test/test_symmeas.jl b/test/test_symmeas.jl
new file mode 100644
index 000000000..ebe035a95
--- /dev/null
+++ b/test/test_symmeas.jl
@@ -0,0 +1,35 @@
+@testitem "Symbolic Measurements" begin
+    using Random
+    using QuantumClifford
+    using Test
+    test_sizes = [10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
+
+    @testset "Symbolic Measurements" begin
+        for n in test_sizes
+            for i in [rand(1:n), 1, n, n÷2+1, n÷2-1, n÷2]
+                for T in [sMRZ,sMZ,sMRX,sMX] # TODO sMY sMRY
+                    for _ in 1:10
+                        apply!(PauliFrame(n,n,n), T(i,1))
+                    end
+                end
+            end
+            for i in [rand(1:n), 1, n, n÷2+1, n÷2-1, n÷2]
+                for T in [sMRZ,sMZ,sMRX,sMX,sMY] # TODO sMRY
+                    for _ in 1:10
+                        s = random_stabilizer(n)
+                        apply!(Register(copy(s),1), T(i,1))
+                    end
+                end
+            end
+        end
+        @test_throws DimensionMismatch SingleQubitOperator(tCNOT,1)
+        @test_throws DimensionMismatch CliffordOperator(sHadamard(5),2)
+        @test_throws ArgumentError CliffordOperator(sHadamard(5),6,compact=true)
+        for T in [sMRZ,sMZ,sMRX,sMX,sMRY,sMY]
+            @test_throws ArgumentError T(-1)
+            @test_throws ArgumentError T(-1,0)
+            @test_throws ArgumentError T(-1,1)
+            @test T(1,0) == T(1)
+        end
+    end
+end

From d3f42d23f62d14b227adff705898d4a343d744b5 Mon Sep 17 00:00:00 2001
From: Feroz Ahmad <feroz.ahmad.email@gmail.com>
Date: Sat, 2 Nov 2024 19:35:34 +0500
Subject: [PATCH 10/15] better bit-wrangling abstraction using less boilerplate
 (#367)

Co-authored-by: Stefan Krastanov <stefan@krastanov.org>
---
 ext/QuantumCliffordGPUExt/apply_noise.jl  | 11 ++----
 ext/QuantumCliffordGPUExt/pauli_frames.jl | 12 ++----
 src/QuantumClifford.jl                    | 24 ++++++++++++
 src/pauli_frames.jl                       | 46 ++++++++---------------
 src/pauli_operator.jl                     | 14 ++++---
 src/project_trace_reset.jl                | 15 +++-----
 src/symbolic_cliffords.jl                 | 21 +++--------
 7 files changed, 68 insertions(+), 75 deletions(-)

diff --git a/ext/QuantumCliffordGPUExt/apply_noise.jl b/ext/QuantumCliffordGPUExt/apply_noise.jl
index 4aa99ed4b..201eb7c60 100644
--- a/ext/QuantumCliffordGPUExt/apply_noise.jl
+++ b/ext/QuantumCliffordGPUExt/apply_noise.jl
@@ -1,16 +1,11 @@
-using QuantumClifford: _div, _mod
+using QuantumClifford: get_bitmask_idxs
 
 #according to https://github.com/JuliaGPU/CUDA.jl/blob/ac1bc29a118e7be56d9edb084a4dea4224c1d707/test/core/device/random.jl#L33
 #CUDA.jl supports calling rand() inside kernel
 function applynoise!(frame::PauliFrameGPU{T},noise::UnbiasedUncorrelatedNoise,i::Int) where {T <: Unsigned}
     p = noise.p
-    lowbit = T(1)
-    ibig = _div(T,i-1)+1
-    ismall = _mod(T,i-1)
-    ismallm = lowbit<<(ismall)
-
-    stab = frame.frame
-    xzs = tab(stab).xzs
+    xzs = tab(frame.frame).xzs
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
     rows = size(stab, 1)
 
     @run_cuda applynoise_kernel(xzs, p, ibig, ismallm, rows) rows
diff --git a/ext/QuantumCliffordGPUExt/pauli_frames.jl b/ext/QuantumCliffordGPUExt/pauli_frames.jl
index 34a0565fd..4d6e03654 100644
--- a/ext/QuantumCliffordGPUExt/pauli_frames.jl
+++ b/ext/QuantumCliffordGPUExt/pauli_frames.jl
@@ -1,3 +1,5 @@
+using QuantumClifford: get_bitmask_idxs
+
 ##############################
 # sMZ
 ##############################
@@ -21,10 +23,7 @@ function apply!(frame::PauliFrameGPU{T}, op::QuantumClifford.sMZ) where {T <: Un
     op.bit == 0 && return frame
     i = op.qubit
     xzs = frame.frame.tab.xzs
-    lowbit = T(1)
-    ibig = QuantumClifford._div(T,i-1)+1
-    ismall = QuantumClifford._mod(T,i-1)
-    ismallm = lowbit<<(ismall)
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
     (@run_cuda apply_sMZ_kernel!(xzs, frame.measurements, op, ibig, ismallm, length(frame)) length(frame))
     return frame
 end
@@ -55,10 +54,7 @@ end
 function apply!(frame::PauliFrameGPU{T}, op::QuantumClifford.sMRZ) where {T <: Unsigned} # TODO sMRX, sMRY
     i = op.qubit
     xzs = frame.frame.tab.xzs
-    lowbit = T(1)
-    ibig = QuantumClifford._div(T,i-1)+1
-    ismall = QuantumClifford._mod(T,i-1)
-    ismallm = lowbit<<(ismall)
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
     (@run_cuda apply_sMRZ_kernel!(xzs, frame.measurements, op, ibig, ismallm, length(frame)) length(frame))
     return frame
 end
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 9d892c958..6b4945b6b 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -912,6 +912,30 @@ function unsafe_bitfindnext_(chunks::AbstractVector{T}, start::Int) where T<:Uns
     return nothing
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Computes bitmask indices for an unsigned integer at index `i`
+within the binary structure of a `Tableau` or `PauliOperator`.
+
+For `Tableau`, the method operates on the `.xzs` field, while
+for `PauliOperator`, it uses the `.xz` field. It calculates
+the following values based on the index `i`:
+
+- `lowbit`, the lowest bit.
+- `ibig`, the index of the word containing the bit.
+- `ismall`, the position of the bit within the word.
+- `ismallm`, a bitmask isolating the specified bit.
+"""
+@inline function get_bitmask_idxs(xzs::AbstractArray{<:Unsigned}, i::Int)
+    T = eltype(xzs)
+    lowbit = T(1)
+    ibig = _div(T, i-1) + 1
+    ismall = _mod(T, i-1)
+    ismallm = lowbit << ismall
+    return lowbit, ibig, ismall, ismallm
+end
+
 """Permute the qubits (i.e., columns) of the tableau in place."""
 function Base.permute!(s::Tableau, perm::AbstractVector)
     for r in 1:size(s,1)
diff --git a/src/pauli_frames.jl b/src/pauli_frames.jl
index 45e9763c8..af84bd63b 100644
--- a/src/pauli_frames.jl
+++ b/src/pauli_frames.jl
@@ -82,12 +82,8 @@ end
 function apply!(frame::PauliFrame, op::sMZ) # TODO sMY, and faster sMX
     op.bit == 0 && return frame
     i = op.qubit
-    xzs = frame.frame.tab.xzs
-    T = eltype(xzs)
-    lowbit = T(1)
-    ibig = _div(T,i-1)+1
-    ismall = _mod(T,i-1)
-    ismallm = lowbit<<(ismall)
+    xzs = tab(frame.frame).xzs
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
 
     @inbounds @simd for f in eachindex(frame)
         should_flip = !iszero(xzs[ibig,f] & ismallm)
@@ -99,12 +95,8 @@ end
 
 function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRY, and faster sMRX
     i = op.qubit
-    xzs = frame.frame.tab.xzs
-    T = eltype(xzs)
-    lowbit = T(1)
-    ibig = _div(T,i-1)+1
-    ismall = _mod(T,i-1)
-    ismallm = lowbit<<(ismall)
+    xzs = tab(frame.frame).xzs
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
 
     if op.bit != 0
         @inbounds @simd for f in eachindex(frame)
@@ -122,45 +114,39 @@ end
 
 function applynoise!(frame::PauliFrame,noise::UnbiasedUncorrelatedNoise,i::Int)
     p = noise.p
-    T = eltype(frame.frame.tab.xzs)
+    xzs = tab(frame.frame).xzs
 
-    lowbit = T(1)
-    ibig = _div(T,i-1)+1
-    ismall = _mod(T,i-1)
-    ismallm = lowbit<<(ismall)
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
     p = p/3
 
     @inbounds @simd for f in eachindex(frame)
         r = rand()
         if  r < p # X error
-            frame.frame.tab.xzs[ibig,f] ⊻= ismallm
+            xzs[ibig,f] ⊻= ismallm
         elseif r < 2p # Z error
-            frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm
+            xzs[end÷2+ibig,f] ⊻= ismallm
         elseif r < 3p # Y error
-            frame.frame.tab.xzs[ibig,f] ⊻= ismallm
-            frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm
+            xzs[ibig,f] ⊻= ismallm
+            xzs[end÷2+ibig,f] ⊻= ismallm
         end
     end
     return frame
 end
 
 function applynoise!(frame::PauliFrame,noise::PauliNoise,i::Int)
-    T = eltype(frame.frame.tab.xzs)
+    xzs = tab(frame.frame).xzs
 
-    lowbit = T(1)
-    ibig = _div(T,i-1)+1
-    ismall = _mod(T,i-1)
-    ismallm = lowbit<<(ismall)
+    lowbit, ibig, ismall, ismallm = get_bitmask_idxs(xzs,i)
 
     @inbounds @simd for f in eachindex(frame)
         r = rand()
         if  r < noise.px # X error
-            frame.frame.tab.xzs[ibig,f] ⊻= ismallm
+            xzs[ibig,f] ⊻= ismallm
         elseif r < noise.px+noise.pz # Z error
-            frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm
+            xzs[end÷2+ibig,f] ⊻= ismallm
         elseif r < noise.px+noise.pz+noise.py # Y error
-            frame.frame.tab.xzs[ibig,f] ⊻= ismallm
-            frame.frame.tab.xzs[end÷2+ibig,f] ⊻= ismallm
+            xzs[ibig,f] ⊻= ismallm
+            xzs[end÷2+ibig,f] ⊻= ismallm
         end
     end
     return frame
diff --git a/src/pauli_operator.jl b/src/pauli_operator.jl
index a843ef10b..2246d4d5f 100644
--- a/src/pauli_operator.jl
+++ b/src/pauli_operator.jl
@@ -87,19 +87,23 @@ macro P_str(a)
     quote _P_str($a) end
 end
 
-Base.getindex(p::PauliOperator{Tₚ,Tᵥ}, i::Int) where {Tₚ, Tᵥₑ<:Unsigned, Tᵥ<:AbstractVector{Tᵥₑ}} = ((p.xz[_div(Tᵥₑ, i-1)+1] & Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))!=0x0)::Bool, ((p.xz[end÷2+_div(Tᵥₑ,i-1)+1] & Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))!=0x0)::Bool
+function Base.getindex(p::PauliOperator{Tₚ,Tᵥ}, i::Int) where {Tₚ, Tᵥₑ<:Unsigned, Tᵥ<:AbstractVector{Tᵥₑ}}
+    _, ibig, _, ismallm = get_bitmask_idxs(p.xz,i)
+    ((p.xz[ibig] & ismallm) != 0x0)::Bool, ((p.xz[end÷2+ibig] & ismallm) != 0x0)::Bool
+end
 Base.getindex(p::PauliOperator{Tₚ,Tᵥ}, r) where {Tₚ, Tᵥₑ<:Unsigned, Tᵥ<:AbstractVector{Tᵥₑ}} = PauliOperator(p.phase[], xbit(p)[r], zbit(p)[r])
 
 function Base.setindex!(p::PauliOperator{Tₚ,Tᵥ}, (x,z)::Tuple{Bool,Bool}, i) where {Tₚ, Tᵥₑ, Tᵥ<:AbstractVector{Tᵥₑ}}
+    _, ibig, _, ismallm = get_bitmask_idxs(p.xz,i)
     if x
-        p.xz[_div(Tᵥₑ,i-1)+1] |= Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1)
+        p.xz[ibig] |= ismallm
     else
-        p.xz[_div(Tᵥₑ,i-1)+1] &= ~(Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))
+        p.xz[ibig] &= ~(ismallm)
     end
     if z
-        p.xz[end÷2+_div(Tᵥₑ,i-1)+1] |= Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1)
+        p.xz[end÷2+ibig] |= ismallm
     else
-        p.xz[end÷2+_div(Tᵥₑ,i-1)+1] &= ~(Tᵥₑ(0x1)<<_mod(Tᵥₑ,i-1))
+        p.xz[end÷2+ibig] &= ~(ismallm)
     end
     p
 end
diff --git a/src/project_trace_reset.jl b/src/project_trace_reset.jl
index 25aed09ce..4f7fa6211 100644
--- a/src/project_trace_reset.jl
+++ b/src/project_trace_reset.jl
@@ -42,17 +42,15 @@ function _generate!(pauli::PauliOperator{Tₚ,Tᵥ}, stabilizer::Stabilizer{Tabl
     xzs = tab(stabilizer).xzs
     xs = @view xzs[1:end÷2,:]
     zs = @view xzs[end÷2+1:end,:]
-    lowbit = Tₘₑ(0x1)
     zerobit = Tₘₑ(0x0)
     px,pz = xview(pauli), zview(pauli)
     used_indices = Int[]
     used = 0
     # remove Xs
     while (i=unsafe_bitfindnext_(px,1); i !== nothing) # TODO awkward notation due to https://github.com/JuliaLang/julia/issues/45499
-        jbig = _div(Tₘₑ,i-1)+1
-        jsmall = lowbit<<_mod(Tₘₑ,i-1)
-        candidate = findfirst(e->e&jsmall!=zerobit, # TODO some form of reinterpret might be faster than equality check
-                              xs[jbig,used+1:end])
+        _, ibig, _, ismallm = get_bitmask_idxs(xzs,i)
+        candidate = findfirst(e->e&ismallm!=zerobit, # TODO some form of reinterpret might be faster than equality check
+                              xs[ibig,used+1:end])
         if isnothing(candidate)
             return nothing
         else
@@ -63,10 +61,9 @@ function _generate!(pauli::PauliOperator{Tₚ,Tᵥ}, stabilizer::Stabilizer{Tabl
     end
     # remove Zs
     while (i=unsafe_bitfindnext_(pz,1); i !== nothing) # TODO awkward notation due to https://github.com/JuliaLang/julia/issues/45499
-        jbig = _div(Tₘₑ,i-1)+1
-        jsmall = lowbit<<_mod(Tₘₑ,i-1)
-        candidate = findfirst(e->e&jsmall!=zerobit, # TODO some form of reinterpret might be faster than equality check
-                              zs[jbig,used+1:end])
+        _, ibig, _, ismallm = get_bitmask_idxs(xzs,i)
+        candidate = findfirst(e->e&ismallm!=zerobit, # TODO some form of reinterpret might be faster than equality check
+                              zs[ibig,used+1:end])
         if isnothing(candidate)
             return nothing
         else
diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl
index 98f0e5d61..03f5f2e69 100644
--- a/src/symbolic_cliffords.jl
+++ b/src/symbolic_cliffords.jl
@@ -415,12 +415,9 @@ LinearAlgebra.inv(op::sInvSQRTZZ) = sSQRTZZ(op.q1, op.q2)
 """Apply a Pauli Z to the `i`-th qubit of state `s`. You should use `apply!(stab,sZ(i))` instead of this."""
 function apply_single_z!(stab::AbstractStabilizer, i)
     s = tab(stab)
-    Tₘₑ = eltype(s.xzs)
-    bigi = _div(Tₘₑ,i-1)+1
-    smalli = _mod(Tₘₑ,i-1)
-    mask = Tₘₑ(0x1)<<smalli
+    _, ibig, _, ismallm = get_bitmask_idxs(s.xzs,i)
     @inbounds @simd for row in 1:size(s.xzs,2)
-        if !iszero(s.xzs[bigi,row] & mask)
+        if !iszero(s.xzs[ibig,row] & ismallm)
             s.phases[row] = (s.phases[row]+0x2)&0x3
         end
     end
@@ -430,12 +427,9 @@ end
 """Apply a Pauli X to the `i`-th qubit of state `s`. You should use `apply!(stab,sX(i))` instead of this."""
 function apply_single_x!(stab::AbstractStabilizer, i)
     s = tab(stab)
-    Tₘₑ = eltype(s.xzs)
-    bigi = _div(Tₘₑ,i-1)+1
-    smalli = _mod(Tₘₑ,i-1)
-    mask = Tₘₑ(0x1)<<smalli
+    _, ibig, _, ismallm = get_bitmask_idxs(s.xzs,i)
     @inbounds @simd for row in 1:size(s.xzs,2)
-        if !iszero(s.xzs[end÷2+bigi,row] & mask)
+        if !iszero(s.xzs[end÷2+ibig,row] & ismallm)
             s.phases[row] = (s.phases[row]+0x2)&0x3
         end
     end
@@ -445,12 +439,9 @@ end
 """Apply a Pauli Y to the `i`-th qubit of state `s`. You should use `apply!(stab,sY(i))` instead of this."""
 function apply_single_y!(stab::AbstractStabilizer, i)
     s = tab(stab)
-    Tₘₑ = eltype(s.xzs)
-    bigi = _div(Tₘₑ,i-1)+1
-    smalli = _mod(Tₘₑ,i-1)
-    mask = Tₘₑ(0x1)<<smalli
+    _, ibig, _, ismallm = get_bitmask_idxs(s.xzs,i)
     @inbounds @simd for row in 1:size(s.xzs,2)
-        if !iszero((s.xzs[bigi,row] & mask) ⊻ (s.xzs[end÷2+bigi,row] & mask))
+        if !iszero((s.xzs[ibig,row] & ismallm) ⊻ (s.xzs[end÷2+ibig,row] & ismallm))
             s.phases[row] = (s.phases[row]+0x2)&0x3
         end
     end

From c3546cd48ba1ea2876687b432cbba48ba101d5fc Mon Sep 17 00:00:00 2001
From: Feroz Ahmad <feroz.ahmad.email@gmail.com>
Date: Sat, 2 Nov 2024 20:10:02 +0500
Subject: [PATCH 11/15] implement sSQRTXX, sInvSQRTXX, sSQRTYY, sInvSQRTYY
 (#368)

---
 src/QuantumClifford.jl    |  2 +-
 src/symbolic_cliffords.jl | 10 ++++++++++
 test/test_symcliff.jl     |  5 +++++
 3 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 6b4945b6b..36901ed23 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -52,7 +52,7 @@ export
     sCNOT, sCPHASE, sSWAP,
     sXCX, sXCY, sXCZ, sYCX, sYCY, sYCZ, sZCX, sZCY, sZCZ,
     sZCrY, sInvZCrY, sSWAPCX, sInvSWAPCX, sCZSWAP, sCXSWAP, sISWAP, sInvISWAP,
-    sSQRTZZ, sInvSQRTZZ,
+    sSQRTZZ, sInvSQRTZZ, sSQRTXX, sInvSQRTXX, sSQRTYY, sInvSQRTYY,
     # Misc Ops
     SparseGate,
     sMX, sMY, sMZ, PauliMeasurement, Reset, sMRX, sMRY, sMRZ,
diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl
index 03f5f2e69..fcff4cb38 100644
--- a/src/symbolic_cliffords.jl
+++ b/src/symbolic_cliffords.jl
@@ -339,6 +339,12 @@ end
 @qubitop2 SQRTZZ    (x1       , x1⊻x2⊻z1 , x2       , x1⊻z2⊻x2 , ~iszero((x1 & z1 & ~x2) | (~x1 & x2 & z2)))
 @qubitop2 InvSQRTZZ (x1       , x1⊻x2⊻z1 , x2       , x1⊻z2⊻x2 , ~iszero((x1 &~z1 & ~x2) | (~x1 & x2 &~z2)))
 
+@qubitop2 SQRTXX    (z1⊻z2⊻x1, z1      , z1⊻x2⊻z2, z2      , ~iszero((~x1 & z1 &~z2) | (~z1 &~x2 & z2)))
+@qubitop2 InvSQRTXX (z1⊻z2⊻x1, z1      , z1⊻x2⊻z2, z2      , ~iszero(( x1 & z1 &~z2) | (~z1 & x2 & z2)))
+
+@qubitop2 SQRTYY    (z1⊻x2⊻z2, x1⊻z2⊻x2, x1⊻z1⊻z2, x1⊻x2⊻z1, ~iszero((~x1 &~z1 & x2 &~z2) | ( x1 &~z1 &~x2 &~z2) | ( x1 &~z1 & x2 & z2) | ( x1 & z1 & x2 &~z2)))
+@qubitop2 InvSQRTYY (z1⊻x2⊻z2, x1⊻z2⊻x2, x1⊻z1⊻z2, x1⊻x2⊻z1, ~iszero(( x1 & z1 &~x2 & z2) | (~x1 & z1 & x2 & z2) | (~x1 & z1 &~x2 &~z2) | (~x1 &~z1 &~x2 & z2)))
+
 #=
 To get the boolean formulas for the phase, it is easiest to first write down the truth table for the phase:
 for i in 0:15
@@ -405,6 +411,10 @@ LinearAlgebra.inv(op::sISWAP)     = sInvISWAP(op.q1, op.q2)
 LinearAlgebra.inv(op::sInvISWAP)  = sISWAP(op.q1, op.q2)
 LinearAlgebra.inv(op::sSQRTZZ)    = sInvSQRTZZ(op.q1, op.q2)
 LinearAlgebra.inv(op::sInvSQRTZZ) = sSQRTZZ(op.q1, op.q2)
+LinearAlgebra.inv(op::sSQRTXX)    = sInvSQRTXX(op.q1, op.q2)
+LinearAlgebra.inv(op::sInvSQRTXX) = sSQRTXX(op.q1, op.q2)
+LinearAlgebra.inv(op::sSQRTYY)    = sInvSQRTYY(op.q1, op.q2)
+LinearAlgebra.inv(op::sInvSQRTYY) = sSQRTYY(op.q1, op.q2)
 
 ##############################
 # Functions that perform direct application of common operators without needing an operator instance
diff --git a/test/test_symcliff.jl b/test/test_symcliff.jl
index 40ed505ac..6bbad80d4 100644
--- a/test/test_symcliff.jl
+++ b/test/test_symcliff.jl
@@ -76,6 +76,7 @@
             @test CliffordOperator(inv(random_op), i) == inv(CliffordOperator(random_op, i))
             @test CliffordOperator(inv(SingleQubitOperator(random_op)), i) == inv(CliffordOperator(random_op, i))
         end
+    end
 
     @testset "Consistency checks with Stim" begin
        # see https://github.com/quantumlib/Stim/blob/main/doc/gates.md
@@ -113,5 +114,9 @@
         @test CliffordOperator(sInvISWAP)  == C"-ZY -YZ IZ ZI"
         @test CliffordOperator(sSQRTZZ)    == C"YZ ZY ZI IZ"
         @test CliffordOperator(sInvSQRTZZ) == C"-YZ -ZY ZI IZ"
+        @test CliffordOperator(sSQRTXX)    == C"XI IX -YX -XY"
+        @test CliffordOperator(sInvSQRTXX) == C"XI IX YX XY"
+        @test CliffordOperator(sSQRTYY)    == C"-ZY -YZ XY YX"
+        @test CliffordOperator(sInvSQRTYY) == C"ZY YZ -XY -YX"
     end
 end

From 56acb6ebf9b3aa75d79a8ddc0d03cd5bb83aecb7 Mon Sep 17 00:00:00 2001
From: Feroz Ahmad <feroz.ahmad.email@gmail.com>
Date: Sun, 3 Nov 2024 05:42:36 +0500
Subject: [PATCH 12/15] fix #191: too restrictive typeassert for
 MixedDestabilizer (#366)

Co-authored-by: Stefan Krastanov <stefan@krastanov.org>
---
 src/QuantumClifford.jl |  4 ++--
 test/test_stabs.jl     | 14 ++++++++++++++
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 36901ed23..6b90af7a8 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -571,8 +571,8 @@ function MixedDestabilizer(stab::Stabilizer{T}; undoperm=true, reportperm=false)
         t[n+r+s+1:end] = sZ # The other logical set in the tableau
     end
     if undoperm
-        t = t[:,invperm(permx[permz])]::T
-        return MixedDestabilizer(t, r+s)::MixedDestabilizer{T}
+        t = t[:,invperm(permx[permz])]
+        return MixedDestabilizer(t, r+s)
     end
     if reportperm
         return (MixedDestabilizer(t, r+s)::MixedDestabilizer{T}, r, permx, permz)
diff --git a/test/test_stabs.jl b/test/test_stabs.jl
index b479934c2..d6b9041de 100644
--- a/test/test_stabs.jl
+++ b/test/test_stabs.jl
@@ -1,4 +1,5 @@
 @testitem "Stabilizers" begin
+    using QuantumClifford
     using QuantumClifford: stab_looks_good, destab_looks_good, mixed_stab_looks_good, mixed_destab_looks_good
     test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
     @testset "Pure and Mixed state initialization" begin
@@ -107,4 +108,17 @@
         @test hcat(copy(tab(s1)), copy(tab(s2))) == T"-YZZY XXYX"
         @test hcat(copy(tab(s1)), copy(tab(s2)), copy(tab(s1)), copy(tab(s2))) == T"YZZYYZZY XXYXXXYX"
     end
+
+    @testset "MixedDestabilizer over subarrays (#191)" begin
+        # Case 1: QuantumClifford.Tableau{Vector{UInt8}, Matrix{UInt64}}
+        n = 6
+        stab = random_stabilizer(n)
+        regular_arr = MixedDestabilizer(stab; undoperm=true)
+        @test isa(regular_arr, MixedDestabilizer)
+        # Case 2: Tableau{SubArray{...}, SubArray{...}, Tuple{Base.Slice{...}}}
+        stab = random_stabilizer(n)
+        substab = @view stab[3:n]
+        md_via_subarr = MixedDestabilizer(substab; undoperm=true)
+        @test isa(md_via_subarr, MixedDestabilizer)
+    end
 end

From 406e5eaf23f3c461906b7e58f72230072571568d Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Sun, 3 Nov 2024 10:27:20 -0500
Subject: [PATCH 13/15] update changelog

---
 CHANGELOG.md | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e324370b0..7bf3a570d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,9 +5,11 @@
 
 # News
 
-## v0.9.14 - 2024-11-02
+## v0.9.14 - 2024-11-03
 
 - **(fix)** `affectedqubits()` on `sMX`, `sMY`, and `sMR*`
+- **(fix)** restrictive type-assert in `MixedDestabilizer` failing on views of tableaux
+- Implementing additional named two-qubit gates: `sSQRTXX, sInvSQRTXX, sSQRTYY, sInvSQRTYY`
 
 ## v0.9.13 - 2024-10-30
 

From 553133622ebe66c3b35675ada9a739027e3f9367 Mon Sep 17 00:00:00 2001
From: Feroz Ahmad <feroz.ahmad.email@gmail.com>
Date: Tue, 5 Nov 2024 09:05:05 +0500
Subject: [PATCH 14/15] add `check_repr_commutation_relation` to test the CSS
 orthogonality condition for 2BGA's Group Algebra with a General Group `G`
 (#403)

Co-authored-by: Stefan Krastanov <github.acc@krastanov.org>
---
 .../QuantumCliffordHeckeExt.jl                    |  5 +++--
 ext/QuantumCliffordHeckeExt/util.jl               | 15 +++++++++++++++
 src/ecc/codes/util.jl                             |  3 +++
 test/test_ecc_base.jl                             |  2 ++
 4 files changed, 23 insertions(+), 2 deletions(-)
 create mode 100644 ext/QuantumCliffordHeckeExt/util.jl

diff --git a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
index 29e9de8ce..b8508ff66 100644
--- a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
+++ b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
@@ -5,14 +5,15 @@ using DocStringExtensions
 import QuantumClifford, LinearAlgebra
 import Hecke: Group, GroupElem, AdditiveGroup, AdditiveGroupElem,
     GroupAlgebra, GroupAlgebraElem, FqFieldElem, representation_matrix, dim, base_ring,
-    multiplication_table, coefficients, abelian_group, group_algebra
+    multiplication_table, coefficients, abelian_group, group_algebra, rand
 import Nemo
 import Nemo: characteristic, matrix_repr, GF, ZZ, lift
 
 import QuantumClifford.ECC: AbstractECC, CSS, ClassicalCode,
     hgp, code_k, code_n, code_s, iscss, parity_checks, parity_checks_x, parity_checks_z, parity_checks_xz,
-    two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes
+    two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes, check_repr_commutation_relation
 
+include("util.jl")
 include("types.jl")
 include("lifted.jl")
 include("lifted_product.jl")
diff --git a/ext/QuantumCliffordHeckeExt/util.jl b/ext/QuantumCliffordHeckeExt/util.jl
new file mode 100644
index 000000000..d8fc65e85
--- /dev/null
+++ b/ext/QuantumCliffordHeckeExt/util.jl
@@ -0,0 +1,15 @@
+"""
+Checks the commutation relation between the left and right representation matrices
+for two randomly-sampled elements `a` and `b` in the group algebra `ℱ[G]` with a general group `G`.
+It verifies the commutation relation that states, `L(a)·R(b) = R(b)·L(a)`. This
+property shows that matrices from the left and right representation sets commute
+with each other, which is an important property related to the CSS orthogonality
+condition.
+"""
+function check_repr_commutation_relation(GA::GroupAlgebra)
+    a, b = rand(GA), rand(GA)
+    # Check commutation relation: L(a)R(b) = R(b)L(a)
+    L_a = representation_matrix(a)
+    R_b = representation_matrix(b, :right)
+    return L_a * R_b == R_b * L_a
+end
diff --git a/src/ecc/codes/util.jl b/src/ecc/codes/util.jl
index 1063785f9..b96a75120 100644
--- a/src/ecc/codes/util.jl
+++ b/src/ecc/codes/util.jl
@@ -6,3 +6,6 @@ function hgp(h₁,h₂)
     hz = hcat(kron(LinearAlgebra.I(n₁), h₂), kron(h₁', LinearAlgebra.I(r₂)))
     hx, hz
 end
+
+"""Implemented in a package extension with Hecke."""
+function check_repr_commutation_relation end
diff --git a/test/test_ecc_base.jl b/test/test_ecc_base.jl
index 9f4a5ec9e..26b00c8f3 100644
--- a/test/test_ecc_base.jl
+++ b/test/test_ecc_base.jl
@@ -1,6 +1,7 @@
 using Test
 using QuantumClifford
 using QuantumClifford.ECC
+using QuantumClifford.ECC: check_repr_commutation_relation
 using InteractiveUtils
 
 import Nemo: GF
@@ -46,6 +47,7 @@ other_lifted_product_codes = []
 # [[882, 24, d≤24]] code from (B1) in Appendix B of [panteleev2021degenerate](@cite)
 l = 63
 GA = group_algebra(GF(2), abelian_group(l))
+@test check_repr_commutation_relation(GA) # TODO use this check more pervasively throughout the test suite
 A = zeros(GA, 7, 7)
 x = gens(GA)[]
 A[LinearAlgebra.diagind(A)] .= x^27

From 40924b804c128b3aa2186cd647b7e3252bbd0180 Mon Sep 17 00:00:00 2001
From: Feroz Ahmad <feroz.ahmad.email@gmail.com>
Date: Tue, 5 Nov 2024 09:07:30 +0500
Subject: [PATCH 15/15] =?UTF-8?q?add=20tests=20and=20examples=20of=20=20C?=
 =?UTF-8?q?=E2=82=98=20=C3=97=20C=E2=82=82=202BGA=20codes=20(#392)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Stefan Krastanov <github.acc@krastanov.org>
---
 docs/src/references.bib                       |  11 ++
 ext/QuantumCliffordHeckeExt/lifted_product.jl |  28 ++++-
 test/test_ecc_2bga.jl                         | 119 ++++++++++++++++++
 3 files changed, 154 insertions(+), 4 deletions(-)
 create mode 100644 test/test_ecc_2bga.jl

diff --git a/docs/src/references.bib b/docs/src/references.bib
index 281ace9c1..d82fc7440 100644
--- a/docs/src/references.bib
+++ b/docs/src/references.bib
@@ -512,3 +512,14 @@ @article{anderson2014fault
   year={2014},
   publisher={APS}
 }
+
+@article{lin2024quantum,
+  title={Quantum two-block group algebra codes},
+  author={Lin, Hsiang-Ku and Pryadko, Leonid P},
+  journal={Physical Review A},
+  volume={109},
+  number={2},
+  pages={022407},
+  year={2024},
+  publisher={APS}
+}
diff --git a/ext/QuantumCliffordHeckeExt/lifted_product.jl b/ext/QuantumCliffordHeckeExt/lifted_product.jl
index ef09eddbb..21cd96481 100644
--- a/ext/QuantumCliffordHeckeExt/lifted_product.jl
+++ b/ext/QuantumCliffordHeckeExt/lifted_product.jl
@@ -153,12 +153,32 @@ code_s(c::LPCode) = size(c.repr(zero(c.GA)), 1) * (size(c.A, 1) * size(c.B, 1) +
 Two-block group algebra (2GBA) codes, which are a special case of lifted product codes
 from two group algebra elements `a` and `b`, used as `1x1` base matrices.
 
+Here is an example of a [[56, 28, 2]] 2BGA code from Table 2 of [lin2024quantum](@cite)
+with direct product of `C₄ x C₂`.
+
+```jldoctest
+julia> import Hecke: group_algebra, GF, abelian_group, gens;
+
+julia> GA = group_algebra(GF(2), abelian_group([14,2]));
+
+julia> x = gens(GA)[1];
+
+julia> s = gens(GA)[2];
+
+julia> A = 1 + x^7
+
+julia> B = 1 + x^7 + s + x^8 + s*x^7 + x
+
+julia> c = two_block_group_algebra_codes(A,B);
+
+julia> code_n(c), code_k(c)
+(56, 28)
+```
+
 See also: [`LPCode`](@ref), [`generalized_bicycle_codes`](@ref), [`bicycle_codes`](@ref)
-""" # TODO doctest example
+"""
 function two_block_group_algebra_codes(a::GroupAlgebraElem, b::GroupAlgebraElem)
-    A = reshape([a], (1, 1))
-    B = reshape([b], (1, 1))
-    LPCode(A, B)
+    LPCode([a;;], [b;;])
 end
 
 """
diff --git a/test/test_ecc_2bga.jl b/test/test_ecc_2bga.jl
new file mode 100644
index 000000000..e284abfab
--- /dev/null
+++ b/test/test_ecc_2bga.jl
@@ -0,0 +1,119 @@
+@testitem "ECC 2BGA" begin
+    using Hecke
+    using Hecke: group_algebra, GF, abelian_group, gens
+    using QuantumClifford.ECC: LPCode, code_k, code_n
+
+    @testset "Reproduce Table 2 lin2024quantum" begin # TODO these tests should probably just use the `two_block_group_algebra_codes` function as that would make them much shorter and simpler
+        # codes taken from Table 2 of [lin2024quantum](@cite)
+
+        # m = 4
+        GA = group_algebra(GF(2), abelian_group([4,2]))
+        x = gens(GA)[1]
+        s = gens(GA)[2]
+        A = [1 + x;;]
+        B = [1 + x + s + x^2 + s*x + s*x^3;;]
+        c = LPCode(A,B)
+        # [[16, 2, 4]] 2BGA code
+        @test code_n(c) == 16 && code_k(c) == 2
+        A = [1 + x;;]
+        B = [1 + x + s + x^2 + s*x + x^3;;]
+        c = LPCode(A,B)
+        # [[16, 4, 4]] 2BGA code
+        @test code_n(c) == 16 && code_k(c) == 4
+        A = [1 + s;;]
+        B = [1 + x + s + x^2 + s*x + x^2;;]
+        c = LPCode(A,B)
+        # [[16, 8, 2]] 2BGA code
+        @test code_n(c) == 16 && code_k(c) == 8
+
+        # m = 6
+        GA = group_algebra(GF(2), abelian_group([6,2]))
+        x = gens(GA)[1]
+        s = gens(GA)[2]
+        A = [1 + x;;]
+        B = [1 + x^3 + s + x^4 + x^2 + s*x;;]
+        c = LPCode(A,B)
+        # [[24, 4, 5]] 2BGA code
+        @test code_n(c) == 24 && code_k(c) == 4
+        A = [1 + x^3;;]
+        B = [1 + x^3 + s + x^4 + s*x^3 + x;;]
+        c = LPCode(A,B)
+        # [[24, 12, 2]] 2BGA code
+        @test code_n(c) == 24 && code_k(c) == 12
+
+        # m = 8
+        GA = group_algebra(GF(2), abelian_group([8,2]))
+        x = gens(GA)[1]
+        s = gens(GA)[2]
+        A = [1 + x^6;;]
+        B = [1 + s*x^7 + s*x^4 + x^6 + s*x^5 + s*x^2;;]
+        c = LPCode(A,B)
+        # [[32, 8, 4]] 2BGA code
+        @test code_n(c) == 32 && code_k(c) == 8
+        A = [1 + s*x^4;;]
+        B = [1 + s*x^7 + s*x^4 + x^6 + x^3 + s*x^2;;]
+        c = LPCode(A,B)
+        # [[32, 16, 2]] 2BGA code
+        @test code_n(c) == 32 && code_k(c) == 16
+
+        # m = 10
+        GA = group_algebra(GF(2), abelian_group([10,2]))
+        x = gens(GA)[1]
+        s = gens(GA)[2]
+        A = [1 + x;;]
+        B = [1 + x^5 + x^6 + s*x^6 + x^7 + s*x^3;;]
+        c = LPCode(A,B)
+        # [[40, 4, 8]] 2BGA code
+        @test code_n(c) == 40 && code_k(c) == 4
+        A = [1 + x^6;;]
+        B = [1 + x^5 + s + x^6 + x + s*x^2;;]
+        c = LPCode(A,B)
+        # [[40, 8, 5]] 2BGA code
+        @test code_n(c) == 40 && code_k(c) == 8
+        A = [1 + x^5;;]
+        B = [1 + x^5 + s + x^6 + s*x^5 + x;;]
+        c = LPCode(A,B)
+        # [[40, 20, 2]] 2BGA code
+        @test code_n(c) == 40 && code_k(c) == 20
+
+        # m = 12
+        GA = group_algebra(GF(2), abelian_group([12,2]))
+        x = gens(GA)[1]
+        s = gens(GA)[2]
+        A = [1 + s*x^10;;]
+        B = [1 + x^3 + s*x^6 + x^4 + x^7 + x^8;;]
+        c = LPCode(A,B)
+        # [[48, 8, 6]] 2BGA code
+        @test code_n(c) == 48 && code_k(c) == 8
+        A = [1 + x^3;;]
+        B = [1 + x^3 + s*x^6 + x^4 + s*x^9 + x^7;;]
+        c = LPCode(A,B)
+        # [[48, 12, 4]] 2BGA code
+        @test code_n(c) == 48 && code_k(c) == 12
+        A = [1 + x^4;;]
+        B = [1 + x^3 + s*x^6 + x^4 + x^7 + s*x^10;;]
+        c = LPCode(A,B)
+        # [[48, 16, 3]] 2BGA code
+        @test code_n(c) == 48 && code_k(c) == 16
+        A = [1 + s*x^6;;]
+        B = [1 + x^3 + s*x^6 + x^4 + s*x^9 + s*x^10;;]
+        c = LPCode(A,B)
+        # [[48, 24, 2]] 2BGA code
+        @test code_n(c) == 48 && code_k(c) == 24
+
+        # m = 14
+        GA = group_algebra(GF(2), abelian_group([14,2]))
+        x = gens(GA)[1]
+        s = gens(GA)[2]
+        A = [1 + x^8;;]
+        B = [1 + x^7 + s + x^8 + x^9 + s*x^4;;]
+        c = LPCode(A,B)
+        # [[56, 8, 7]] 2BGA code
+        @test code_n(c) == 56 && code_k(c) == 8
+        A = [1 + x^7;;]
+        B = [1 + x^7 + s + x^8 + s*x^7 + x;;]
+        c = LPCode(A,B)
+        # [[56, 28, 2]] 2BGA code
+        @test code_n(c) == 56 && code_k(c) == 28
+    end
+end