diff --git a/pytket/conanfile.py b/pytket/conanfile.py
index c50c6c9909..bc553d5788 100644
--- a/pytket/conanfile.py
+++ b/pytket/conanfile.py
@@ -38,7 +38,7 @@ def requirements(self):
         self.requires("pybind11_json/0.2.14")
         self.requires("symengine/0.13.0")
         self.requires("tkassert/0.3.4@tket/stable")
-        self.requires("tket/1.3.53@tket/stable")
+        self.requires("tket/1.3.54@tket/stable")
         self.requires("tklog/0.3.3@tket/stable")
         self.requires("tkrng/0.3.3@tket/stable")
         self.requires("tktokenswap/0.3.9@tket/stable")
diff --git a/tket/CMakeLists.txt b/tket/CMakeLists.txt
index 58fa4797f8..683c351eef 100644
--- a/tket/CMakeLists.txt
+++ b/tket/CMakeLists.txt
@@ -169,9 +169,11 @@ target_sources(tket
         src/Circuit/SubcircuitFinder.cpp
         src/Circuit/ThreeQubitConversion.cpp
         src/Circuit/ToffoliBox.cpp
+        src/Clifford/APState.cpp
         src/Clifford/ChoiMixTableau.cpp
         src/Clifford/SymplecticTableau.cpp
         src/Clifford/UnitaryTableau.cpp
+        src/Converters/APStateConverters.cpp
         src/Converters/ChoiMixTableauConverters.cpp
         src/Converters/Gauss.cpp
         src/Converters/PauliGraphConverters.cpp
@@ -324,6 +326,7 @@ target_sources(tket
         include/tket/Circuit/StatePreparation.hpp
         include/tket/Circuit/ThreeQubitConversion.hpp
         include/tket/Circuit/ToffoliBox.hpp
+        include/tket/Clifford/APState.hpp
         include/tket/Clifford/ChoiMixTableau.hpp
         include/tket/Clifford/SymplecticTableau.hpp
         include/tket/Clifford/UnitaryTableau.hpp
diff --git a/tket/conanfile.py b/tket/conanfile.py
index 2602d319f1..dcc80f7bb7 100644
--- a/tket/conanfile.py
+++ b/tket/conanfile.py
@@ -23,7 +23,7 @@
 
 class TketConan(ConanFile):
     name = "tket"
-    version = "1.3.53"
+    version = "1.3.54"
     package_type = "library"
     license = "Apache 2"
     homepage = "https://github.com/CQCL/tket"
diff --git a/tket/include/tket/Clifford/APState.hpp b/tket/include/tket/Clifford/APState.hpp
new file mode 100644
index 0000000000..3e49111205
--- /dev/null
+++ b/tket/include/tket/Clifford/APState.hpp
@@ -0,0 +1,369 @@
+// Copyright 2019-2024 Cambridge Quantum Computing
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+#include "tket/OpType/OpType.hpp"
+#include "tket/Utils/MatrixAnalysis.hpp"
+
+namespace tket {
+
+/**
+ * APState class, giving a unique form for a (possibly mixed) Clifford state,
+ * through which we can track global phase when appropriate.
+ *
+ * The "affine with phases" form of a Clifford state from ZX calculus (see
+ * Kissinger & van de Wetering, "Picturing Quantum Software") represents n-qubit
+ * Clifford states uniquely with the following data:
+ * - A binary (n,n) matrix A.
+ * - A binary (n) vector B.
+ * - A symmetric, zero-diagonal, binary (n,n) matrix E.
+ * - A n vector P of integers mod 4 describing S gates.
+ *
+ * This gives a canonical statevector (up to a normalisation scalar):
+ * \f[\sum_{x, Ax=B} i^{\Phi(x)} |x>\f]
+ *
+ * This canonical statevector fixes a reference phase from which we can track
+ * global phase with an additional parameter.
+ *
+ * We can generalise this to mixed qubit states by adding a binary (n,n) matrix
+ * C, with which we now have the canonical density matrix (up to a normalisation
+ * scalar):
+ * \f[\sum_{x_1, x_2; Ax_1 = B = Ax_2, Cx_1 = Cx_2}
+ *         i^{x_1^T E x_1 + P^T x_1} |x_1><x_2| (-i)^{x_2^T E x_2 + P^T x_2}\f]
+ * where the inner products are calculated in Z4.
+ *
+ * We can encode this via the ZX-calculus as:
+ * - A green spider for each qubit q, with phase given by P(q)pi/2.
+ * - The green spiders are connected via Hadamard edges according to E.
+ * - For each i, a red spider with phase B(i)pi, connected to the green spiders
+ *   according to the row A(i, -).
+ * - For each j, a discard connected to a red spider, connected to the green
+ *   spiders according to the row C(j, -).
+ *
+ * Several operations on this data leave the state unchanged:
+ * - We can freely add rows simultaneously in A and B whilst preserving the
+ *   state, and similarly combine rows in C or add rows from A to C.
+ * - If a green spider q is connected to precisely one red spider (a single 1
+ *   appears in the column A(-, q), and nothing in column C(-, q)),
+ *   transformations exist that allow us to set P(q), E(q, -), and E(-, q) to
+ *   zero, by some sequence of bialgebra and local complementation rules.
+ * - We can perform a local complementation around any discard j - concerning
+ *   its neighbours (those q with C(j, q) = 1), we may increment each P(q) by 1
+ *   and flip each E(q, q').
+ * However, whilst a normal form exists (generalising "reduced AP-form"),
+ * attempting to maintain that form massively increases the number of cases we
+ * need to consider for each gate application. We instead just offer a method
+ * that reduces a given APState to the normal form which can be used for
+ * comparison checks.
+ *
+ * To mimic all behaviour that ChoiMixTableau supports, we consider the gate set
+ * {CZ, S, V, Init, PostSelect, Collapse}, where V is the only difficult one to
+ * derive (always comes down to a sequence of local complementations in the ZX
+ * picture).
+ */
+class APState {
+ public:
+  /**
+   * A binary (n,n) matrix used in describing the subspace of computational
+   * basis states in the support of the state.
+   */
+  MatrixXb A;
+
+  /**
+   * A binary n vector used in describing the subspace of computaitonal basis
+   * states in the support of the state.
+   */
+  VectorXb B;
+
+  /**
+   * A binary (n,n) matrix describing the incoherent constraints, relating the
+   * ket and bra sides of the density matrix representation of the state.
+   */
+  MatrixXb C;
+
+  /**
+   * A symmetric, zero diagonal matrix whose entries indicate CZs between
+   * qubits.
+   */
+  MatrixXb E;
+
+  /**
+   * A vector indicating S^{P(i)} on qubit i.
+   */
+  Eigen::VectorXi P;
+
+  /**
+   * The global phase term (in half-turns).
+   */
+  Expr phase;
+
+  /**
+   * Constructs a state in AP form from given data.
+   */
+  APState(
+      const MatrixXb& A_, const VectorXb& B_, const MatrixXb& C_,
+      const MatrixXb& E_, const Eigen::VectorXi& P_, const Expr& phase_);
+
+  /**
+   * Constructs the state |0>^{x n_qubits} in AP form.
+   */
+  APState(unsigned n_qubits);
+
+  /**
+   * Constructs the state in AP form from a given statevector.
+   */
+  APState(const Eigen::VectorXcd& sv);
+
+  /**
+   * Constructs the state in AP form from a given density matrix.
+   */
+  APState(const Eigen::MatrixXcd& dm);
+
+  /**
+   * Other required constructors
+   */
+  APState(const APState& other) = default;
+  APState(APState&& other) = default;
+  APState& operator=(const APState& other) = default;
+  APState& operator=(APState& other) = default;
+
+  /**
+   * Verifies the internal correctness of the data structure. Throws an
+   * exception if the data structure is invalid.
+   */
+  void verify() const;
+
+  /**
+   * Calculates the statevector of the state.
+   *
+   * Throws an exception if C is non-zero (i.e. it represents a mixed state).
+   */
+  Eigen::VectorXcd to_statevector() const;
+
+  /**
+   * Calculates the density matrix of the state.
+   */
+  Eigen::MatrixXcd to_density_matrix() const;
+
+  /**
+   * Applies a CZ gate to the state. O(1).
+   */
+  void apply_CZ(unsigned ctrl, unsigned trgt);
+
+  /**
+   * Applies an S gate to the state. O(1).
+   */
+  void apply_S(unsigned q);
+
+  /**
+   * Applies a V gate to the state. O(n^2) wrt number of qubits in the state.
+   */
+  void apply_V(unsigned q);
+
+  /**
+   * Applies an X gate to the state, faster than applying V twice. O(n) wrt
+   * number of qubits in the state.
+   */
+  void apply_X(unsigned q);
+
+  /**
+   * Applies an unparameterised Clifford gate to the state on the chosen qubits.
+   * O(n^2) wrt number of qubits in the state.
+   */
+  void apply_gate(OpType type, const std::vector<unsigned>& qbs);
+
+  /**
+   * Initialises a new qubit in the |0> state. Returns the index of the new
+   * qubit. O(n) wrt number of qubits in the state.
+   */
+  unsigned init_qubit();
+
+  /**
+   * Post-selects the chosen qubit to <0|, removing it from the state. Moves the
+   * final qubit into the place of q, returning its old index. O(n^3) wrt
+   * number of qubits in the state.
+   */
+  unsigned post_select(unsigned q);
+
+  /**
+   * Collapses the given qubit in the Z basis. O(n^3) wrt number of qubits in
+   * the state.
+   */
+  void collapse_qubit(unsigned q);
+
+  /**
+   * Discards the given qubit, removing it from the state. Moves the final qubit
+   * into the place of q, returning its old index. O(n^3) wrt number of qubits
+   * in the state.
+   */
+  unsigned discard_qubit(unsigned q);
+
+  /**
+   * Manipulates the state into the following normal form:
+   * - A is in reduced row-echelon form. We refer to the leading columns of each
+   * row of A as "leading qubits".
+   * - C is in reduced row-echelon form, and furthermore is zero in any column
+   * of a leading qubit. We refer to the leading columns of each row of C as
+   * "mixed qubits".
+   * - Each entry of E is either between a mixed qubit and a free qubit, or
+   * between two free qubits.
+   * - For each leading or mixed qubit, the index of P is zero.
+   *
+   * Takes time O(n^4) wrt number of qubits in the state.
+   */
+  void normal_form();
+
+  bool operator==(const APState& other) const;
+};
+
+/**
+ * A wrapper for APState to provide the same interface as ChoiMixTableau wrt
+ * qubit indexing and distinction between input vs output.
+ *
+ * When applying gates, the methods automatically transpose anything occurring
+ * over the input subspace according to the CJ-isomorphism.
+ */
+class ChoiAPState {
+ public:
+  enum class TableauSegment { Input, Output };
+  typedef std::pair<Qubit, TableauSegment> col_key_t;
+  typedef boost::bimap<col_key_t, unsigned> tableau_col_index_t;
+
+  /**
+   * The internal APState.
+   */
+  APState ap_;
+
+  /**
+   * Map between column indices and the corresponding qubit ID and type.
+   */
+  tableau_col_index_t col_index_;
+
+  /**
+   * Construct the identity unitary over n qubits (given default qubit names).
+   */
+  explicit ChoiAPState(unsigned n);
+  /**
+   * Construct the identity unitary over specific qubits.
+   */
+  explicit ChoiAPState(const qubit_vector_t& qbs);
+  /**
+   * Construct the APState from the underlying binary matrices.
+   * Qubits are given default names and mapped such that the first columns are
+   * inputs and the last columns are outputs.
+   */
+  explicit ChoiAPState(
+      const MatrixXb& A, const VectorXb& b, const MatrixXb& C,
+      const MatrixXb& E, const Eigen::VectorXi& P, const Expr& phase = 0,
+      unsigned n_ins = 0);
+  /**
+   * Other required constructors
+   */
+  ChoiAPState(const ChoiAPState& other) = default;
+  ChoiAPState(ChoiAPState&& other) = default;
+  ChoiAPState& operator=(const ChoiAPState& other) = default;
+  ChoiAPState& operator=(ChoiAPState& other) = default;
+
+  /**
+   * Get the total number of qubits/boundaries (sum of inputs and outputs) of
+   * the process.
+   */
+  unsigned get_n_boundaries() const;
+
+  /**
+   * Get the number of boundaries representing inputs to the process.
+   */
+  unsigned get_n_inputs() const;
+
+  /**
+   * Get the number of boundaries representing outputs to the process.
+   */
+  unsigned get_n_outputs() const;
+
+  /**
+   * Get all qubit names present in the input segment.
+   */
+  qubit_vector_t input_qubits() const;
+
+  /**
+   * Get all qubit names present in the output segment.
+   */
+  qubit_vector_t output_qubits() const;
+
+  /**
+   * Transforms the state according to consuming a Clifford gate at either end
+   * of the circuit. Applying a gate to the inputs really consumes its transpose
+   * via the CJ isomorphism. This method handles the transposition internally by
+   * reversing the order of basic self-transpose gates. For multi-qubit gates,
+   * the qubits applied must either be all inputs or all outputs.
+   */
+  void apply_gate(
+      OpType type, const qubit_vector_t& qbs,
+      TableauSegment seg = TableauSegment::Output);
+
+  /**
+   * Initialises a new qubit in the |0> state (or <0| for inputs). O(n) wrt
+   * number of qubits in the state.
+   */
+  void init_qubit(const Qubit& qb, TableauSegment seg = TableauSegment::Output);
+
+  /**
+   * Post-selects the chosen qubit to <0| (or applies |0> on inputs), removing
+   * it from the state. This DOES NOT CHECK whether or not the post-selection
+   * would be successful (e.g. we can post-select <0| on a |1> state). O(n^3)
+   * wrt number of qubits in the state.
+   */
+  void post_select(
+      const Qubit& qb, TableauSegment seg = TableauSegment::Output);
+
+  /**
+   * Discards the given qubit, removing it from the state. O(n^3) wrt number of
+   * qubits in the state.
+   */
+  void discard_qubit(
+      const Qubit& qb, TableauSegment seg = TableauSegment::Output);
+
+  /**
+   * Collapses the given qubit in the Z basis. O(n^3) wrt number of qubits in
+   * the state.
+   */
+  void collapse_qubit(
+      const Qubit& qb, TableauSegment seg = TableauSegment::Output);
+
+  /**
+   * Permutes columns into a canonical order: chosen segment first, subordered
+   * in ILO.
+   */
+  void canonical_column_order(TableauSegment first = TableauSegment::Input);
+
+  /**
+   * Reduces the underlying APState to its normal form wrt the current ordering
+   * of qubits in its columns.
+   */
+  void normal_form();
+
+  /**
+   * Renames qubits.
+   */
+  void rename_qubits(const qubit_map_t& qmap, TableauSegment seg);
+
+  bool operator==(const ChoiAPState& other) const;
+};
+
+JSON_DECL(APState)
+JSON_DECL(ChoiAPState::TableauSegment)
+JSON_DECL(ChoiAPState)
+
+}  // namespace tket
\ No newline at end of file
diff --git a/tket/include/tket/Converters/Converters.hpp b/tket/include/tket/Converters/Converters.hpp
index a3de8c9146..2ccdd3cc31 100644
--- a/tket/include/tket/Converters/Converters.hpp
+++ b/tket/include/tket/Converters/Converters.hpp
@@ -15,6 +15,7 @@
 #pragma once
 
 #include "tket/Circuit/Circuit.hpp"
+#include "tket/Clifford/APState.hpp"
 #include "tket/Clifford/ChoiMixTableau.hpp"
 #include "tket/Clifford/UnitaryTableau.hpp"
 #include "tket/PauliGraph/PauliGraph.hpp"
@@ -130,6 +131,73 @@ std::pair<Circuit, qubit_map_t> cm_tableau_to_unitary_extension_circuit(
     const std::vector<Qubit> &post_names = {},
     CXConfigType cx_config = CXConfigType::Snake);
 
+/**
+ * Construct an APState for a given circuit.
+ * Assumes all input qubits are initialised in the |0> state.
+ * Will throw an exception if it contains non-Clifford gates.
+ */
+APState circuit_to_apstate(const Circuit &circ);
+
+/**
+ * Construct a circuit producing the state described by the APState.
+ * Uses the standard circuit form implied by the ZX description of reduced
+ * AP-form (layered as X-H-CX-(Collapse-CX)*-CZ-S).
+ */
+Circuit apstate_to_circuit(const APState &ap);
+
+/**
+ * Construct a ChoiAPState for a given circuit.
+ * Will incorporate qubit initialisations and discarding into the circuit.
+ * Will throw an exception if it contains non-Clifford gates.
+ */
+ChoiAPState circuit_to_choi_apstate(const Circuit &circ);
+
+/**
+ * Constructs a circuit producing the same effect as a ChoiAPState.
+ * Uses the same synthesis method as cm_tableau_to_exact_circuit.
+ */
+std::pair<Circuit, qubit_map_t> choi_apstate_to_exact_circuit(
+    ChoiAPState ap, CXConfigType cx_config = CXConfigType::Snake);
+
+/**
+ * Constructs a unitary circuit whose stabilizer group contains all the
+ * stabilizers of the ChoiAPState and possibly more. See
+ * cm_tableau_to_unitary_extension_circuit for more details.
+ */
+std::pair<Circuit, qubit_map_t> choi_apstate_to_unitary_extension_circuit(
+    ChoiAPState ap, const std::vector<Qubit> &init_names = {},
+    const std::vector<Qubit> &post_names = {},
+    CXConfigType cx_config = CXConfigType::Snake);
+
+/**
+ * Convert a SymplecticTableau describing a stabiliser state to AP-form. Since
+ * tableau methods do not track global phase, we set it to 0.
+ */
+APState tableau_to_apstate(SymplecticTableau tab);
+
+/**
+ * Convert an APState describing a stabiliser state into a tableau form,
+ * discarding the global phase information. This allows us to reuse tableau
+ * synthesis methods to synthesise APState objects (the result would need
+ * re-simulating as an APState to compare and deduce the required global phase).
+ */
+SymplecticTableau apstate_to_tableau(APState ap);
+
+/**
+ * Convert a ChoiMixTableau describing a stabiliser state to AP-form. Since
+ * tableau methods do not track global phase, we set it to 0.
+ */
+ChoiAPState cm_tableau_to_choi_apstate(const ChoiMixTableau &tab);
+
+/**
+ * Convert a ChoiAPState describing a stabiliser state into a tableau form,
+ * discarding the global phase information. This allows us to reuse tableau
+ * synthesis methods to synthesise ChoiAPState objects (the result would need
+ * re-simulating as a ChoiAPState to compare and deduce the required global
+ * phase).
+ */
+ChoiMixTableau choi_apstate_to_cm_tableau(const ChoiAPState &ap);
+
 PauliGraph circuit_to_pauli_graph(const Circuit &circ);
 
 /**
diff --git a/tket/src/Clifford/APState.cpp b/tket/src/Clifford/APState.cpp
new file mode 100644
index 0000000000..c9d96d241b
--- /dev/null
+++ b/tket/src/Clifford/APState.cpp
@@ -0,0 +1,1473 @@
+// Copyright 2019-2024 Cambridge Quantum Computing
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "tket/Clifford/APState.hpp"
+
+#include <boost/foreach.hpp>
+#include <tkassert/Assert.hpp>
+
+#include "tket/OpType/OpTypeInfo.hpp"
+
+namespace tket {
+
+/*************************** APState IMPLEMENTATION ***************************/
+
+APState::APState(
+    const MatrixXb& A_, const VectorXb& B_, const MatrixXb& C_,
+    const MatrixXb& E_, const Eigen::VectorXi& P_, const Expr& phase_)
+    : A(A_), B(B_), C(C_), E(E_), P(P_), phase(phase_) {
+  verify();
+}
+
+APState::APState(unsigned n_qubits)
+    : A(MatrixXb::Identity(n_qubits, n_qubits)),
+      B(VectorXb::Zero(n_qubits)),
+      C(MatrixXb::Zero(n_qubits, n_qubits)),
+      E(MatrixXb::Zero(n_qubits, n_qubits)),
+      P(Eigen::VectorXi::Zero(n_qubits)),
+      phase(0) {}
+
+unsigned clifford_phase(const Complex& c) {
+  unsigned res = 0;
+  unsigned n_possibles = 0;
+  if (c.real() > EPS) {
+    res = 0;
+    ++n_possibles;
+  }
+  if (c.imag() > EPS) {
+    res = 1;
+    ++n_possibles;
+  }
+  if (c.real() < -EPS) {
+    res = 2;
+    ++n_possibles;
+  }
+  if (c.imag() < -EPS) {
+    res = 3;
+    ++n_possibles;
+  }
+  TKET_ASSERT(n_possibles == 1);
+  return res;
+}
+
+APState::APState(const Eigen::VectorXcd& sv) {
+  unsigned n_qbs = 0;
+  while (sv.size() > (1 << n_qbs)) ++n_qbs;
+  TKET_ASSERT(sv.size() == (1 << n_qbs));
+
+  // Find non-zero entries as a vector space and offset
+  unsigned z0 = 0;
+  for (unsigned x = 0; x < sv.size(); ++x) {
+    if (sv(x) != 0.) {
+      z0 = x;
+      break;
+    }
+  }
+  std::vector<unsigned> offsets;
+  unsigned n_non_zero = 0;
+  for (unsigned x = 1; x < sv.size(); ++x) {
+    if (sv(z0 ^ x) != 0.) {
+      ++n_non_zero;
+      if (n_non_zero == (1u << offsets.size())) offsets.push_back(x);
+    }
+  }
+
+  // Find A as the dual space
+  MatrixXb offset_mat(offsets.size(), n_qbs);
+  for (unsigned r = 0; r < offsets.size(); ++r) {
+    unsigned off = offsets.at(r);
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      // Binary encoding of offsets in reverse order to guarantee free qubits
+      // are the later ones, meaning we produce A in row echelon form
+      offset_mat(r, c) = (off & (1 << c)) != 0;
+    }
+  }
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(offset_mat);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      offset_mat(op.second, j) ^= offset_mat(op.first, j);
+    }
+  }
+  std::map<unsigned, unsigned> mat_leaders;
+  A = MatrixXb::Zero(n_qbs, n_qbs);
+  for (unsigned c = 0; c < n_qbs; ++c) {
+    bool free_qubit = false;
+    for (unsigned r = 0; r < offsets.size(); ++r) {
+      if (offset_mat(r, c)) {
+        auto inserted = mat_leaders.insert({r, c});
+        if (inserted.second) {
+          free_qubit = true;
+          break;
+        } else {
+          // Reverse bit orderings back to normal
+          A(n_qbs - 1 - c, n_qbs - 1 - inserted.first->second) = true;
+        }
+      }
+    }
+    A(n_qbs - 1 - c, n_qbs - 1 - c) = !free_qubit;
+  }
+
+  VectorXb z0_vec(n_qbs);
+  for (unsigned i = 0; i < n_qbs; ++i) {
+    z0_vec(i) = ((z0 & (1 << (n_qbs - 1 - i))) != 0);
+  }
+  B = A * z0_vec;
+
+  E = MatrixXb::Zero(n_qbs, n_qbs);
+  P = Eigen::VectorXi::Zero(n_qbs);
+  unsigned neutral_z0 = z0;  // Index with 0 for all free qubits
+  std::map<unsigned, unsigned> offset_for_free;
+  for (const std::pair<const unsigned, unsigned>& row_qfree : mat_leaders) {
+    unsigned offset = 0;
+    for (unsigned i = 0; i < n_qbs; ++i)
+      if (offset_mat(row_qfree.first, i)) offset += (1 << i);
+    unsigned qfree = n_qbs - 1 - row_qfree.second;
+    offset_for_free.insert({qfree, offset});
+    if ((neutral_z0 & (1 << row_qfree.second)) != 0)
+      neutral_z0 = neutral_z0 ^ offset;
+  }
+  for (const std::pair<const unsigned, unsigned>& qfree_offset :
+       offset_for_free) {
+    unsigned qfree = qfree_offset.first;
+    unsigned offset = qfree_offset.second;
+    auto complex_phase = sv(neutral_z0 ^ offset) / sv(neutral_z0);
+    P(qfree) = clifford_phase(complex_phase);
+    for (const std::pair<const unsigned, unsigned>& qfree_offset2 :
+         offset_for_free) {
+      if (qfree_offset == qfree_offset2) break;  // Only go up to solved phases
+      unsigned qfree2 = qfree_offset2.first;
+      unsigned offset2 = qfree_offset2.second;
+      auto complex_phase_cz =
+          sv(neutral_z0 ^ offset ^ offset2) / sv(neutral_z0);
+      E(qfree, qfree2) = E(qfree2, qfree) =
+          ((unsigned)(clifford_phase(complex_phase_cz) - P(qfree) - P(qfree2)) %
+           4) == 2;
+    }
+  }
+
+  phase = Expr(std::arg(sv(neutral_z0)) / PI);
+}
+
+APState::APState(const Eigen::MatrixXcd& dm) {
+  Eigen::VectorXcd dm_as_vec = dm.reshaped();
+  APState pure_double(dm_as_vec);
+  unsigned n_qbs = pure_double.A.cols() / 2;
+  A = MatrixXb::Zero(n_qbs, n_qbs);
+  B = VectorXb::Zero(n_qbs);
+  C = MatrixXb::Zero(n_qbs, n_qbs);
+
+  MatrixXb full_mat(2 * n_qbs, 2 * n_qbs + 1);
+  full_mat.block(0, 0, 2 * n_qbs, 2 * n_qbs) = pure_double.A;
+  full_mat.col(2 * n_qbs) = pure_double.B;
+
+  // This full matrix is generated by AI, IA, and CC. Gaussian elimination puts
+  // zeros at the bottom, followed by IA.
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(full_mat);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < 2 * n_qbs + 1; ++j) {
+      full_mat(op.second, j) ^= full_mat(op.first, j);
+    }
+  }
+
+  unsigned first_zero = 0;
+  for (unsigned r = n_qbs; r > 0;) {
+    --r;
+    bool is_zero = true;
+    // By symmetry of AI, IA, and CC, the bottom non-empty row must have some
+    // component on the right side
+    for (unsigned c = n_qbs; c < 2 * n_qbs; ++c) {
+      if (full_mat(r, c)) {
+        is_zero = false;
+        break;
+      }
+    }
+    if (!is_zero) {
+      first_zero = r + 1;
+      break;
+    }
+  }
+  unsigned first_right = 0;
+  for (unsigned r = first_zero; r > 0;) {
+    --r;
+    bool is_right_only = true;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (full_mat(r, c)) {
+        is_right_only = false;
+        break;
+      }
+    }
+    if (!is_right_only) {
+      first_right = r + 1;
+      break;
+    }
+  }
+  A.block(0, 0, first_zero - first_right, n_qbs) =
+      full_mat.block(first_right, n_qbs, first_zero - first_right, n_qbs);
+  B.head(first_zero - first_right) =
+      full_mat.block(first_right, 2 * n_qbs, first_zero - first_right, 1);
+
+  // Flip the column order and reduce the remaining rows to obtain CC above AI;
+  // get the row combinations from the reordered matrix but apply them to the
+  // matrix with the correct column ordering
+  MatrixXb remaining_rows(first_right, 2 * n_qbs);
+  for (unsigned r = 0; r < first_right; ++r) {
+    for (unsigned c = 0; c < 2 * n_qbs; ++c) {
+      remaining_rows(r, c) = full_mat(r, 2 * n_qbs - 1 - c);
+    }
+  }
+  row_ops = gaussian_elimination_row_ops(remaining_rows);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < 2 * n_qbs; ++j) {
+      full_mat(op.second, j) ^= full_mat(op.first, j);
+    }
+  }
+
+  unsigned first_left = 0;
+  for (unsigned r = first_right; r > 0;) {
+    --r;
+    bool is_left_only = true;
+    for (unsigned c = n_qbs; c < 2 * n_qbs; ++c) {
+      if (full_mat(r, c)) {
+        is_left_only = false;
+        break;
+      }
+    }
+    if (!is_left_only) {
+      first_left = r + 1;
+      break;
+    }
+  }
+  C.block(0, 0, first_left, n_qbs) = full_mat.block(0, 0, first_left, n_qbs);
+
+  /**
+   * Recall that pure_double.A is generated by AI, IA, and CC. The statevector
+   * constructor builds this in normal form, so it is in reduced row-echelon
+   * form. The first block of qubits may have leaders for leaders or mixed
+   * qubits in the density matrix, but the second block only has them for true
+   * leaders in the density matrix. This means pure_double.E and pure_double.P
+   * are not just the direct sum of two copies of the true E and P, but the
+   * bottom segments contain them exactly.
+   */
+  E = pure_double.E.block(n_qbs, n_qbs, n_qbs, n_qbs);
+  P = pure_double.P.tail(n_qbs);
+
+  phase = 0;
+}
+
+void APState::verify() const {
+  // Check dimensions all agree
+  unsigned n_qubits = A.rows();
+  TKET_ASSERT(A.cols() == n_qubits);
+  TKET_ASSERT(B.size() == n_qubits);
+  TKET_ASSERT(C.rows() == n_qubits);
+  TKET_ASSERT(C.cols() == n_qubits);
+  TKET_ASSERT(E.rows() == n_qubits);
+  TKET_ASSERT(E.cols() == n_qubits);
+  TKET_ASSERT(P.size() == n_qubits);
+
+  for (unsigned r = 0; r < n_qubits; ++r) {
+    for (unsigned c = 0; c < r; ++c) {
+      // Check E is symmetric
+      TKET_ASSERT(E(r, c) == E(c, r));
+    }
+    // Check diagonals of E are zero
+    TKET_ASSERT(!E(r, r));
+  }
+}
+
+VectorXb z2_mult(const MatrixXb& M, const VectorXb& v) {
+  VectorXb res = VectorXb::Zero(M.rows());
+  for (unsigned i = 0; i < M.cols(); ++i) {
+    if (v(i)) {
+      for (unsigned j = 0; j < M.rows(); ++j) res(j) ^= M(j, i);
+    }
+  }
+  return res;
+}
+
+Eigen::VectorXcd APState::to_statevector() const {
+  unsigned n_qubits = A.cols();
+  Eigen::VectorXcd sv = Eigen::VectorXcd::Zero(1 << n_qubits);
+  unsigned n_terms = 0;
+  Complex g_phase = std::exp(i_ * PI * eval_expr(phase).value());
+  for (unsigned x = 0; x < (1u << n_qubits); ++x) {
+    VectorXb x_binary = VectorXb::Zero(n_qubits);
+    for (unsigned i = 0; i < n_qubits; ++i) {
+      unsigned mask = 1 << (n_qubits - 1 - i);
+      x_binary(i) = ((x & mask) != 0);
+    }
+
+    if (z2_mult(A, x_binary) == B) {
+      ++n_terms;
+      unsigned i_phases = 0;
+      for (unsigned q = 0; q < n_qubits; ++q) {
+        if (x_binary(q)) {
+          i_phases += P(q);
+          for (unsigned q2 = q; q2 < n_qubits; ++q2) {
+            if (E(q, q2) && x_binary(q2)) i_phases += 2;
+          }
+        }
+      }
+      switch (i_phases % 4) {
+        case 0: {
+          sv(x) = g_phase;
+          break;
+        }
+        case 1: {
+          sv(x) = i_ * g_phase;
+          break;
+        }
+        case 2: {
+          sv(x) = -g_phase;
+          break;
+        }
+        case 3: {
+          sv(x) = -i_ * g_phase;
+          break;
+        }
+        default: {
+          TKET_ASSERT(false);
+        }
+      }
+    }
+  }
+  return pow(n_terms, -0.5) * sv;
+}
+
+Eigen::MatrixXcd APState::to_density_matrix() const {
+  unsigned n_qubits = A.cols();
+  Eigen::MatrixXcd dm = Eigen::MatrixXcd::Zero(1 << n_qubits, 1 << n_qubits);
+  for (unsigned x = 0; x < (1u << n_qubits); ++x) {
+    VectorXb x_binary = VectorXb::Zero(n_qubits);
+    for (unsigned i = 0; i < n_qubits; ++i) {
+      unsigned mask = 1 << (n_qubits - 1 - i);
+      x_binary(i) = ((x & mask) != 0);
+    }
+
+    if (z2_mult(A, x_binary) == B) {
+      MatrixXb Cx = z2_mult(C, x_binary);
+      for (unsigned y = 0; y < (1u << n_qubits); ++y) {
+        VectorXb y_binary = VectorXb::Zero(n_qubits);
+        for (unsigned i = 0; i < n_qubits; ++i) {
+          unsigned mask = 1 << (n_qubits - 1 - i);
+          y_binary(i) = ((y & mask) != 0);
+        }
+
+        if ((z2_mult(A, y_binary) == B) && (z2_mult(C, y_binary) == Cx)) {
+          unsigned i_phases = 0;
+          for (unsigned q = 0; q < n_qubits; ++q) {
+            if (x_binary(q)) {
+              i_phases += P(q);
+              for (unsigned q2 = q; q2 < n_qubits; ++q2) {
+                if (E(q, q2) && x_binary(q2)) i_phases += 2;
+              }
+            }
+            if (y_binary(q)) {
+              i_phases -= P(q);
+              for (unsigned q2 = q; q2 < n_qubits; ++q2) {
+                if (E(q, q2) && y_binary(q2)) i_phases += 2;
+              }
+            }
+          }
+          switch (i_phases % 4) {
+            case 0: {
+              dm(x, y) = 1;
+              break;
+            }
+            case 1: {
+              dm(x, y) = i_;
+              break;
+            }
+            case 2: {
+              dm(x, y) = -1;
+              break;
+            }
+            case 3: {
+              dm(x, y) = -i_;
+              break;
+            }
+            default: {
+              TKET_ASSERT(false);
+            }
+          }
+        }
+      }
+    }
+  }
+  return dm / dm.trace();
+}
+
+void APState::apply_CZ(unsigned ctrl, unsigned trgt) {
+  E(ctrl, trgt) ^= true;
+  E(trgt, ctrl) ^= true;
+}
+
+void APState::apply_S(unsigned q) { P(q) += 1; }
+
+void APState::apply_V(unsigned q) {
+  unsigned n_qbs = A.cols();
+  std::list<unsigned> A_rows;
+  std::list<unsigned> C_rows;
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    if (A(r, q)) A_rows.push_back(r);
+    if (C(r, q)) C_rows.push_back(r);
+  }
+  if ((unsigned)P(q) % 2 == 0) {
+    bool a = ((unsigned)P(q) % 4) == 2;
+    if (A_rows.empty()) {
+      for (unsigned q2 = 0; q2 < n_qbs; ++q2) {
+        if (E(q, q2)) {
+          // Update local phase on neighbours
+          P(q2) += a ? 3 : 1;
+          // Local complementation between neighbours
+          for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) {
+            if (E(q, q3)) {
+              E(q2, q3) ^= true;
+              E(q3, q2) ^= true;
+            }
+          }
+          // Add connections to all C_rows
+          for (unsigned r : C_rows) {
+            C(r, q2) ^= true;
+          }
+        }
+      }
+      // Global phase
+      phase += a ? .25 : -.25;
+    } else {
+      // Choose one of the red spiders to remove
+      unsigned r = A_rows.back();
+      A_rows.pop_back();
+      // Stratify neighbourhoods of r (in A) and q (in E)
+      std::list<unsigned> just_r;
+      std::list<unsigned> just_q;
+      std::list<unsigned> both;
+      for (unsigned q2 = 0; q2 < n_qbs; ++q2) {
+        if (q2 == q) continue;
+        if (A(r, q2)) {
+          if (E(q, q2))
+            both.push_back(q2);
+          else
+            just_r.push_back(q2);
+        } else if (E(q, q2))
+          just_q.push_back(q2);
+      }
+      // Update A and B
+      for (unsigned ar : A_rows) {
+        A(ar, q) = false;
+        for (unsigned q2 : just_r) A(ar, q2) ^= true;
+        for (unsigned q2 : both) A(ar, q2) ^= true;
+        B(ar) ^= B(r);
+      }
+      // Update C
+      for (unsigned cr : C_rows) {
+        C(cr, q) = false;
+        for (unsigned q2 : just_r) C(cr, q2) ^= true;
+        for (unsigned q2 : both) C(cr, q2) ^= true;
+      }
+      // Update E and P
+      for (unsigned rn : just_r) {
+        // Complementation within just_r
+        for (unsigned rn2 : just_r) {
+          E(rn, rn2) ^= true;
+        }
+        // Reset diagonal
+        E(rn, rn) = false;
+        // Complementation between just_r and just_q
+        for (unsigned qn : just_q) {
+          E(rn, qn) ^= true;
+          E(qn, rn) ^= true;
+        }
+        // Connect to q
+        E(rn, q) = true;
+        E(q, rn) = true;
+        // Local phases
+        P(rn) += (a ^ B(r)) ? 1 : 3;
+      }
+      for (unsigned bn : both) {
+        // Complementation within both
+        for (unsigned bn2 : both) {
+          E(bn, bn2) ^= true;
+        }
+        // Reset diagonal
+        E(bn, bn) = false;
+        // Complementation between both and just_q
+        for (unsigned qn : just_q) {
+          E(bn, qn) ^= true;
+          E(qn, bn) ^= true;
+        }
+        // Local phases
+        P(bn) += a ? 3 : 1;
+      }
+      for (unsigned qn : just_q) {
+        // Remove connection to q
+        E(q, qn) = false;
+        E(qn, q) = false;
+        // Local phases
+        P(qn) += B(r) ? 2 : 0;
+      }
+      // Local phase on q
+      P(q) = B(r) ? 1 : 3;
+      // Global phase
+      if (B(r)) phase += a ? .5 : 1.5;
+      // Remove row r from A and B
+      for (unsigned i = 0; i < n_qbs; ++i) {
+        A(r, i) = false;
+      }
+      B(r) = false;
+    }
+  } else {
+    bool a = ((unsigned)P(q) % 4) == 3;
+    if (!A_rows.empty()) {
+      // Choose one of the red spiders to remove
+      unsigned r = A_rows.back();
+      A_rows.pop_back();
+      // Stratify neighbourhoods of r (in A) and q (in E)
+      std::list<unsigned> just_r;
+      std::list<unsigned> just_q;
+      std::list<unsigned> both;
+      for (unsigned q2 = 0; q2 < n_qbs; ++q2) {
+        if (q2 == q) continue;
+        if (A(r, q2)) {
+          if (E(q, q2))
+            both.push_back(q2);
+          else
+            just_r.push_back(q2);
+        } else if (E(q, q2))
+          just_q.push_back(q2);
+      }
+      // Update A and B
+      for (unsigned ar : A_rows) {
+        A(ar, q) = false;
+        for (unsigned q2 : just_r) A(ar, q2) ^= true;
+        for (unsigned q2 : both) A(ar, q2) ^= true;
+        B(ar) ^= B(r);
+      }
+      // Update C
+      for (unsigned cr : C_rows) {
+        C(cr, q) = false;
+        for (unsigned q2 : just_r) C(cr, q2) ^= true;
+        for (unsigned q2 : both) C(cr, q2) ^= true;
+      }
+      // Update E and P
+      for (unsigned rn : just_r) {
+        // Complementation between just_r and just_q
+        for (unsigned qn : just_q) {
+          E(rn, qn) ^= true;
+          E(qn, rn) ^= true;
+        }
+        // Complementation between just_r and both
+        for (unsigned bn : both) {
+          E(rn, bn) ^= true;
+          E(bn, rn) ^= true;
+        }
+        // Connect to q
+        E(rn, q) = true;
+        E(q, rn) = true;
+        // Local phases
+        P(rn) += a ? 2 : 0;
+      }
+      for (unsigned bn : both) {
+        // Complementation between both and just_q
+        for (unsigned qn : just_q) {
+          E(bn, qn) ^= true;
+          E(qn, bn) ^= true;
+        }
+        // Local phases
+        P(bn) += (a ^ B(r)) ? 0 : 2;
+      }
+      for (unsigned qn : just_q) {
+        // Remove connection to q
+        E(qn, q) = false;
+        E(q, qn) = false;
+        // Local phases
+        P(qn) += B(r) ? 2 : 0;
+      }
+      // Local phase on q
+      P(q) = B(r) ? 1 : 3;
+      // Global phase
+      if (a && B(r)) phase += 1;
+      // Remove row r from A and B
+      for (unsigned i = 0; i < n_qbs; ++i) {
+        A(r, i) = false;
+      }
+      B(r) = false;
+    } else if (!C_rows.empty()) {
+      // Choose one of the discards to pivot around
+      unsigned d = C_rows.back();
+      C_rows.pop_back();
+      // Stratify neighbourhoods of d (in C) and q (in E)
+      std::list<unsigned> just_d;
+      std::list<unsigned> just_q;
+      std::list<unsigned> both;
+      for (unsigned q2 = 0; q2 < n_qbs; ++q2) {
+        if (q2 == q) continue;
+        if (C(d, q2)) {
+          if (E(q, q2))
+            both.push_back(q2);
+          else
+            just_d.push_back(q2);
+        } else if (E(q, q2))
+          just_q.push_back(q2);
+      }
+      // Update C
+      for (unsigned cr : C_rows) {
+        C(cr, q) = false;
+        for (unsigned q2 : just_d) C(cr, q2) ^= true;
+        for (unsigned q2 : both) C(cr, q2) ^= true;
+      }
+      // Update E and P
+      for (unsigned dn : just_d) {
+        // Complementation between just_d and just_q
+        for (unsigned qn : just_q) {
+          E(dn, qn) ^= true;
+          E(qn, dn) ^= true;
+        }
+        // Complementation between just_d and both
+        for (unsigned bn : both) {
+          E(dn, bn) ^= true;
+          E(bn, dn) ^= true;
+        }
+        // Connect to q
+        E(dn, q) = true;
+        E(q, dn) = true;
+        // Remove connection with d
+        C(d, dn) = false;
+        // Local phases
+        P(dn) += a ? 2 : 0;
+      }
+      for (unsigned bn : both) {
+        // Complementation between both and just_q
+        for (unsigned qn : just_q) {
+          E(bn, qn) ^= true;
+          E(qn, bn) ^= true;
+        }
+        // Local phases
+        P(bn) += a ? 0 : 2;
+      }
+      for (unsigned qn : just_q) {
+        // Connect to d
+        C(d, qn) = true;
+        // Remove connection to q
+        E(qn, q) = false;
+        E(q, qn) = false;
+        // No local phase change
+      }
+      // Local phase on q
+      P(q) = 3;
+      // No global phase change
+    } else {
+      VectorXb new_row = VectorXb::Zero(n_qbs);
+      new_row(q) = true;
+      for (unsigned q2 = 0; q2 < n_qbs; ++q2) {
+        if (E(q, q2)) {
+          // Connect to new red spider
+          new_row(q2) = true;
+          // Local phase on neighbours
+          P(q2) += a ? 1 : 3;
+          // Local complementation between neighbours
+          for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) {
+            if (E(q, q3)) {
+              E(q2, q3) ^= true;
+              E(q3, q2) ^= true;
+            }
+          }
+          // Remove connection with q
+          E(q, q2) = false;
+          E(q2, q) = false;
+        }
+      }
+      // Reset local phase on q
+      P(q) = 0;
+      // Global phase
+      phase += a ? 1.5 : 0;
+      // Add new_row to A and B
+      MatrixXb combined_mat = MatrixXb::Zero(n_qbs + 1, n_qbs + 1);
+      combined_mat.block(0, 0, n_qbs, n_qbs) = A;
+      combined_mat.block(0, n_qbs, n_qbs, 1) = B;
+      combined_mat.block(n_qbs, 0, 1, n_qbs) = new_row.transpose();
+      combined_mat(n_qbs, n_qbs) = a;
+      std::vector<std::pair<unsigned, unsigned>> row_ops =
+          gaussian_elimination_row_ops(combined_mat);
+      for (const std::pair<unsigned, unsigned>& op : row_ops) {
+        for (unsigned j = 0; j < n_qbs + 1; ++j) {
+          combined_mat(op.second, j) ^= combined_mat(op.first, j);
+        }
+      }
+      A = combined_mat.block(0, 0, n_qbs, n_qbs);
+      B = combined_mat.block(0, n_qbs, n_qbs, 1);
+    }
+  }
+}
+
+void APState::apply_X(unsigned q) {
+  // Push through the local phase
+  phase += P(q) * .5;
+  P(q) = -P(q);
+  // Pushing through the CZs adds Zs on neighbours
+  for (unsigned q2 = 0; q2 < E.cols(); ++q2) {
+    if (E(q, q2)) P(q2) += 2;
+  }
+  // Pushing through adjacency matrix adds Xs onto connected reds
+  for (unsigned r = 0; r < A.rows(); ++r) {
+    if (A(r, q)) B(r) ^= true;
+  }
+}
+
+void APState::apply_gate(OpType type, const std::vector<unsigned>& qbs) {
+  switch (type) {
+    case OpType::Z: {
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(0));
+      break;
+    }
+    case OpType::X: {
+      apply_X(qbs.at(0));
+      break;
+    }
+    case OpType::Y: {
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(0));
+      apply_X(qbs.at(0));
+      phase += .5;
+      break;
+    }
+    case OpType::S: {
+      apply_S(qbs.at(0));
+      break;
+    }
+    case OpType::Sdg: {
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(0));
+      break;
+    }
+    case OpType::V: {
+      apply_V(qbs.at(0));
+      break;
+    }
+    case OpType::SX: {
+      apply_V(qbs.at(0));
+      phase += .25;
+      break;
+    }
+    case OpType::Vdg: {
+      apply_V(qbs.at(0));
+      apply_X(qbs.at(0));
+      phase += .5;
+      break;
+    }
+    case OpType::SXdg: {
+      apply_V(qbs.at(0));
+      apply_X(qbs.at(0));
+      phase += .25;
+      break;
+    }
+    case OpType::H: {
+      apply_S(qbs.at(0));
+      apply_V(qbs.at(0));
+      apply_S(qbs.at(0));
+      break;
+    }
+    case OpType::CX: {
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_S(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_S(qbs.at(1));
+      break;
+    }
+    case OpType::CY: {
+      apply_V(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_X(qbs.at(1));
+      phase += .5;
+      break;
+    }
+    case OpType::CZ: {
+      apply_CZ(qbs.at(0), qbs.at(1));
+      break;
+    }
+    case OpType::ZZMax: {
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      phase -= .25;
+      break;
+    }
+    case OpType::ECR: {
+      apply_S(qbs.at(0));
+      apply_X(qbs.at(0));
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_S(qbs.at(1));
+      phase += -.25;
+      break;
+    }
+    case OpType::ISWAPMax: {
+      apply_V(qbs.at(0));
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_V(qbs.at(0));
+      apply_V(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_V(qbs.at(0));
+      apply_V(qbs.at(1));
+      apply_S(qbs.at(1));
+      phase += 1.;
+      break;
+    }
+    case OpType::SWAP: {
+      apply_S(qbs.at(0));
+      apply_V(qbs.at(0));
+      apply_S(qbs.at(0));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_S(qbs.at(0));
+      apply_V(qbs.at(0));
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_S(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_S(qbs.at(0));
+      apply_V(qbs.at(0));
+      apply_S(qbs.at(0));
+      apply_S(qbs.at(1));
+      apply_V(qbs.at(1));
+      apply_S(qbs.at(1));
+      apply_CZ(qbs.at(0), qbs.at(1));
+      apply_S(qbs.at(0));
+      apply_V(qbs.at(0));
+      apply_S(qbs.at(0));
+      break;
+    }
+    case OpType::BRIDGE: {
+      apply_S(qbs.at(2));
+      apply_V(qbs.at(2));
+      apply_S(qbs.at(2));
+      apply_CZ(qbs.at(0), qbs.at(2));
+      apply_S(qbs.at(2));
+      apply_V(qbs.at(2));
+      apply_S(qbs.at(2));
+      break;
+    }
+    case OpType::noop: {
+      break;
+    }
+    case OpType::Reset: {
+      unsigned q = qbs.at(0);
+      discard_qubit(q);
+      unsigned new_qb = init_qubit();
+      // Swap qubits to preserve the ordering
+      if (q != new_qb) {
+        unsigned n_qbs = A.cols();
+        for (unsigned r = 0; r < n_qbs; ++r) {
+          bool temp = A(r, q);
+          A(r, q) = A(r, new_qb);
+          A(r, new_qb) = temp;
+        }
+        for (unsigned r = 0; r < n_qbs; ++r) {
+          bool temp = C(r, q);
+          C(r, q) = C(r, new_qb);
+          C(r, new_qb) = temp;
+        }
+        for (unsigned r = 0; r < n_qbs; ++r) {
+          bool temp = E(r, q);
+          E(r, q) = E(r, new_qb);
+          E(r, new_qb) = temp;
+        }
+        for (unsigned c = 0; c < n_qbs; ++c) {
+          bool temp = E(q, c);
+          E(q, c) = E(new_qb, c);
+          E(new_qb, c) = temp;
+        }
+        int temp = P(q);
+        P(q) = P(new_qb);
+        P(new_qb) = temp;
+      }
+      break;
+    }
+    case OpType::Collapse: {
+      collapse_qubit(qbs.at(0));
+      break;
+    }
+    case OpType::Phase: {
+      throw std::logic_error(
+          "OpType::Phase cannot be applied via APState::apply_gate");
+    }
+    default: {
+      throw BadOpType(
+          "Cannot be applied to a APState: not a Clifford gate", type);
+    }
+  }
+}
+
+unsigned APState::init_qubit() {
+  unsigned n_qbs = A.cols();
+  A.conservativeResize(n_qbs + 1, n_qbs + 1);
+  A.col(n_qbs) = VectorXb::Zero(n_qbs + 1);
+  A.row(n_qbs) = VectorXb::Zero(n_qbs + 1);
+  A(n_qbs, n_qbs) = true;
+  B.conservativeResize(n_qbs + 1);
+  B(n_qbs) = false;
+  C.conservativeResize(n_qbs + 1, n_qbs + 1);
+  C.col(n_qbs) = VectorXb::Zero(n_qbs + 1);
+  C.row(n_qbs) = VectorXb::Zero(n_qbs + 1);
+  E.conservativeResize(n_qbs + 1, n_qbs + 1);
+  E.col(n_qbs) = VectorXb::Zero(n_qbs + 1);
+  E.row(n_qbs) = VectorXb::Zero(n_qbs + 1);
+  P.conservativeResize(n_qbs + 1);
+  P(n_qbs) = 0;
+  return n_qbs;
+}
+
+unsigned APState::post_select(unsigned q) {
+  unsigned n_qbs = A.cols();
+  if (q >= n_qbs)
+    throw std::invalid_argument("APState: post-selecting a non-existent qubit");
+  MatrixXb AB = MatrixXb::Zero(n_qbs, n_qbs);
+  AB.block(0, 0, n_qbs, n_qbs) = A;
+  AB.col(q) = A.col(n_qbs - 1);
+  AB.col(n_qbs - 1) = B;
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(AB);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      AB(op.second, j) ^= AB(op.first, j);
+    }
+  }
+  A.resize(n_qbs - 1, n_qbs - 1);
+  A = AB.block(0, 0, n_qbs - 1, n_qbs - 1);
+  B.resize(n_qbs - 1);
+  B = AB.block(0, n_qbs - 1, n_qbs - 1, 1);
+  C.col(q) = C.col(n_qbs - 1);
+  row_ops = gaussian_elimination_row_ops(C);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      C(op.second, j) ^= C(op.first, j);
+    }
+  }
+  C.conservativeResize(n_qbs - 1, n_qbs - 1);
+  E.col(q) = E.col(n_qbs - 1);
+  E.row(q) = E.row(n_qbs - 1);
+  E.conservativeResize(n_qbs - 1, n_qbs - 1);
+  P(q) = P(n_qbs - 1);
+  P.conservativeResize(n_qbs - 1);
+  return n_qbs - 1;
+}
+
+void APState::collapse_qubit(unsigned q) {
+  unsigned n_qbs = A.cols();
+  MatrixXb C_ext = MatrixXb::Zero(n_qbs + 1, n_qbs);
+  C_ext.block(0, 0, n_qbs, n_qbs) = C;
+  C_ext(n_qbs, q) = 1;
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(C_ext);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      C_ext(op.second, j) ^= C_ext(op.first, j);
+    }
+  }
+  C = C_ext.block(0, 0, n_qbs, n_qbs);
+}
+
+unsigned APState::discard_qubit(unsigned q) {
+  collapse_qubit(q);
+  apply_V(q);
+  collapse_qubit(q);
+  return post_select(q);
+}
+
+void APState::normal_form() {
+  unsigned n_qbs = A.cols();
+  // Get A into reduced row-echelon form
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(A);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      A(op.second, j) ^= A(op.first, j);
+    }
+    B(op.second) ^= B(op.first);
+  }
+  // Identify leading qubits
+  std::map<unsigned, unsigned> leader_to_row;
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    bool leader_found = false;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (A(r, c)) {
+        leader_found = true;
+        leader_to_row.insert({c, r});
+        break;
+      }
+    }
+    if (!leader_found) break;
+  }
+  // Remove leaders from C
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    for (const std::pair<const unsigned, unsigned>& lr : leader_to_row) {
+      if (C(r, lr.first)) {
+        for (unsigned c = 0; c < n_qbs; ++c) {
+          C(r, c) ^= A(lr.second, c);
+        }
+      }
+    }
+  }
+  // Remove leaders from E
+  for (const std::pair<const unsigned, unsigned>& lr : leader_to_row) {
+    for (unsigned q2 = lr.first; q2 < n_qbs; ++q2) {
+      if (E(lr.first, q2)) {
+        E(lr.first, q2) = false;
+        E(q2, lr.first) = false;
+        for (unsigned q3 = lr.first + 1; q3 < n_qbs; ++q3) {
+          if (A(lr.second, q3)) {
+            E(q2, q3) ^= true;
+            E(q3, q2) ^= true;
+          }
+        }
+        P(q2) += (B(lr.second) ^ A(lr.second, q2)) ? 2 : 0;
+      }
+    }
+  }
+  // Remove leaders from P
+  for (const std::pair<const unsigned, unsigned>& lr : leader_to_row) {
+    if ((unsigned)P(lr.first) % 2 == 1) {
+      // Complementation within neighbours under A (excluding the leader)
+      for (unsigned q2 = lr.first + 1; q2 < n_qbs; ++q2) {
+        if (A(lr.second, q2)) {
+          for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) {
+            if (A(lr.second, q3)) {
+              E(q2, q3) ^= true;
+              E(q3, q2) ^= true;
+            }
+          }
+          // Local phases of neighbours
+          P(q2) += (P(lr.first) % 4) + (B(lr.second) ? 2 : 0);
+        }
+      }
+    } else if ((unsigned)P(lr.first) % 4 == 2) {
+      // Local phases of neighbours
+      for (unsigned q2 = lr.first + 1; q2 < n_qbs; ++q2) {
+        if (A(lr.second, q2)) P(q2) += 2;
+      }
+    }
+    // Global phase
+    if (B(lr.second)) phase += P(lr.first) * .5;
+    // Reset P
+    P(lr.first) = 0;
+  }
+  // Get C into reduced row-echelon form
+  row_ops = gaussian_elimination_row_ops(C);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      C(op.second, j) ^= C(op.first, j);
+    }
+  }
+  // Identify mixed qubits
+  std::map<unsigned, unsigned> mixed_to_row;
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    bool mixed_found = false;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (C(r, c)) {
+        mixed_found = true;
+        mixed_to_row.insert({c, r});
+        break;
+      }
+    }
+    if (!mixed_found) break;
+  }
+  // Remove E connections between mixed qubits
+  for (auto mr1 = mixed_to_row.begin(); mr1 != mixed_to_row.end(); ++mr1) {
+    for (auto mr2 = mixed_to_row.begin(); mr2 != mr1; ++mr2) {
+      if (E(mr1->first, mr2->first)) {
+        // Local complementation along the sum of their rows in C
+        // mr2->first < mr1->first, and they are the first entries in their rows
+        for (unsigned i = mr2->first; i < n_qbs; ++i) {
+          if (C(mr1->second, i) ^ C(mr2->second, i)) {
+            for (unsigned j = i + 1; j < n_qbs; ++j) {
+              if (C(mr1->second, j) ^ C(mr2->second, j)) {
+                E(i, j) ^= true;
+                E(j, i) ^= true;
+              }
+            }
+            // Add local phases around neighbourhood
+            P(i) += 1;
+          }
+        }
+      }
+    }
+  }
+  // Remove mixed qubits from P
+  for (const std::pair<const unsigned, unsigned>& mr : mixed_to_row) {
+    unsigned Pm = (unsigned)P(mr.first) % 4;
+    if (Pm % 2 == 1) {
+      // Complementation within row of C (including the mixed qubit)
+      for (unsigned q2 = mr.first; q2 < n_qbs; ++q2) {
+        if (C(mr.second, q2)) {
+          for (unsigned q3 = q2 + 1; q3 < n_qbs; ++q3) {
+            if (C(mr.second, q3)) {
+              E(q2, q3) ^= true;
+              E(q3, q2) ^= true;
+            }
+          }
+          // Local phases (we saved to Pm earlier so we can reset P(mr.first)
+          // here)
+          P(q2) += -Pm;
+        }
+      }
+    } else if (Pm == 2) {
+      // Local phases within row of C
+      for (unsigned q2 = mr.first; q2 < n_qbs; ++q2) {
+        if (C(mr.second, q2)) {
+          P(q2) += 2;
+        }
+      }
+    }
+    // No global phase change
+  }
+}
+
+bool APState::operator==(const APState& other) const {
+  for (unsigned i = 0; i < P.size(); ++i) {
+    if (((unsigned)P(i) % 4) != ((unsigned)other.P(i) % 4)) return false;
+  }
+  return (A == other.A) && (B == other.B) && (C == other.C) && (E == other.E) &&
+         (phase == other.phase);
+}
+
+void to_json(nlohmann::json& j, const APState& aps) {
+  j["nqubits"] = aps.A.cols();
+  j["A"] = aps.A;
+  j["B"] = aps.B;
+  j["C"] = aps.C;
+  j["E"] = aps.E;
+  j["P"] = aps.P;
+  j["phase"] = aps.phase;
+}
+
+void from_json(const nlohmann::json& j, APState& aps) {
+  unsigned n_qbs = j.at("nqubits").get<unsigned>();
+  MatrixXb A(n_qbs, n_qbs);
+  VectorXb B(n_qbs);
+  MatrixXb C(n_qbs, n_qbs);
+  MatrixXb E(n_qbs, n_qbs);
+  Eigen::VectorXi P(n_qbs);
+  from_json(j.at("A"), A);
+  from_json(j.at("B"), B);
+  from_json(j.at("C"), C);
+  from_json(j.at("E"), E);
+  from_json(j.at("P"), P);
+  Expr phase = j.at("phase").get<Expr>();
+  aps = APState(A, B, C, E, P, phase);
+}
+
+/************************* ChoiAPState IMPLEMENTATION *************************/
+
+static APState id_aps(unsigned n) {
+  MatrixXb A(2 * n, 2 * n);
+  A << MatrixXb::Identity(n, n), MatrixXb::Identity(n, n),
+      MatrixXb::Zero(n, 2 * n);
+  return APState(
+      A, VectorXb::Zero(2 * n), MatrixXb::Zero(2 * n, 2 * n),
+      MatrixXb::Zero(2 * n, 2 * n), Eigen::VectorXi::Zero(2 * n), 0);
+}
+
+ChoiAPState::ChoiAPState(unsigned n) : ap_(id_aps(n)), col_index_() {
+  for (unsigned i = 0; i < n; ++i) {
+    col_index_.insert({{Qubit(i), TableauSegment::Input}, i});
+    col_index_.insert({{Qubit(i), TableauSegment::Output}, n + i});
+  }
+}
+
+ChoiAPState::ChoiAPState(const qubit_vector_t& qbs)
+    : ap_(id_aps(qbs.size())), col_index_() {
+  unsigned n = qbs.size();
+  unsigned i = 0;
+  for (const Qubit& qb : qbs) {
+    col_index_.insert({{qb, TableauSegment::Input}, i});
+    col_index_.insert({{qb, TableauSegment::Output}, n + i});
+    ++i;
+  }
+}
+
+ChoiAPState::ChoiAPState(
+    const MatrixXb& A, const VectorXb& B, const MatrixXb& C, const MatrixXb& E,
+    const Eigen::VectorXi& P, const Expr& phase, unsigned n_ins)
+    : ap_(A, B, C, E, P, phase), col_index_() {
+  unsigned n_qbs = A.cols();
+  if (n_ins > n_qbs)
+    throw std::invalid_argument(
+        "Number of inputs of a ChoiAPState cannot be larger than the number of "
+        "qubits");
+  for (unsigned i = 0; i < n_ins; ++i) {
+    col_index_.insert({{Qubit(i), TableauSegment::Input}, i});
+  }
+  for (unsigned i = 0; i < n_qbs - n_ins; ++i) {
+    col_index_.insert({{Qubit(i), TableauSegment::Output}, n_ins + i});
+  }
+}
+
+unsigned ChoiAPState::get_n_boundaries() const { return ap_.A.cols(); }
+
+unsigned ChoiAPState::get_n_inputs() const { return input_qubits().size(); }
+
+unsigned ChoiAPState::get_n_outputs() const { return output_qubits().size(); }
+
+qubit_vector_t ChoiAPState::input_qubits() const {
+  qubit_vector_t ins;
+  BOOST_FOREACH (
+      tableau_col_index_t::left_const_reference entry, col_index_.left) {
+    if (entry.first.second == TableauSegment::Input)
+      ins.push_back(entry.first.first);
+  }
+  return ins;
+}
+
+qubit_vector_t ChoiAPState::output_qubits() const {
+  qubit_vector_t outs;
+  BOOST_FOREACH (
+      tableau_col_index_t::left_const_reference entry, col_index_.left) {
+    if (entry.first.second == TableauSegment::Output)
+      outs.push_back(entry.first.first);
+  }
+  return outs;
+}
+
+void ChoiAPState::apply_gate(
+    OpType type, const qubit_vector_t& qbs, TableauSegment seg) {
+  std::vector<unsigned> u_qbs;
+  for (const Qubit& q : qbs) {
+    u_qbs.push_back(col_index_.left.at(col_key_t{q, seg}));
+  }
+  if (seg == TableauSegment::Output) {
+    ap_.apply_gate(type, u_qbs);
+  }
+  if (seg == TableauSegment::Input) {
+    switch (type) {
+      case OpType::Z:
+      case OpType::X:
+      case OpType::S:
+      case OpType::Sdg:
+      case OpType::SX:
+      case OpType::V:
+      case OpType::SXdg:
+      case OpType::Vdg:
+      case OpType::H:
+      case OpType::CX:
+      case OpType::CZ:
+      case OpType::ZZMax:
+      case OpType::ISWAPMax:
+      case OpType::SWAP:
+      case OpType::BRIDGE:
+      case OpType::noop: {
+        ap_.apply_gate(type, u_qbs);
+        break;
+      }
+      case OpType::Y: {
+        ap_.apply_gate(OpType::Y, u_qbs);
+        ap_.phase += 1.;
+        break;
+      }
+      case OpType::CY: {
+        ap_.apply_V(u_qbs.at(1));
+        ap_.apply_X(u_qbs.at(1));
+        ap_.apply_CZ(u_qbs.at(0), u_qbs.at(1));
+        ap_.apply_V(u_qbs.at(1));
+        break;
+      }
+      case OpType::ECR: {
+        ap_.apply_S(u_qbs.at(1));
+        ap_.apply_V(u_qbs.at(1));
+        ap_.apply_S(u_qbs.at(1));
+        ap_.apply_CZ(u_qbs.at(0), u_qbs.at(1));
+        ap_.apply_V(u_qbs.at(1));
+        ap_.apply_S(u_qbs.at(1));
+        ap_.apply_X(u_qbs.at(0));
+        ap_.apply_S(u_qbs.at(0));
+        ap_.phase += .25;
+        break;
+      }
+      case OpType::Reset: {
+        unsigned q = u_qbs.at(0);
+        unsigned removed_q = ap_.post_select(q);
+        tableau_col_index_t::right_iterator it =
+            col_index_.right.find(removed_q);
+        col_key_t removed_key = it->second;
+        col_index_.right.erase(it);
+        it = col_index_.right.find(u_qbs.at(0));
+        col_index_.right.erase(it);
+        col_index_.insert({removed_key, q});
+        unsigned new_q = ap_.init_qubit();
+        ap_.apply_V(new_q);
+        ap_.collapse_qubit(new_q);
+        col_index_.insert({{qbs.at(0), seg}, new_q});
+        break;
+      }
+      case OpType::Phase: {
+        throw std::logic_error(
+            "OpType::Phase cannot be applied via ChoiAPState::apply_gate");
+      }
+      default: {
+        throw BadOpType(
+            "Cannot be applied to a ChoiAPState: not a Clifford gate", type);
+      }
+    }
+  }
+}
+
+void ChoiAPState::init_qubit(const Qubit& qb, TableauSegment seg) {
+  unsigned u_qb = ap_.init_qubit();
+  col_index_.insert({{qb, seg}, u_qb});
+}
+
+void ChoiAPState::post_select(const Qubit& qb, TableauSegment seg) {
+  tableau_col_index_t::left_iterator lit =
+      col_index_.left.find(col_key_t{qb, seg});
+  unsigned u_qb = lit->second;
+  unsigned removed_u_qb = ap_.post_select(u_qb);
+  col_index_.left.erase(lit);
+  tableau_col_index_t::right_iterator rit = col_index_.right.find(removed_u_qb);
+  if (rit != col_index_.right.end()) {
+    // Ignore if post-selected last column so none others need changing
+    col_key_t removed_key = rit->second;
+    col_index_.right.erase(rit);
+    col_index_.insert({removed_key, u_qb});
+  }
+}
+
+void ChoiAPState::discard_qubit(const Qubit& qb, TableauSegment seg) {
+  collapse_qubit(qb, seg);
+  apply_gate(OpType::V, {qb}, seg);
+  collapse_qubit(qb, seg);
+  post_select(qb, seg);
+}
+
+void ChoiAPState::collapse_qubit(const Qubit& qb, TableauSegment seg) {
+  unsigned u_qb = col_index_.left.at(col_key_t{qb, seg});
+  ap_.collapse_qubit(u_qb);
+}
+
+void ChoiAPState::canonical_column_order(TableauSegment first) {
+  std::set<Qubit> ins;
+  std::set<Qubit> outs;
+  BOOST_FOREACH (
+      tableau_col_index_t::left_const_reference entry, col_index_.left) {
+    if (entry.first.second == TableauSegment::Input)
+      ins.insert(entry.first.first);
+    else
+      outs.insert(entry.first.first);
+  }
+  tableau_col_index_t new_index;
+  unsigned i = 0;
+  if (first == TableauSegment::Input) {
+    for (const Qubit& q : ins) {
+      new_index.insert({{q, TableauSegment::Input}, i});
+      ++i;
+    }
+  }
+  for (const Qubit& q : outs) {
+    new_index.insert({{q, TableauSegment::Output}, i});
+    ++i;
+  }
+  if (first == TableauSegment::Output) {
+    for (const Qubit& q : ins) {
+      new_index.insert({{q, TableauSegment::Input}, i});
+      ++i;
+    }
+  }
+  MatrixXb A = MatrixXb::Zero(i, i);
+  MatrixXb C = MatrixXb::Zero(i, i);
+  // E requires reordering both columns and rows
+  // Reorder columns to get Etemp, then reorder rows for E
+  MatrixXb Etemp = MatrixXb::Zero(i, i);
+  MatrixXb E = MatrixXb::Zero(i, i);
+  Eigen::VectorXi P = Eigen::VectorXi(i);
+  for (unsigned j = 0; j < i; ++j) {
+    col_key_t key = new_index.right.at(j);
+    unsigned c = col_index_.left.at(key);
+    A.col(j) = ap_.A.col(c);
+    C.col(j) = ap_.C.col(c);
+    Etemp.col(j) = ap_.E.col(c);
+    P(j) = ap_.P(c);
+  }
+  for (unsigned j = 0; j < i; ++j) {
+    col_key_t key = new_index.right.at(j);
+    unsigned c = col_index_.left.at(key);
+    E.row(j) = Etemp.row(c);
+  }
+  ap_ = APState(A, ap_.B, C, E, P, ap_.phase);
+  col_index_ = new_index;
+}
+
+void ChoiAPState::normal_form() { ap_.normal_form(); }
+
+void ChoiAPState::rename_qubits(const qubit_map_t& qmap, TableauSegment seg) {
+  tableau_col_index_t new_index;
+  BOOST_FOREACH (
+      tableau_col_index_t::left_const_reference entry, col_index_.left) {
+    auto found = qmap.find(entry.first.first);
+    if (entry.first.second == seg && found != qmap.end())
+      new_index.insert({{found->second, seg}, entry.second});
+    else
+      new_index.insert({entry.first, entry.second});
+  }
+  col_index_ = new_index;
+}
+
+bool ChoiAPState::operator==(const ChoiAPState& other) const {
+  return (col_index_ == other.col_index_) && (ap_ == other.ap_);
+}
+
+void to_json(nlohmann::json& j, const ChoiAPState::TableauSegment& seg) {
+  j = (seg == ChoiAPState::TableauSegment::Input) ? "In" : "Out";
+}
+
+void from_json(const nlohmann::json& j, ChoiAPState::TableauSegment& seg) {
+  const std::string str_seg = j.get<std::string>();
+  seg = (str_seg == "In") ? ChoiAPState::TableauSegment::Input
+                          : ChoiAPState::TableauSegment::Output;
+}
+
+void to_json(nlohmann::json& j, const ChoiAPState& cap) {
+  j["aps"] = cap.ap_;
+  std::vector<ChoiAPState::col_key_t> qbs;
+  for (unsigned i = 0; i < cap.get_n_boundaries(); ++i) {
+    qbs.push_back(cap.col_index_.right.at(i));
+  }
+  j["qubits"] = qbs;
+}
+
+void from_json(const nlohmann::json& j, ChoiAPState& cap) {
+  j.at("aps").get_to(cap.ap_);
+  std::vector<ChoiAPState::col_key_t> qbs =
+      j.at("qubits").get<std::vector<ChoiAPState::col_key_t>>();
+  if ((unsigned)qbs.size() != (unsigned)cap.ap_.A.cols())
+    throw std::invalid_argument(
+        "Number of qubits in json ChoiAPState does not match APState "
+        "size.");
+  cap.col_index_.clear();
+  for (unsigned i = 0; i < qbs.size(); ++i) {
+    cap.col_index_.insert({qbs.at(i), i});
+  }
+}
+
+}  // namespace tket
\ No newline at end of file
diff --git a/tket/src/Clifford/ChoiMixTableau.cpp b/tket/src/Clifford/ChoiMixTableau.cpp
index ef0b08410e..9a1b1b6271 100644
--- a/tket/src/Clifford/ChoiMixTableau.cpp
+++ b/tket/src/Clifford/ChoiMixTableau.cpp
@@ -71,7 +71,7 @@ ChoiMixTableau::ChoiMixTableau(
     col_index_.insert({{Qubit(i), TableauSegment::Input}, i});
   }
   for (unsigned i = 0; i < n_bounds - n_ins; ++i) {
-    col_index_.insert({{Qubit(i), TableauSegment::Output}, i});
+    col_index_.insert({{Qubit(i), TableauSegment::Output}, n_ins + i});
   }
 }
 
diff --git a/tket/src/Converters/APStateConverters.cpp b/tket/src/Converters/APStateConverters.cpp
new file mode 100644
index 0000000000..d8070ff35c
--- /dev/null
+++ b/tket/src/Converters/APStateConverters.cpp
@@ -0,0 +1,437 @@
+// Copyright 2019-2024 Cambridge Quantum Computing
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "tket/Converters/Converters.hpp"
+
+namespace tket {
+
+APState circuit_to_apstate(const Circuit& circ) {
+  APState aps(circ.n_qubits());
+  std::map<UnitID, unsigned> qb_ordering;
+  for (const Qubit& q : circ.all_qubits())
+    qb_ordering.insert({q, (unsigned)qb_ordering.size()});
+  for (const Command& com : circ) {
+    auto args = com.get_args();
+    std::vector<unsigned> qbs;
+    for (const UnitID& q : args) qbs.push_back(qb_ordering.at(q));
+    aps.apply_gate(com.get_op_ptr()->get_type(), qbs);
+  }
+  return aps;
+}
+
+Circuit apstate_to_circuit(const APState& ap) {
+  unsigned n_qbs = ap.A.cols();
+  Circuit circ(n_qbs);
+  circ.qubit_create_all();
+
+  MatrixXb ABgauss(n_qbs, n_qbs + 1);
+  ABgauss.block(0, 0, n_qbs, n_qbs) = ap.A;
+  ABgauss.col(n_qbs) = ap.B;
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(ABgauss);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs + 1; ++j) {
+      ABgauss(op.second, j) ^= ABgauss(op.first, j);
+    }
+  }
+  std::map<unsigned, unsigned> leader_to_row;
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    bool empty_row = true;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (ABgauss(r, c)) {
+        leader_to_row.insert({c, r});
+        empty_row = false;
+        break;
+      }
+    }
+    if (empty_row) break;
+  }
+
+  for (unsigned q = 0; q < n_qbs; ++q) {
+    auto lr = leader_to_row.find(q);
+    if (lr == leader_to_row.end())
+      circ.add_op<unsigned>(OpType::H, {q});
+    else if (ABgauss(lr->second, n_qbs))
+      circ.add_op<unsigned>(OpType::X, {q});
+  }
+  for (const std::pair<const unsigned, unsigned>& lr : leader_to_row) {
+    for (unsigned ctrl = lr.first + 1; ctrl < n_qbs; ++ctrl) {
+      if (ABgauss(lr.second, ctrl))
+        circ.add_op<unsigned>(OpType::CX, {ctrl, lr.first});
+    }
+  }
+  for (unsigned d = 0; d < n_qbs; ++d) {
+    std::optional<unsigned> first = std::nullopt;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (ap.C(d, c)) {
+        if (first) {
+          circ.add_op<unsigned>(OpType::CX, {c, *first});
+        } else {
+          first = c;
+        }
+      }
+    }
+    if (first) {
+      circ.add_op<unsigned>(OpType::Collapse, {*first});
+      for (unsigned c = n_qbs; c > *first + 1;) {
+        --c;
+        if (ap.C(d, c)) circ.add_op<unsigned>(OpType::CX, {c, *first});
+      }
+    }
+  }
+  for (unsigned q1 = 0; q1 < n_qbs; ++q1) {
+    for (unsigned q2 = q1 + 1; q2 < n_qbs; ++q2) {
+      if (ap.E(q1, q2)) circ.add_op<unsigned>(OpType::CZ, {q1, q2});
+    }
+    switch ((unsigned)ap.P(q1) % 4) {
+      case 1: {
+        circ.add_op<unsigned>(OpType::S, {q1});
+        break;
+      }
+      case 2: {
+        circ.add_op<unsigned>(OpType::Z, {q1});
+        break;
+      }
+      case 3: {
+        circ.add_op<unsigned>(OpType::Sdg, {q1});
+        break;
+      }
+      default: {
+        break;
+      }
+    }
+  }
+  circ.add_phase(ap.phase);
+  return circ;
+}
+
+ChoiAPState circuit_to_choi_apstate(const Circuit& circ) {
+  ChoiAPState ap(circ.all_qubits());
+  for (const Qubit& q : circ.created_qubits()) {
+    ap.post_select(q, ChoiAPState::TableauSegment::Input);
+  }
+  for (const Command& com : circ) {
+    auto args = com.get_args();
+    qubit_vector_t qbs = {args.begin(), args.end()};
+    ap.apply_gate(com.get_op_ptr()->get_type(), qbs);
+  }
+  ap.rename_qubits(
+      circ.implicit_qubit_permutation(), ChoiAPState::TableauSegment::Output);
+  for (const Qubit& q : circ.discarded_qubits()) {
+    ap.discard_qubit(q);
+  }
+  ap.canonical_column_order();
+  ap.normal_form();
+  return ap;
+}
+
+std::pair<Circuit, qubit_map_t> choi_apstate_to_exact_circuit(
+    ChoiAPState ap, CXConfigType cx_config) {
+  ChoiMixTableau tab = choi_apstate_to_cm_tableau(ap);
+  std::pair<Circuit, qubit_map_t> res =
+      cm_tableau_to_exact_circuit(tab, cx_config);
+  ChoiAPState reconstructed = circuit_to_choi_apstate(res.first);
+  reconstructed.rename_qubits(res.second, ChoiAPState::TableauSegment::Output);
+  ap.canonical_column_order();
+  ap.normal_form();
+  reconstructed.canonical_column_order();
+  reconstructed.normal_form();
+  res.first.add_phase(ap.ap_.phase - reconstructed.ap_.phase);
+  return res;
+}
+
+std::pair<Circuit, qubit_map_t> choi_apstate_to_unitary_extension_circuit(
+    ChoiAPState ap, const std::vector<Qubit>& init_names,
+    const std::vector<Qubit>& post_names, CXConfigType cx_config) {
+  ChoiMixTableau tab = choi_apstate_to_cm_tableau(ap);
+  std::pair<Circuit, qubit_map_t> res = cm_tableau_to_unitary_extension_circuit(
+      tab, init_names, post_names, cx_config);
+  ChoiAPState reconstructed = circuit_to_choi_apstate(res.first);
+  for (const Qubit& q : init_names)
+    reconstructed.post_select(q, ChoiAPState::TableauSegment::Input);
+  for (const Qubit& q : post_names)
+    reconstructed.post_select(q, ChoiAPState::TableauSegment::Output);
+  reconstructed.rename_qubits(res.second, ChoiAPState::TableauSegment::Output);
+  ap.canonical_column_order();
+  ap.normal_form();
+  reconstructed.canonical_column_order();
+  reconstructed.normal_form();
+  res.first.add_phase(ap.ap_.phase - reconstructed.ap_.phase);
+  return res;
+}
+
+APState tableau_to_apstate(SymplecticTableau tab) {
+  unsigned n_qbs = tab.get_n_qubits();
+  unsigned n_rows = tab.get_n_rows();
+  if (n_rows > n_qbs)
+    throw std::logic_error(
+        "tableau_to_apstate requires a tableau with up to n commuting rows for "
+        "n qubits");
+  MatrixXb fullmat = MatrixXb::Zero(n_rows, 2 * n_qbs);
+  /**
+   * Gaussian elimination by the x matrix first ensures the bottom rows are only
+   * Zs, i.e. describing rows of A. Reversing the columns of the x matrix
+   * guarantees that each row has an X on at most one free qubit, simplifying
+   * the code for finding E and P
+   */
+  for (unsigned c = 0; c < n_qbs; ++c) {
+    fullmat.col(c) = tab.xmat.col(n_qbs - 1 - c);
+  }
+  fullmat.block(0, n_qbs, n_rows, n_qbs) = tab.zmat;
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(fullmat);
+  for (const std::pair<unsigned, unsigned>& op : row_ops)
+    tab.row_mult(op.first, op.second);
+
+  MatrixXb A = MatrixXb::Zero(n_qbs, n_qbs);
+  VectorXb B = VectorXb::Zero(n_qbs);
+  MatrixXb C = MatrixXb::Zero(n_qbs, n_qbs);
+  MatrixXb E = MatrixXb::Zero(n_qbs, n_qbs);
+  Eigen::VectorXi P = Eigen::VectorXi::Zero(n_qbs);
+
+  unsigned n_leading = 0;
+  for (unsigned r = n_rows; r > 0;) {
+    --r;
+    bool is_a_row = true;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (tab.xmat(r, c)) {
+        is_a_row = false;
+        break;
+      }
+    }
+    if (!is_a_row) break;
+    ++n_leading;
+  }
+  A.block(0, 0, n_leading, n_qbs) =
+      tab.zmat.block(n_rows - n_leading, 0, n_leading, n_qbs);
+  B.head(n_leading) = tab.phase.tail(n_leading);
+  std::map<unsigned, unsigned> leader_to_row;
+  for (unsigned r = 0; r < n_leading; ++r) {
+    for (unsigned c = r; c < n_qbs; ++c) {
+      if (A(r, c)) {
+        leader_to_row.insert({c, r});
+        break;
+      }
+    }
+  }
+
+  /**
+   * Each free qubit q is after all leaders connected to it, by reduced
+   * row-echelon of A. Therefore Gaussian elimination of the x matrix in reverse
+   * order gives rows corresponding to the columns of A of free qubits (plus the
+   * free qubit itself).
+   *
+   * Then the corresponding row q of the z matrix is the free qubit's row of E,
+   * plus the rows of any E for any mixed qubit connected to q in C, plus an
+   * extra flip at q if P(q) is odd. We can first look at the column for each
+   * mixed qubit to identify their rows of E and subtract from this matrix to
+   * leave the interactions between free qubits and their local phases.
+   */
+
+  // Start by identifying the free and mixed qubits
+  std::map<unsigned, unsigned> free_to_row;
+  for (unsigned r = 0; r < n_rows - n_leading; ++r) {
+    for (unsigned c = n_qbs - r; c > 0;) {
+      --c;
+      if (tab.xmat(r, c)) {
+        free_to_row.insert({c, r});
+        break;
+      }
+    }
+  }
+  std::map<unsigned, unsigned> mixed_to_row;
+  for (unsigned q = 0; q < n_qbs; ++q) {
+    if (leader_to_row.find(q) == leader_to_row.end() &&
+        free_to_row.find(q) == free_to_row.end()) {
+      unsigned r = mixed_to_row.size();
+      mixed_to_row.insert({q, r});
+      C(r, q) = true;
+    }
+  }
+  // Identify C and the mixed rows of E by looking at what mixed qubits appear
+  // in the rows of each free qubit
+  for (const std::pair<const unsigned, unsigned>& fr : free_to_row) {
+    for (const std::pair<const unsigned, unsigned>& mr : mixed_to_row) {
+      if (tab.xmat(fr.second, mr.first)) C(mr.second, fr.first) = true;
+      if (tab.zmat(fr.second, mr.first)) {
+        E(mr.first, fr.first) = true;
+        E(fr.first, mr.first) = true;
+      }
+    }
+  }
+  // Identify connections in E between free qubits
+  for (auto fr1 = free_to_row.begin(); fr1 != free_to_row.end(); ++fr1) {
+    for (auto fr2 = free_to_row.begin(); fr2 != fr1; ++fr2) {
+      unsigned n_shared_mixed = 0;
+      for (const std::pair<const unsigned, unsigned>& mr : mixed_to_row) {
+        if (C(mr.second, fr1->first) && C(mr.second, fr2->first))
+          ++n_shared_mixed;
+      }
+      if (tab.zmat(fr1->second, fr2->first) ^ (n_shared_mixed % 2 == 1)) {
+        E(fr1->first, fr2->first) = true;
+        E(fr2->first, fr1->first) = true;
+      }
+    }
+  }
+  // Identify P
+  for (const std::pair<const unsigned, unsigned>& fr : free_to_row) {
+    unsigned n_mixed_in_c_and_e = 0;
+    for (const std::pair<const unsigned, unsigned>& mr : mixed_to_row) {
+      if (C(mr.second, fr.first) && E(mr.first, fr.first)) ++n_mixed_in_c_and_e;
+    }
+    if (tab.zmat(fr.second, fr.first) ^ (n_mixed_in_c_and_e % 2 == 1))
+      P(fr.first) += 1;
+    if (tab.phase(fr.second) ^ (n_mixed_in_c_and_e % 2 == 1)) P(fr.first) += 2;
+  }
+
+  return APState(A, B, C, E, P, 0);
+}
+
+SymplecticTableau apstate_to_tableau(APState ap) {
+  unsigned n_qbs = ap.A.cols();
+  // Want A and C in reduced row-echelon form to identify leaders and mixed
+  // qubits, but don't need the rest in normal form
+  std::vector<std::pair<unsigned, unsigned>> row_ops =
+      gaussian_elimination_row_ops(ap.A);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      ap.A(op.second, j) ^= ap.A(op.first, j);
+    }
+    ap.B(op.second) ^= ap.B(op.first);
+  }
+  std::map<unsigned, unsigned> leader_to_row;
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    bool leader_found = false;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (ap.A(r, c)) {
+        leader_found = true;
+        leader_to_row.insert({c, r});
+        break;
+      }
+    }
+    if (!leader_found) break;
+  }
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    for (const std::pair<const unsigned, unsigned>& lr : leader_to_row) {
+      if (ap.C(r, lr.first)) {
+        for (unsigned c = 0; c < n_qbs; ++c) {
+          ap.C(r, c) ^= ap.A(lr.second, c);
+        }
+      }
+    }
+  }
+  row_ops = gaussian_elimination_row_ops(ap.C);
+  for (const std::pair<unsigned, unsigned>& op : row_ops) {
+    for (unsigned j = 0; j < n_qbs; ++j) {
+      ap.C(op.second, j) ^= ap.C(op.first, j);
+    }
+  }
+  std::map<unsigned, unsigned> mixed_to_row;
+  for (unsigned r = 0; r < n_qbs; ++r) {
+    bool mixed_found = false;
+    for (unsigned c = 0; c < n_qbs; ++c) {
+      if (ap.C(r, c)) {
+        mixed_found = true;
+        mixed_to_row.insert({c, r});
+        break;
+      }
+    }
+    if (!mixed_found) break;
+  }
+  unsigned n_rows = n_qbs - mixed_to_row.size();
+  MatrixXb xmat = MatrixXb::Zero(n_rows, n_qbs);
+  MatrixXb zmat = MatrixXb::Zero(n_rows, n_qbs);
+  VectorXb phase = VectorXb::Zero(n_rows);
+
+  // One stabiliser per leader, with Z on that qubit and Z on every neighbour in
+  // A
+  unsigned n_leading = leader_to_row.size();
+  zmat.block(0, 0, n_leading, n_qbs) = ap.A.block(0, 0, n_leading, n_qbs);
+  phase.head(n_leading) = ap.B.head(n_leading);
+
+  // One stabiliser per free, with X on that qubit, every neighbour in C, and
+  // the odd neighbourhood of this set in A; pushing this through the phase
+  // polynomial adds Zs
+  unsigned r = n_leading;
+  for (unsigned q = 0; q < n_qbs; ++q) {
+    if (leader_to_row.find(q) == leader_to_row.end() &&
+        mixed_to_row.find(q) == mixed_to_row.end()) {
+      // Calculate the Xs
+      xmat(r, q) = true;
+      for (const std::pair<const unsigned, unsigned>& mr : mixed_to_row) {
+        xmat(r, mr.first) = ap.C(mr.second, q);
+      }
+      for (const std::pair<const unsigned, unsigned>& lr : leader_to_row) {
+        unsigned n_mixed_in_between = 0;
+        for (const std::pair<const unsigned, unsigned>& mr : mixed_to_row) {
+          if (ap.A(lr.second, mr.first) && ap.C(mr.second, q))
+            ++n_mixed_in_between;
+        }
+        xmat(r, lr.first) = ap.A(lr.second, q) ^ (n_mixed_in_between % 2 == 1);
+      }
+      // Push through the phase polynomial to calculate the Zs
+      for (unsigned q2 = 0; q2 < n_qbs; ++q2) {
+        if (xmat(r, q2)) {
+          // Pushing an X through each CZ creates a Z
+          for (unsigned q3 = 0; q3 < n_qbs; ++q3) {
+            if (ap.E(q2, q3)) {
+              zmat(r, q3) ^= true;
+              // If both qubits of a CZ have Xs, we will need to reorder the X
+              // and Z on one to match the other
+              if (q2 < q3 && xmat(r, q3)) phase(r) ^= true;
+            }
+          }
+          // Pushing an X through the local phase
+          zmat(r, q2) ^= ((unsigned)ap.P(q2) % 2 == 1);
+          phase(r) ^= ((unsigned)ap.P(q2) % 4 > 1);
+        }
+      }
+      ++r;
+    }
+  }
+
+  return SymplecticTableau(xmat, zmat, phase);
+}
+
+ChoiAPState cm_tableau_to_choi_apstate(const ChoiMixTableau& tab) {
+  ChoiAPState ap(0);
+  ap.ap_ = tableau_to_apstate(tab.tab_);
+  BOOST_FOREACH (
+      ChoiMixTableau::tableau_col_index_t::left_const_reference entry,
+      tab.col_index_.left) {
+    ChoiAPState::TableauSegment seg =
+        (entry.first.second == ChoiMixTableau::TableauSegment::Input)
+            ? ChoiAPState::TableauSegment::Input
+            : ChoiAPState::TableauSegment::Output;
+    ap.col_index_.insert({{entry.first.first, seg}, entry.second});
+  }
+  return ap;
+}
+
+ChoiMixTableau choi_apstate_to_cm_tableau(const ChoiAPState& ap) {
+  ChoiMixTableau tab(0);
+  tab.tab_ = apstate_to_tableau(ap.ap_);
+  BOOST_FOREACH (
+      ChoiAPState::tableau_col_index_t::left_const_reference entry,
+      ap.col_index_.left) {
+    ChoiMixTableau::TableauSegment seg =
+        (entry.first.second == ChoiAPState::TableauSegment::Input)
+            ? ChoiMixTableau::TableauSegment::Input
+            : ChoiMixTableau::TableauSegment::Output;
+    tab.col_index_.insert({{entry.first.first, seg}, entry.second});
+  }
+  return tab;
+}
+
+}  // namespace tket
\ No newline at end of file
diff --git a/tket/test/CMakeLists.txt b/tket/test/CMakeLists.txt
index 44437e6df2..cbaefa17cb 100644
--- a/tket/test/CMakeLists.txt
+++ b/tket/test/CMakeLists.txt
@@ -103,6 +103,7 @@ add_executable(test-tket
     src/Simulation/test_CircuitSimulator.cpp
     src/Simulation/test_PauliExpBoxUnitaryCalculator.cpp
     src/test_AASRoute.cpp
+    src/test_APState.cpp
     src/test_ArchitectureAwareSynthesis.cpp
     src/test_Architectures.cpp
     src/test_Assertion.cpp
diff --git a/tket/test/src/test_APState.cpp b/tket/test/src/test_APState.cpp
new file mode 100644
index 0000000000..064c3830df
--- /dev/null
+++ b/tket/test/src/test_APState.cpp
@@ -0,0 +1,844 @@
+// Copyright 2019-2024 Cambridge Quantum Computing
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include <catch2/catch_test_macros.hpp>
+
+#include "Simulation/ComparisonFunctions.hpp"
+#include "testutil.hpp"
+#include "tket/Circuit/Simulation/CircuitSimulator.hpp"
+#include "tket/Clifford/APState.hpp"
+#include "tket/Converters/Converters.hpp"
+
+namespace tket {
+namespace test_APState {
+
+void test_apply_gate(
+    const MatrixXb& A, const VectorXb& B, const MatrixXb& E,
+    const Eigen::VectorXi& P, OpType ot, const std::vector<unsigned>& args) {
+  MatrixXb C = MatrixXb::Zero(A.cols(), A.cols());
+  APState ap(A, B, C, E, P, 0);
+  ap.verify();
+  auto sv_before = ap.to_statevector();
+
+  Circuit circ(A.cols());
+  circ.add_op<unsigned>(ot, args);
+  auto gate_u = tket_sim::get_unitary(circ);
+
+  ap.apply_gate(ot, args);
+
+  ap.verify();
+  auto sv_after = ap.to_statevector();
+
+  CHECK(tket_sim::compare_statevectors_or_unitaries(
+      gate_u * sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL));
+}
+
+void test_apply_gate_dm(
+    const MatrixXb& A, const VectorXb& B, const MatrixXb& C, const MatrixXb& E,
+    const Eigen::VectorXi& P, OpType ot, const std::vector<unsigned>& args) {
+  APState ap(A, B, C, E, P, 0);
+  ap.verify();
+  auto dm_before = ap.to_density_matrix();
+
+  Circuit circ(A.cols());
+  circ.add_op<unsigned>(ot, args);
+  auto gate_u = tket_sim::get_unitary(circ);
+
+  ap.apply_gate(ot, args);
+
+  ap.verify();
+  auto dm_after = ap.to_density_matrix();
+
+  CHECK(dm_after.isApprox(gate_u * dm_before * gate_u.adjoint(), EPS));
+}
+
+void test_apply_gate_dm_input(
+    const MatrixXb& A, const VectorXb& B, const MatrixXb& C, const MatrixXb& E,
+    const Eigen::VectorXi& P, OpType ot, const qubit_vector_t& args) {
+  ChoiAPState ap(A, B, C, E, P, 0, A.cols());
+  ap.ap_.verify();
+  auto dm_before = ap.ap_.to_density_matrix();
+
+  Circuit circ(A.cols());
+  circ.add_op<Qubit>(ot, args);
+  auto gate_u = tket_sim::get_unitary(circ);
+
+  ap.apply_gate(ot, args, ChoiAPState::TableauSegment::Input);
+
+  ap.ap_.verify();
+  auto dm_after = ap.ap_.to_density_matrix();
+
+  CHECK(dm_after.isApprox(
+      gate_u.transpose() * dm_before * gate_u.conjugate(), EPS));
+}
+
+SCENARIO("Normal form") {
+  GIVEN("Make A reduced row-echelon form (pure state)") {
+    MatrixXb A = MatrixXb::Zero(4, 4);
+    VectorXb B = VectorXb::Zero(4);
+    MatrixXb C = MatrixXb::Zero(4, 4);
+    MatrixXb E = MatrixXb::Zero(4, 4);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+    A(0, 2) = A(0, 3) = true;
+    A(1, 0) = A(1, 1) = A(1, 2) = true;
+    A(2, 0) = A(2, 1) = A(2, 2) = true;
+    A(3, 0) = A(3, 1) = true;
+    APState ap(A, B, C, E, P, 0);
+    auto sv_before = ap.to_statevector();
+    ap.normal_form();
+    auto sv_after = ap.to_statevector();
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL));
+    MatrixXb corrA = MatrixXb::Zero(4, 4);
+    corrA(0, 0) = corrA(0, 1) = true;
+    corrA(1, 2) = true;
+    corrA(2, 3) = true;
+    CHECK(ap.A == corrA);
+  }
+  GIVEN("Make A and C reduced row-echelon form") {
+    MatrixXb A = MatrixXb::Zero(4, 4);
+    VectorXb B = VectorXb::Zero(4);
+    MatrixXb C = MatrixXb::Zero(4, 4);
+    MatrixXb E = MatrixXb::Zero(4, 4);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+    A(0, 2) = A(0, 3) = true;
+    A(1, 0) = A(1, 1) = A(1, 2) = true;
+    C(0, 0) = C(0, 1) = C(0, 2) = true;
+    C(1, 0) = true;
+    C(3, 2) = true;
+    APState ap(A, B, C, E, P, 0);
+    auto dm_before = ap.to_density_matrix();
+    ap.normal_form();
+    auto dm_after = ap.to_density_matrix();
+    CHECK(dm_after.isApprox(dm_before, EPS));
+    MatrixXb corrA = MatrixXb::Zero(4, 4);
+    MatrixXb corrC = MatrixXb::Zero(4, 4);
+    corrA(0, 0) = corrA(0, 1) = corrA(0, 3) = true;
+    corrA(1, 2) = corrA(1, 3) = true;
+    corrC(0, 1) = true;
+    corrC(1, 3) = true;
+    CHECK(ap.A == corrA);
+    CHECK(ap.C == corrC);
+  }
+  GIVEN("Removing leaders from E and P (pure state)") {
+    MatrixXb A = MatrixXb::Zero(5, 5);
+    VectorXb B = VectorXb::Zero(5);
+    MatrixXb C = MatrixXb::Zero(5, 5);
+    MatrixXb E = MatrixXb::Zero(5, 5);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(5);
+    A(0, 0) = A(0, 2) = A(0, 3) = true;
+    A(1, 1) = A(1, 2) = true;
+    E(0, 1) = E(1, 0) = true;
+    E(0, 3) = E(3, 0) = true;
+    E(0, 4) = E(4, 0) = true;
+    for (unsigned b = 0; b < 2; ++b) {
+      for (unsigned p = 0; p < 4; ++p) {
+        B(0) = (b == 1);
+        P(0) = p;
+        APState ap(A, B, C, E, P, 0);
+        auto sv_before = ap.to_statevector();
+        ap.normal_form();
+        auto sv_after = ap.to_statevector();
+        // Just check using statevector; too much changes in each case to nicely
+        // test the matrices
+        CHECK(tket_sim::compare_statevectors_or_unitaries(
+            sv_before, sv_after, tket_sim::MatrixEquivalence::EQUAL));
+      }
+    }
+  }
+  GIVEN("Removing mixed qubits from E and P") {
+    MatrixXb A = MatrixXb::Zero(5, 5);
+    VectorXb B = VectorXb::Zero(5);
+    MatrixXb C = MatrixXb::Zero(5, 5);
+    MatrixXb E = MatrixXb::Zero(5, 5);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(5);
+    C(0, 0) = C(0, 2) = C(0, 3) = true;
+    C(1, 1) = true;
+    E(0, 1) = E(1, 0) = true;
+    E(0, 2) = E(2, 0) = true;
+    E(0, 4) = E(4, 0) = true;
+    for (unsigned p = 0; p < 4; ++p) {
+      P(0) = p;
+      APState ap(A, B, C, E, P, 0);
+      auto dm_before = ap.to_density_matrix();
+      ap.normal_form();
+      auto dm_after = ap.to_density_matrix();
+      // Just check using statevector; too much changes in each case to nicely
+      // test the matrices
+      CHECK(dm_after.isApprox(dm_before, EPS));
+    }
+  }
+}
+
+SCENARIO("CZ cases") {
+  GIVEN("CZ on free qubits") {
+    MatrixXb A = MatrixXb::Zero(4, 4);
+    VectorXb B = VectorXb::Zero(4);
+    MatrixXb E = MatrixXb::Zero(4, 4);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+    A(0, 0) = A(0, 1) = A(0, 3) = true;
+    E(2, 3) = E(3, 2) = true;
+    test_apply_gate(A, B, E, P, OpType::CZ, {1, 2});
+  }
+  GIVEN("CZ on one leading qubit and connected free") {
+    for (unsigned b = 0; b < 2; ++b) {
+      MatrixXb A = MatrixXb::Zero(5, 5);
+      VectorXb B = VectorXb::Zero(5);
+      MatrixXb E = MatrixXb::Zero(5, 5);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(5);
+      A(0, 0) = A(0, 2) = A(0, 3) = true;
+      A(1, 1) = A(1, 3) = A(1, 4) = true;
+      B(1) = (b == 1);
+      test_apply_gate(A, B, E, P, OpType::CZ, {0, 3});
+    }
+  }
+  GIVEN("CZ on one leading qubit and unconnected free") {
+    for (unsigned b = 0; b < 2; ++b) {
+      MatrixXb A = MatrixXb::Zero(5, 5);
+      VectorXb B = VectorXb::Zero(5);
+      MatrixXb E = MatrixXb::Zero(5, 5);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(5);
+      A(0, 0) = A(0, 2) = A(0, 3) = true;
+      A(1, 1) = A(1, 3) = A(1, 4) = true;
+      B(1) = (b == 1);
+      test_apply_gate(A, B, E, P, OpType::CZ, {0, 4});
+    }
+  }
+  GIVEN("CZ on leading qubits") {
+    for (unsigned b1 = 0; b1 < 2; ++b1) {
+      for (unsigned b2 = 0; b2 < 2; ++b2) {
+        MatrixXb A = MatrixXb::Zero(8, 8);
+        VectorXb B = VectorXb::Zero(8);
+        MatrixXb E = MatrixXb::Zero(8, 8);
+        Eigen::VectorXi P = Eigen::VectorXi::Zero(8);
+        A(0, 0) = A(0, 2) = A(0, 3) = A(0, 4) = A(0, 5) = true;
+        A(1, 1) = A(1, 4) = A(1, 5) = A(1, 6) = A(1, 7) = true;
+        B(0) = (b1 == 1);
+        B(1) = (b2 == 1);
+        test_apply_gate(A, B, E, P, OpType::CZ, {0, 1});
+      }
+    }
+  }
+  GIVEN("CZ on mixed state") {
+    for (unsigned b1 = 0; b1 < 2; ++b1) {
+      for (unsigned b2 = 0; b2 < 2; ++b2) {
+        MatrixXb A = MatrixXb::Zero(5, 5);
+        VectorXb B = VectorXb::Zero(5);
+        MatrixXb C = MatrixXb::Zero(5, 5);
+        MatrixXb E = MatrixXb::Zero(5, 5);
+        Eigen::VectorXi P = Eigen::VectorXi::Zero(5);
+        A(0, 1) = A(0, 3) = A(0, 4) = true;
+        A(1, 0) = A(1, 2) = A(1, 3) = true;
+        B(0) = (b1 == 1);
+        B(1) = (b1 == 1);
+        C(0, 0) = C(0, 2) = true;
+        E(0, 3) = E(3, 0) = true;
+        E(1, 3) = E(3, 1) = true;
+        P(0) = 3;
+        P(1) = 1;
+        test_apply_gate_dm(A, B, C, E, P, OpType::CZ, {0, 1});
+      }
+    }
+  }
+}
+
+SCENARIO("S cases") {
+  GIVEN("S on free qubit") {
+    MatrixXb A = MatrixXb::Zero(3, 3);
+    VectorXb B = VectorXb::Zero(3);
+    MatrixXb E = MatrixXb::Zero(3, 3);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+    A(0, 0) = A(0, 1) = A(0, 2) = true;
+    E(1, 2) = E(2, 1) = true;
+    test_apply_gate(A, B, E, P, OpType::S, {2});
+  }
+  GIVEN("S on leading qubit") {
+    for (unsigned b = 0; b < 2; ++b) {
+      MatrixXb A = MatrixXb::Zero(4, 4);
+      VectorXb B = VectorXb::Zero(4);
+      MatrixXb E = MatrixXb::Zero(4, 4);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+      A(0, 0) = A(0, 1) = A(0, 2) = true;
+      B(0) = (b == 1);
+      E(1, 3) = E(3, 1) = true;
+      test_apply_gate(A, B, E, P, OpType::S, {0});
+    }
+  }
+  GIVEN("S on disconnected leading qubit") {
+    for (unsigned b = 0; b < 2; ++b) {
+      MatrixXb A = MatrixXb::Zero(1, 1);
+      VectorXb B = VectorXb::Zero(1);
+      MatrixXb E = MatrixXb::Zero(1, 1);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(1);
+      A(0, 0) = true;
+      B(0) = (b == 1);
+      test_apply_gate(A, B, E, P, OpType::S, {0});
+    }
+  }
+  GIVEN("S on a mixed state") {
+    MatrixXb A = MatrixXb::Zero(3, 3);
+    VectorXb B = VectorXb::Zero(3);
+    MatrixXb C = MatrixXb::Zero(3, 3);
+    MatrixXb E = MatrixXb::Zero(3, 3);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+    A(0, 0) = A(0, 1) = A(0, 2) = true;
+    C(0, 1) = C(0, 2) = true;
+    E(1, 2) = E(2, 1) = true;
+    test_apply_gate_dm(A, B, C, E, P, OpType::S, {2});
+  }
+}
+
+SCENARIO("V cases") {
+  GIVEN("V on leading qubit") {
+    for (unsigned b = 0; b < 2; ++b) {
+      MatrixXb A = MatrixXb::Zero(4, 4);
+      VectorXb B = VectorXb::Zero(4);
+      MatrixXb E = MatrixXb::Zero(4, 4);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+      A(0, 0) = A(0, 2) = A(0, 3) = true;
+      A(1, 1) = A(1, 3) = true;
+      B(0) = (b == 1);
+      test_apply_gate(A, B, E, P, OpType::V, {0});
+    }
+  }
+  GIVEN("V on free qubit with some leading") {
+    for (unsigned b = 0; b < 2; ++b) {
+      for (unsigned p = 0; p < 4; ++p) {
+        MatrixXb A = MatrixXb::Zero(9, 9);
+        VectorXb B = VectorXb::Zero(9);
+        MatrixXb E = MatrixXb::Zero(9, 9);
+        Eigen::VectorXi P = Eigen::VectorXi::Zero(9);
+        A(0, 0) = A(0, 2) = A(0, 4) = A(0, 5) = A(0, 7) = true;
+        A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = A(1, 5) = A(1, 6) = true;
+        B(1) = (b == 1);
+        E(4, 5) = E(5, 4) = true;
+        E(4, 6) = E(6, 4) = true;
+        E(4, 7) = E(7, 4) = true;
+        E(4, 8) = E(8, 4) = true;
+        P(4) = p;
+        test_apply_gate(A, B, E, P, OpType::V, {4});
+      }
+    }
+  }
+  GIVEN("V on free qubit with some earlier connected free") {
+    for (unsigned p1 = 0; p1 < 4; ++p1) {
+      for (unsigned p2 = 0; p2 < 4; ++p2) {
+        MatrixXb A = MatrixXb::Zero(9, 9);
+        VectorXb B = VectorXb::Zero(9);
+        MatrixXb E = MatrixXb::Zero(9, 9);
+        Eigen::VectorXi P = Eigen::VectorXi::Zero(9);
+        A(0, 0) = true;
+        A(1, 1) = A(1, 4) = A(1, 6) = true;
+        A(2, 2) = A(2, 4) = A(2, 7) = true;
+        A(3, 3) = A(3, 4) = A(3, 8) = A(3, 8) = true;
+        E(4, 5) = E(5, 4) = true;
+        E(4, 7) = E(7, 4) = true;
+        E(4, 8) = E(8, 4) = true;
+        E(5, 6) = E(6, 5) = true;
+        E(5, 7) = E(7, 5) = true;
+        E(5, 8) = E(8, 5) = true;
+        P(4) = p1;
+        P(5) = p2;
+        test_apply_gate(A, B, E, P, OpType::V, {5});
+      }
+    }
+  }
+  GIVEN("V on free qubit with no earlier connected free") {
+    for (unsigned p = 0; p < 4; ++p) {
+      MatrixXb A = MatrixXb::Zero(4, 4);
+      VectorXb B = VectorXb::Zero(4);
+      MatrixXb E = MatrixXb::Zero(4, 4);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+      A(0, 0) = A(0, 2) = true;
+      E(1, 2) = E(2, 1) = true;
+      E(1, 3) = E(3, 1) = true;
+      P(1) = p;
+      test_apply_gate(A, B, E, P, OpType::V, {1});
+    }
+  }
+  GIVEN("V on disconnected free qubit") {
+    for (unsigned p = 0; p < 4; ++p) {
+      MatrixXb A = MatrixXb::Zero(1, 1);
+      VectorXb B = VectorXb::Zero(1);
+      MatrixXb E = MatrixXb::Zero(1, 1);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(1, 1);
+      P(0) = p;
+      test_apply_gate(A, B, E, P, OpType::V, {0});
+    }
+  }
+  GIVEN("V on a qubit involved in A and C") {
+    for (unsigned p = 0; p < 4; ++p) {
+      MatrixXb A = MatrixXb::Zero(7, 7);
+      VectorXb B = VectorXb::Zero(7);
+      MatrixXb C = MatrixXb::Zero(7, 7);
+      MatrixXb E = MatrixXb::Zero(7, 7);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(7);
+      A(0, 0) = A(0, 2) = true;
+      A(1, 0) = A(1, 1) = A(1, 2) = A(1, 3) = A(1, 4) = true;
+      C(0, 0) = C(0, 2) = C(0, 3) = true;
+      C(1, 0) = C(1, 1) = C(1, 2) = C(1, 5) = C(1, 6) = true;
+      E(0, 3) = E(3, 0) = true;
+      E(0, 4) = E(4, 0) = true;
+      E(0, 5) = E(5, 0) = true;
+      E(0, 6) = E(6, 0) = true;
+      P(0) = p;
+      test_apply_gate_dm(A, B, C, E, P, OpType::V, {0});
+    }
+  }
+  GIVEN("V on a mixed state with zero A") {
+    for (unsigned p = 0; p < 4; ++p) {
+      MatrixXb A = MatrixXb::Zero(7, 7);
+      VectorXb B = VectorXb::Zero(7);
+      MatrixXb C = MatrixXb::Zero(7, 7);
+      MatrixXb E = MatrixXb::Zero(7, 7);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(7);
+      C(0, 0) = C(0, 2) = true;
+      C(1, 0) = C(1, 1) = C(1, 2) = C(1, 3) = C(1, 4) = true;
+      C(2, 0) = C(2, 2) = C(0, 3) = true;
+      C(3, 0) = C(3, 1) = C(3, 2) = C(3, 5) = C(3, 6) = true;
+      E(0, 3) = E(3, 0) = true;
+      E(0, 4) = E(4, 0) = true;
+      E(0, 5) = E(5, 0) = true;
+      E(0, 6) = E(6, 0) = true;
+      P(0) = p;
+      test_apply_gate_dm(A, B, C, E, P, OpType::V, {0});
+    }
+  }
+}
+
+SCENARIO("Qubit Reset") {
+  GIVEN("Reset a qubit with a local state") {
+    // Qubit 0 is in the |0> state
+    MatrixXb A = MatrixXb::Zero(3, 3);
+    VectorXb B = VectorXb::Zero(3);
+    MatrixXb C = MatrixXb::Zero(3, 3);
+    MatrixXb E = MatrixXb::Zero(3, 3);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+    A(0, 0) = true;
+    A(1, 1) = A(1, 2) = true;
+    B(1) = true;
+    P(2) = 1;
+    APState correct(A, B, C, E, P, 0);
+    // correct is already in normal form
+    for (unsigned s = 0; s < 6; ++s) {
+      // s = 0,1,2,3: XY basis states
+      // s = 4: |0>
+      // s = 5: |1>
+      A(0, 0) = (s < 4) ? false : true;
+      P(0) = (s < 4) ? s : 0;
+      B(0) = (s == 5);
+      APState ap(A, B, C, E, P, 0);
+      ap.apply_gate(OpType::Reset, {0});
+      // Check up to global phase
+      ap.normal_form();
+      ap.phase = correct.phase;
+      CHECK(ap == correct);
+    }
+  }
+  GIVEN("Reset one side of a Bell state") {
+    MatrixXb A = MatrixXb::Zero(3, 3);
+    VectorXb B = VectorXb::Zero(3);
+    MatrixXb C = MatrixXb::Zero(3, 3);
+    MatrixXb E = MatrixXb::Zero(3, 3);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+    A(0, 0) = A(0, 1) = true;
+    APState ap(A, B, C, E, P, 0);
+    ap.apply_gate(OpType::Reset, {0});
+    ap.normal_form();
+    // Qubit 0 ends in |0>
+    A(0, 1) = false;
+    // Qubit 1 ends in maximally-mixed state
+    C(0, 1) = true;
+    APState correct(A, B, C, E, P, 0);
+    CHECK(ap == correct);
+  }
+  GIVEN("Reset on a normal form state") {
+    MatrixXb A = MatrixXb::Zero(4, 4);
+    VectorXb B = VectorXb::Zero(4);
+    MatrixXb C = MatrixXb::Zero(4, 4);
+    MatrixXb E = MatrixXb::Zero(4, 4);
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+    A(0, 0) = A(0, 1) = A(0, 3) = true;
+    B(0) = true;
+    C(0, 1) = C(0, 2) = true;
+    E(1, 3) = E(3, 1) = true;
+    E(2, 3) = E(3, 2) = true;
+    P(2) = 1;
+    P(3) = 2;
+    APState ap(A, B, C, E, P, 0);
+    WHEN("Apply to Qubit 0") {
+      ap.apply_gate(OpType::Reset, {0});
+      ap.normal_form();
+      // Qubit 0 in state |0>
+      A(0, 1) = A(0, 3) = false;
+      B(0) = false;
+      // A row becomes a C row (combine with other row for gaussian form)
+      C(1, 2) = C(1, 3) = true;
+      // More gaussian steps
+      C(0, 2) = false;
+      C(0, 3) = true;
+      // LC about C(1, -) to remove P(2)
+      E(2, 3) = E(3, 2) = false;
+      P(2) = 0;
+      P(3) = 1;
+      APState correct(A, B, C, E, P, {0});
+      // Check equality up to global phase
+      ap.phase = correct.phase;
+      CHECK(ap == correct);
+    }
+    WHEN("Apply to Qubit 1") {
+      ap.apply_gate(OpType::Reset, {1});
+      ap.normal_form();
+      // Correct form verified by hand
+      MatrixXb A = MatrixXb::Zero(4, 4);
+      VectorXb B = VectorXb::Zero(4);
+      MatrixXb C = MatrixXb::Zero(4, 4);
+      MatrixXb E = MatrixXb::Zero(4, 4);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+      A(0, 1) = true;
+      C(0, 0) = C(0, 3) = true;
+      C(1, 2) = true;
+      E(0, 3) = E(3, 0) = true;
+      E(2, 3) = E(3, 2) = true;
+      P(3) = 2;
+      APState correct(A, B, C, E, P, {0});
+      // Check equality up to global phase
+      ap.phase = correct.phase;
+      CHECK(ap == correct);
+    }
+    WHEN("Apply to Qubit 2") {
+      ap.apply_gate(OpType::Reset, {2});
+      ap.normal_form();
+      // Correct form verified by hand
+      MatrixXb A = MatrixXb::Zero(4, 4);
+      VectorXb B = VectorXb::Zero(4);
+      MatrixXb C = MatrixXb::Zero(4, 4);
+      MatrixXb E = MatrixXb::Zero(4, 4);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+      A(0, 0) = A(0, 1) = A(0, 3) = true;
+      A(1, 2) = true;
+      B(0) = true;
+      C(0, 1) = true;
+      C(1, 3) = true;
+      APState correct(A, B, C, E, P, {0});
+      // Check equality up to global phase
+      ap.phase = correct.phase;
+      CHECK(ap == correct);
+    }
+    WHEN("Apply to Qubit 3") {
+      ap.apply_gate(OpType::Reset, {3});
+      ap.normal_form();
+      // Correct form verified by hand
+      MatrixXb A = MatrixXb::Zero(4, 4);
+      VectorXb B = VectorXb::Zero(4);
+      MatrixXb C = MatrixXb::Zero(4, 4);
+      MatrixXb E = MatrixXb::Zero(4, 4);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+      A(0, 3) = true;
+      C(0, 0) = C(0, 2) = true;
+      C(1, 1) = C(1, 2) = true;
+      P(2) = 1;
+      APState correct(A, B, C, E, P, {0});
+      // Check equality up to global phase
+      ap.phase = correct.phase;
+      CHECK(ap == correct);
+    }
+  }
+}
+
+SCENARIO("Gate encodings") {
+  std::list<std::pair<OpType, std::vector<unsigned>>> test_gates = {
+      {OpType::Z, {0}},       {OpType::X, {0}},
+      {OpType::Y, {0}},       {OpType::S, {0}},
+      {OpType::Sdg, {0}},     {OpType::V, {0}},
+      {OpType::Vdg, {0}},     {OpType::SX, {0}},
+      {OpType::SXdg, {0}},    {OpType::H, {0}},
+      {OpType::CX, {0, 1}},   {OpType::CY, {0, 1}},
+      {OpType::CZ, {0, 1}},   {OpType::ZZMax, {0, 1}},
+      {OpType::ECR, {0, 1}},  {OpType::ISWAPMax, {0, 1}},
+      {OpType::SWAP, {0, 1}}, {OpType::BRIDGE, {0, 1, 2}},
+      {OpType::noop, {0}},
+  };
+  GIVEN("Check Z actions") {
+    for (const auto& com : test_gates) {
+      MatrixXb A = MatrixXb::Identity(3, 3);
+      VectorXb B = VectorXb::Zero(3);
+      MatrixXb E = MatrixXb::Zero(3, 3);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+      test_apply_gate(A, B, E, P, com.first, com.second);
+    }
+  }
+  GIVEN("Check X actions") {
+    for (const auto& com : test_gates) {
+      MatrixXb A = MatrixXb::Zero(3, 3);
+      VectorXb B = VectorXb::Zero(3);
+      MatrixXb E = MatrixXb::Zero(3, 3);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+      test_apply_gate(A, B, E, P, com.first, com.second);
+    }
+  }
+  GIVEN("Check Z actions on inputs of ChoiAPState") {
+    for (const auto& com : test_gates) {
+      MatrixXb A = MatrixXb::Identity(3, 3);
+      VectorXb B = VectorXb::Zero(3);
+      MatrixXb C = MatrixXb::Zero(3, 3);
+      MatrixXb E = MatrixXb::Zero(3, 3);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+      qubit_vector_t qbs;
+      for (unsigned q : com.second) qbs.push_back(Qubit(q));
+      test_apply_gate_dm_input(A, B, C, E, P, com.first, qbs);
+    }
+  }
+  GIVEN("Check X actions on inputs of ChoiAPState") {
+    for (const auto& com : test_gates) {
+      MatrixXb A = MatrixXb::Zero(3, 3);
+      VectorXb B = VectorXb::Zero(3);
+      MatrixXb C = MatrixXb::Zero(3, 3);
+      MatrixXb E = MatrixXb::Zero(3, 3);
+      Eigen::VectorXi P = Eigen::VectorXi::Zero(3);
+      qubit_vector_t qbs;
+      for (unsigned q : com.second) qbs.push_back(Qubit(q));
+      test_apply_gate_dm_input(A, B, C, E, P, com.first, qbs);
+    }
+  }
+}
+
+SCENARIO("Loading from a statevector") {
+  MatrixXb A = MatrixXb::Zero(4, 4);
+  VectorXb B = VectorXb::Zero(4);
+  MatrixXb C = MatrixXb::Zero(4, 4);
+  MatrixXb E = MatrixXb::Zero(4, 4);
+  Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+  A(0, 0) = A(0, 2) = A(0, 3) = true;
+  A(1, 1) = A(1, 2) = true;
+  B(0) = true;
+  E(2, 3) = E(3, 2) = true;
+  P(2) = 1;
+  P(3) = 2;
+  APState ap(A, B, C, E, P, 0);
+  auto sv = ap.to_statevector();
+  APState reconstructed(sv);
+  auto sv2 = reconstructed.to_statevector();
+  CHECK(tket_sim::compare_statevectors_or_unitaries(
+      sv, sv2, tket_sim::MatrixEquivalence::EQUAL));
+}
+
+SCENARIO("Loading from a density matrix") {
+  MatrixXb A = MatrixXb::Zero(4, 4);
+  VectorXb B = VectorXb::Zero(4);
+  MatrixXb C = MatrixXb::Zero(4, 4);
+  MatrixXb E = MatrixXb::Zero(4, 4);
+  Eigen::VectorXi P = Eigen::VectorXi::Zero(4);
+  A(0, 0) = A(0, 2) = A(0, 3) = true;
+  C(0, 1) = C(0, 2) = true;
+  C(1, 0) = true;
+  B(0) = true;
+  E(2, 3) = E(3, 2) = true;
+  P(0) = 3;
+  P(1) = 2;
+  P(2) = 1;
+  P(3) = 2;
+  APState ap(A, B, C, E, P, 0);
+  auto dm = ap.to_density_matrix();
+  APState reconstructed(dm);
+  auto dm2 = reconstructed.to_density_matrix();
+  CHECK(dm2.isApprox(dm, EPS));
+  // This state is mixed, so only check equality of normal forms up to phase
+  ap.normal_form();
+  reconstructed.normal_form();
+  reconstructed.phase = ap.phase;
+  CHECK(ap == reconstructed);
+  THEN("Test serialisation") {
+    nlohmann::json j_ap = ap;
+    APState ap2(0);
+    j_ap.get_to(ap2);
+    REQUIRE(ap == ap2);
+  }
+}
+
+SCENARIO("Converting from/to a circuit") {
+  GIVEN("A pure circuit in the standard AP form") {
+    Circuit circ(4);
+    circ.qubit_create_all();
+    circ.add_op<unsigned>(OpType::H, {2});
+    circ.add_op<unsigned>(OpType::H, {3});
+    circ.add_op<unsigned>(OpType::CX, {2, 0});
+    circ.add_op<unsigned>(OpType::CX, {2, 1});
+    circ.add_op<unsigned>(OpType::CX, {3, 1});
+    circ.add_op<unsigned>(OpType::S, {2});
+    circ.add_op<unsigned>(OpType::Z, {3});
+    APState ap = circuit_to_apstate(circ);
+    auto sv_circ = tket_sim::get_statevector(circ);
+    auto sv_ap = ap.to_statevector();
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL));
+    ap.normal_form();
+    sv_ap = ap.to_statevector();
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL));
+    Circuit reconstructed = apstate_to_circuit(ap);
+    CHECK(circ == reconstructed);
+  }
+  GIVEN("A generic pure circuit") {
+    Circuit circ(4);
+    circ.qubit_create_all();
+    circ.add_op<unsigned>(OpType::V, {0});
+    circ.add_op<unsigned>(OpType::CX, {0, 1});
+    circ.add_op<unsigned>(OpType::CY, {1, 3});
+    circ.add_op<unsigned>(OpType::H, {3});
+    circ.add_op<unsigned>(OpType::ZZMax, {2, 3});
+    APState ap = circuit_to_apstate(circ);
+    auto sv_circ = tket_sim::get_statevector(circ);
+    auto sv_ap = ap.to_statevector();
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL));
+    Circuit reconstructed = apstate_to_circuit(ap);
+    auto sv_rec = tket_sim::get_statevector(reconstructed);
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_circ, sv_rec, tket_sim::MatrixEquivalence::EQUAL));
+  }
+  GIVEN("Initialisations, collapses, discards and post-selections") {
+    Circuit circ(5);
+    circ.qubit_create(Qubit(1));
+    circ.qubit_create(Qubit(2));
+    circ.add_op<unsigned>(OpType::H, {4});
+    circ.add_op<unsigned>(OpType::Collapse, {4});
+    circ.add_op<unsigned>(OpType::CX, {4, 1});
+    circ.add_op<unsigned>(OpType::CX, {4, 2});
+    circ.add_op<unsigned>(OpType::CX, {4, 3});
+    circ.add_op<unsigned>(OpType::H, {4});
+    circ.add_op<unsigned>(OpType::H, {1});
+    circ.add_op<unsigned>(OpType::V, {2});
+    circ.add_op<unsigned>(OpType::CX, {1, 2});
+    circ.add_op<unsigned>(OpType::CX, {1, 0});
+    circ.qubit_discard(Qubit(0));
+    ChoiAPState ap = circuit_to_choi_apstate(circ);
+    ap.post_select(Qubit(3), ChoiAPState::TableauSegment::Output);
+    ap.canonical_column_order();
+    ap.normal_form();
+    // Define correct form from a hand calculation
+    MatrixXb A = MatrixXb::Zero(6, 6);
+    VectorXb B = VectorXb::Zero(6);
+    MatrixXb C = MatrixXb::Zero(6, 6);
+    C(0, 0) = C(0, 3) = true;
+    C(1, 1) = true;
+    MatrixXb E = MatrixXb::Zero(6, 6);
+    E(1, 2) = E(2, 1) = true;
+    E(1, 4) = E(4, 1) = true;
+    E(1, 5) = E(5, 1) = true;
+    E(3, 4) = E(4, 3) = true;
+    Eigen::VectorXi P = Eigen::VectorXi::Zero(6);
+    P(3) = P(4) = 3;
+    // Ignore phase by setting them to match
+    APState correct(A, B, C, E, P, ap.ap_.phase);
+    CHECK(ap.ap_ == correct);
+    std::pair<Circuit, qubit_map_t> res_uni =
+        choi_apstate_to_unitary_extension_circuit(ap, {Qubit(1)}, {Qubit(0)});
+    // Rebuild state by initialising, post-selecting, etc.
+    ChoiAPState res_ap = circuit_to_choi_apstate(res_uni.first);
+    qubit_map_t perm;
+    for (const std::pair<const Qubit, Qubit>& p : res_uni.second)
+      perm.insert({p.second, p.first});
+    res_ap.rename_qubits(perm, ChoiAPState::TableauSegment::Output);
+    // Post-select/initialise
+    res_ap.post_select(Qubit(1), ChoiAPState::TableauSegment::Input);
+    res_ap.post_select(Qubit(0), ChoiAPState::TableauSegment::Output);
+    // Collapsing q[4] in X basis as per circ
+    res_ap.apply_gate(
+        OpType::H, {Qubit(4)}, ChoiAPState::TableauSegment::Output);
+    res_ap.collapse_qubit(Qubit(4), ChoiAPState::TableauSegment::Output);
+    res_ap.apply_gate(
+        OpType::H, {Qubit(4)}, ChoiAPState::TableauSegment::Output);
+    // Discarding q[0] also removes Z row for q[0], so recreate this by
+    // XCollapse at input
+    res_ap.apply_gate(
+        OpType::H, {Qubit(0)}, ChoiAPState::TableauSegment::Input);
+    res_ap.collapse_qubit(Qubit(0), ChoiAPState::TableauSegment::Input);
+    res_ap.apply_gate(
+        OpType::H, {Qubit(0)}, ChoiAPState::TableauSegment::Input);
+    res_ap.canonical_column_order();
+    res_ap.normal_form();
+    // Mixed state, so only guaranteed up to phase
+    res_ap.ap_.phase = ap.ap_.phase;
+    CHECK(res_ap == ap);
+    THEN("Serialize and deserialize") {
+      nlohmann::json j_ap = ap;
+      ChoiAPState ap2({});
+      j_ap.get_to(ap2);
+      CHECK(ap == ap2);
+    }
+    THEN("Check conversion to/from a tableau") {
+      ChoiMixTableau tab = choi_apstate_to_cm_tableau(ap);
+      ChoiMixTableau tab2 = circuit_to_cm_tableau(circ);
+      tab2.post_select(Qubit(3), ChoiMixTableau::TableauSegment::Output);
+      tab.canonical_column_order();
+      tab.gaussian_form();
+      tab2.canonical_column_order();
+      tab2.gaussian_form();
+      REQUIRE(tab == tab2);
+      ChoiAPState ap2 = cm_tableau_to_choi_apstate(tab);
+      ap2.canonical_column_order();
+      ap2.normal_form();
+      // Converting to a tableau drops phase, so ignore this in equivalence
+      // check
+      ap2.ap_.phase = ap.ap_.phase;
+      REQUIRE(ap == ap2);
+    }
+  }
+}
+
+SCENARIO("Converting from/to a tableau") {
+  GIVEN("Check pure state up to global phase using circuit") {
+    Circuit circ(8);
+    circ.qubit_create_all();
+    circ.add_op<unsigned>(OpType::X, {1});
+    circ.add_op<unsigned>(OpType::X, {5});
+    circ.add_op<unsigned>(OpType::H, {2});
+    circ.add_op<unsigned>(OpType::H, {4});
+    circ.add_op<unsigned>(OpType::H, {6});
+    circ.add_op<unsigned>(OpType::H, {7});
+    circ.add_op<unsigned>(OpType::CX, {2, 1});
+    circ.add_op<unsigned>(OpType::CX, {4, 0});
+    circ.add_op<unsigned>(OpType::CX, {4, 3});
+    circ.add_op<unsigned>(OpType::CX, {6, 0});
+    circ.add_op<unsigned>(OpType::CX, {6, 1});
+    circ.add_op<unsigned>(OpType::CX, {7, 5});
+    circ.add_op<unsigned>(OpType::CZ, {2, 6});
+    circ.add_op<unsigned>(OpType::CZ, {4, 6});
+    circ.add_op<unsigned>(OpType::CZ, {4, 7});
+    circ.add_op<unsigned>(OpType::CZ, {6, 7});
+    circ.add_op<unsigned>(OpType::S, {2});
+    circ.add_op<unsigned>(OpType::Sdg, {4});
+    circ.add_op<unsigned>(OpType::Z, {7});
+    ChoiMixTableau cmt = circuit_to_cm_tableau(circ);
+    ChoiAPState ap = cm_tableau_to_choi_apstate(cmt);
+    auto sv_circ = tket_sim::get_statevector(circ);
+    auto sv_ap = ap.ap_.to_statevector();
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_circ, sv_ap, tket_sim::MatrixEquivalence::EQUAL));
+    ChoiMixTableau cmt2 = choi_apstate_to_cm_tableau(ap);
+    auto [circ2, perm] = cm_tableau_to_exact_circuit(cmt2);
+    qubit_map_t inv;
+    for (const std::pair<const Qubit, Qubit>& qp : perm)
+      inv.insert({qp.second, qp.first});
+    circ2.permute_boundary_output(inv);
+    auto sv_circ2 = tket_sim::get_statevector(circ2);
+    CHECK(tket_sim::compare_statevectors_or_unitaries(
+        sv_circ, sv_circ2,
+        tket_sim::MatrixEquivalence::EQUAL_UP_TO_GLOBAL_PHASE));
+  }
+}
+
+}  // namespace test_APState
+}  // namespace tket
\ No newline at end of file