Repository: cda-tum/mqt-qmap (language: python)

Issue:
Clifford synthesis
Adding the functionality to synthesize clifford circuits from tableaus and unoptimized input circuits, with some structural changes.

First, Adding LogicBlocks as a dependency and the logging functionality from there.
Introduction of Cliffordresults class to contain all Tableaus and Resulting Circuits.
Introduction of Tableaus as data structure to represent Clifford circuits.

ToDos:

- [x] Check integration of Logicblocks for sneaky errors 
- [x] Check functionality of Cliffordsynthesiszer
- [x] Write Tests for new Classes
- [x] Add ability to load a tableau from a tableau string as in qiskit (https://qiskit.org/documentation/stubs/qiskit.quantum_info.Clifford.html)
- [x] Add Python bindings


Discussion / hints:

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index fc94410d6..fb9346ba3 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -117,7 +117,7 @@ repos:
     rev: "v2.2.2"
     hooks:
       - id: codespell
-        args: ["-L", "wille,linz", "--skip", "*.ipynb"]
+        args: ["-L", "wille,linz,applys", "--skip", "*.ipynb"]
 
   # Clang-format the C++ part of the code base automatically
   - repo: https://github.com/pre-commit/mirrors-clang-format
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 71c1ac5c0..5a76cf041 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -65,6 +65,9 @@ if(NOT Z3_FOUND)
   message(WARNING "Did not find Z3. Exact library and other depending target will not be available")
 endif()
 
+set(THREADS_PREFER_PTHREAD_FLAG ON)
+find_package(Threads REQUIRED)
+
 # add main library code
 add_subdirectory(src)
 
diff --git a/include/cliffordsynthesis/Tableau.hpp b/include/cliffordsynthesis/Tableau.hpp
new file mode 100644
index 000000000..97e01a3be
--- /dev/null
+++ b/include/cliffordsynthesis/Tableau.hpp
@@ -0,0 +1,132 @@
+//
+// This file is part of the MQT QMAP library released under the MIT license.
+// See README.md or go to https://github.com/cda-tum/qmap for more information.
+//
+
+#pragma once
+
+#include "QuantumComputation.hpp"
+#include "utils/logging.hpp"
+
+#include <cstdint>
+#include <fstream>
+#include <limits>
+#include <ostream>
+#include <utility>
+#include <vector>
+
+namespace cs {
+class Tableau {
+  using EntryType   = std::uint8_t;
+  using RowType     = std::vector<EntryType>;
+  using TableauType = std::vector<RowType>;
+  std::size_t nQubits{};
+  TableauType tableau;
+
+public:
+  Tableau() = default;
+  explicit Tableau(const qc::QuantumComputation& qc, std::size_t begin = 0,
+                   std::size_t end = std::numeric_limits<std::size_t>::max());
+  explicit Tableau(const std::size_t nQubits) : nQubits(nQubits) {
+    createDiagonalTableau(nQubits);
+  }
+  explicit Tableau(const std::string& description) {
+    fromString(description);
+    nQubits = tableau.size();
+  }
+
+  [[nodiscard]] RowType operator[](const std::size_t index) {
+    return tableau[index];
+  }
+  [[nodiscard]] RowType operator[](const std::size_t index) const {
+    return tableau[index];
+  }
+
+  [[nodiscard]] RowType at(const std::size_t index) {
+    return tableau.at(index);
+  }
+
+  [[nodiscard]] auto getQubitCount() const { return nQubits; }
+
+  void dump(const std::string& filename) const;
+
+  void dump(std::ostream& of) const;
+
+  void import(const std::string& filename);
+  void import(std::istream& is);
+
+  template <std::size_t N>
+  void populateTableauFrom(const std::bitset<N> bv, const std::size_t nQ,
+                           const std::size_t column) {
+    assert(column <= 2 * nQ);
+    assert(nQ <= nQubits);
+    assert(nQ <= N);
+    for (std::size_t i = 0U; i < nQ; ++i) {
+      if (bv[i]) {
+        tableau[i][column] = 1;
+      } else {
+        tableau[i][column] = 0;
+      }
+    }
+  }
+  void populateTableauFrom(const std::uint64_t bv, const std::size_t nQ,
+                           const std::size_t column) {
+    populateTableauFrom(std::bitset<64>(bv), nQ, column);
+  }
+
+  void applyGate(const qc::Operation* gate);
+  void applyH(std::size_t target);
+  void applyS(std::size_t target);
+  void applySdag(std::size_t target);
+  void applySx(std::size_t target);
+  void applySxdag(std::size_t target);
+  void applyX(std::size_t target);
+  void applyY(std::size_t target);
+  void applyZ(std::size_t target);
+  void applyCX(std::size_t control, std::size_t target);
+  void applyCY(std::size_t control, std::size_t target);
+  void applyCZ(std::size_t control, std::size_t target);
+  void applySwap(std::size_t q1, std::size_t q2);
+
+  friend bool operator==(const Tableau& lhs, const Tableau& rhs) {
+    return lhs.tableau == rhs.tableau;
+  }
+  friend bool operator!=(const Tableau& lhs, const Tableau& rhs) {
+    return !(lhs == rhs);
+  }
+
+  void createDiagonalTableau(std::size_t nQ);
+
+  friend std::ostream& operator<<(std::ostream& os, const Tableau& dt) {
+    os << dt.toString();
+    return os;
+  }
+  friend std::istream& operator>>(std::istream& is, Tableau& dt) {
+    if (is.good()) {
+      std::stringstream ss;
+      ss << is.rdbuf();
+      dt.fromString(ss.str());
+    }
+    return is;
+  }
+
+  [[nodiscard]] std::string toString() const;
+  void                      fromString(const std::string& str);
+
+  template <std::size_t N>
+  [[nodiscard]] std::bitset<N> getBVFrom(const std::size_t column) const {
+    assert(column <= 2 * nQubits);
+    assert(nQubits <= N);
+    std::bitset<N> bv;
+    for (std::size_t i = 0U; i < nQubits; ++i) {
+      if (tableau[i][column] == 1U) {
+        bv[i] = 1;
+      }
+    }
+    return bv;
+  }
+  [[nodiscard]] std::uint64_t getBVFrom(const std::size_t column) const {
+    return getBVFrom<64>(column).to_ullong();
+  }
+};
+} // namespace cs
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 97f92be38..e145746dc 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,23 +4,8 @@ add_subdirectory("${PROJECT_SOURCE_DIR}/extern/qfr" "extern/qfr")
 add_subdirectory("${PROJECT_SOURCE_DIR}/extern/LogicBlocks" "extern/LogicBlocks")
 
 # macro to add core libraries
-macro(ADD_QMAP_LIBRARY libname srcfile)
-  set(lib ${PROJECT_NAME}_${libname}_lib)
-
-  # main project library
-  add_library(
-    ${lib}
-    ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/Architecture.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/Mapper.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/${libname}/${srcfile}.cpp
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/utils.hpp
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/Architecture.hpp
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/Mapper.hpp
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/MappingResults.hpp
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/configuration
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/${libname}/${srcfile}.hpp
-    ${${PROJECT_NAME}_SOURCE_DIR}/include/heuristic/unique_priority_queue.hpp)
+macro(ADD_INTERNAL_LIBRARY libname srcfile)
+  set(lib ${libname})
 
   # set include directories
   target_include_directories(${lib}
@@ -44,7 +29,7 @@ macro(ADD_QMAP_LIBRARY libname srcfile)
   if(MSVC)
     target_compile_options(${lib} PUBLIC /utf-8)
   else()
-    target_compile_options(${lib} PUBLIC -Wall -Wextra $<$<CONFIG:DEBUG>:-Og>)
+    target_compile_options(${lib} PUBLIC -Wall -Wextra $<$<CONFIG:DEBUG>:-g3>)
     if(BINDINGS)
       # adjust visibility settings for building Python bindings
       target_compile_options(${lib} PUBLIC -fvisibility=hidden)
@@ -76,10 +61,58 @@ macro(ADD_QMAP_LIBRARY libname srcfile)
   # add MQT alias
   add_library(MQT::${lib} ALIAS ${lib})
 endmacro()
+# macro to add core libraries
+macro(ADD_QMAP_LIBRARY libname srcfile)
+  set(lib ${PROJECT_NAME}_${libname}_lib)
 
+  # main project library
+  add_library(
+    ${lib}
+    ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Architecture.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Mapper.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/${libname}/${srcfile}.cpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/utils.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/Architecture.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/Mapper.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/MappingResults.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/configuration
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/${libname}/${srcfile}.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/heuristic/unique_priority_queue.hpp)
+
+  add_internal_library(${lib} ${srcfile})
+endmacro()
+# macro to add core libraries
+macro(ADD_SYNTHESIS_LIBRARY libname srcfile)
+  set(lib ${PROJECT_NAME}_${libname}_lib)
+
+  # main project library
+  add_library(
+    ${lib}
+    ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/Architecture.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/cliffordsynthesis/Tableau.cpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/utils.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/Architecture.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/cliffordsynthesis/Tableau.hpp
+    ${${PROJECT_NAME}_SOURCE_DIR}/include/configuration/AvailableArchitecture.hpp)
+
+  add_internal_library(${lib} ${srcfile})
+endmacro()
+
+# exact mapper project library
 add_qmap_library(exact ExactMapper)
+
 # add LogicBlocks library
 target_link_libraries(${PROJECT_NAME}_exact_lib PUBLIC logicblocks::Logic logicblocks::Encodings)
 
+add_synthesis_library(cliffordsynthesis CliffordSynthesizer)
+
+# add LogicBlocks library
+target_link_libraries(${PROJECT_NAME}_cliffordsynthesis_lib PUBLIC logicblocks::Logic
+                                                                   logicblocks::Encodings)
+# add Pthreads
+target_link_libraries(${PROJECT_NAME}_cliffordsynthesis_lib PRIVATE Threads::Threads)
+
 # heuristic mapper project library
 add_qmap_library(heuristic HeuristicMapper)
diff --git a/src/cliffordsynthesis/Tableau.cpp b/src/cliffordsynthesis/Tableau.cpp
new file mode 100644
index 000000000..d92821ff2
--- /dev/null
+++ b/src/cliffordsynthesis/Tableau.cpp
@@ -0,0 +1,325 @@
+//
+// This file is part of the MQT QMAP library released under the MIT license.
+// See README.md or go to https://github.com/cda-tum/qmap for more information.
+//
+
+#include "cliffordsynthesis/Tableau.hpp"
+
+#include "utils.hpp"
+
+namespace cs {
+void Tableau::dump(const std::string& filename) const {
+  auto of = std::ofstream(filename);
+  if (!of.good()) {
+    util::fatal("Error opening file " + filename);
+  }
+  dump(of);
+}
+
+void Tableau::dump(std::ostream& of) const { of << *this; }
+
+void Tableau::import(const std::string& filename) {
+  auto is = std::ifstream(filename);
+  if (!is.good()) {
+    util::fatal("Error opening file " + filename);
+  }
+  import(is);
+}
+
+void Tableau::import(std::istream& is) {
+  tableau.clear();
+
+  std::string              line;
+  std::vector<std::string> data{};
+  char                     delimiter = '|';
+
+  while (std::getline(is, line)) {
+    if (line.find('|', 0) == std::string::npos) {
+      delimiter = ';';
+    }
+    tableau.emplace_back();
+    ::parse_line(line, delimiter, {'\"'}, {'\\', '\r', '\n', '\t'}, data);
+    for (const auto& datum : data) {
+      if (datum.empty()) {
+        continue;
+      }
+      tableau.back().emplace_back(static_cast<EntryType>(std::stoul(datum)));
+    }
+  }
+}
+
+void Tableau::applyGate(const qc::Operation* const gate) {
+  if (gate->getNcontrols() > 1U) {
+    util::fatal("Tableau::applyGate: Only operations with up to one control "
+                "are supported.");
+  }
+
+  const auto target = static_cast<std::size_t>(gate->getTargets().at(0U));
+
+  // non-controlled gates
+  if (!gate->isControlled()) {
+    switch (gate->getType()) {
+    case qc::OpType::H:
+      applyH(target);
+      break;
+    case qc::OpType::S:
+      applyS(target);
+      break;
+    case qc::OpType::Sdag:
+      applySdag(target);
+      break;
+    case qc::OpType::SX:
+      applySx(target);
+      break;
+    case qc::OpType::SXdag:
+      applySxdag(target);
+      break;
+    case qc::OpType::X:
+      applyX(target);
+      break;
+    case qc::OpType::Y:
+      applyY(target);
+      break;
+    case qc::OpType::Z:
+      applyZ(target);
+      break;
+    case qc::OpType::SWAP: {
+      const auto target2 = static_cast<std::size_t>(gate->getTargets().at(1U));
+      applySwap(target, target2);
+      break;
+    }
+    default:
+      // unsupported non-controlled gate type
+      util::fatal("Tableau::applyGate: Unsupported non-controlled gate type " +
+                  qc::toString(gate->getType()));
+    }
+  } else {
+    const auto control =
+        static_cast<std::size_t>((*gate->getControls().begin()).qubit);
+    switch (gate->getType()) {
+    case qc::OpType::X:
+      applyCX(control, target);
+      break;
+    case qc::OpType::Y:
+      applyCY(control, target);
+      break;
+    case qc::OpType::Z:
+      applyCZ(control, target);
+      break;
+    default:
+      // unsupported controlled gate type
+      util::fatal("Tableau::applyGate: Unsupported controlled gate type " +
+                  qc::toString(gate->getType()));
+    }
+  }
+}
+
+void Tableau::createDiagonalTableau(const std::size_t nQ) {
+  nQubits = nQ;
+  tableau.clear();
+  tableau.resize(nQubits);
+  for (std::size_t i = 0U; i < nQubits; ++i) {
+    tableau[i].resize((2U * nQubits) + 1U);
+    for (std::size_t j = nQubits; j < (2U * nQubits); ++j) {
+      tableau[i][j] = 0;
+      if (i == (j - nQubits)) {
+        tableau[i][j] = 1;
+      }
+    }
+  }
+}
+
+std::string Tableau::toString() const {
+  std::stringstream ss;
+  for (const auto& row : tableau) {
+    if (row.size() != tableau.back().size()) {
+      util::fatal("Tableau is not rectangular.");
+      return "";
+    }
+    for (const auto& s : row) {
+      ss << std::to_string(s) << ';';
+    }
+    ss << "\n";
+  }
+  return ss.str();
+}
+
+void Tableau::fromString(const std::string& str) {
+  std::stringstream ss(str);
+  std::string       line;
+  std::getline(ss, line);
+  if (line.empty()) {
+    return;
+  }
+  const auto  rStabilizer = std::regex("([\\+-])([IYZX]+)");
+  std::smatch m;
+  if (std::regex_search(line, rStabilizer)) {
+    // string is a list of stabilizers
+    auto iter = line.cbegin();
+    while (std::regex_search(iter, line.cend(), m, rStabilizer)) {
+      std::string s = m.str(0U);
+      RowType     row;
+
+      for (const auto c : s) {
+        if (c == 'I' || c == 'Z') {
+          row.push_back(0);
+        } else if (c == 'X' || c == 'Y') {
+          row.push_back(1);
+        }
+      }
+      for (const auto c : s) {
+        if (c == 'I' || c == 'X') {
+          row.push_back(0);
+        } else if (c == 'Y' || c == 'Z') {
+          row.push_back(1);
+        }
+      }
+      if (s[0U] == '-') {
+        row.push_back(1);
+      } else {
+        row.push_back(0);
+      }
+      tableau.push_back(row);
+      iter = m[0].second;
+    }
+  } else {
+    // assume string is a semicolon separated binary matrix
+    ss = std::stringstream(str);
+    import(ss);
+  }
+}
+
+void Tableau::applyH(const std::size_t target) {
+  assert(target < nQubits);
+  for (auto i = 0U; i < nQubits; ++i) {
+    tableau[i][2U * nQubits] ^=
+        (tableau[i][target] & tableau[i][target + nQubits]);
+    std::swap(tableau[i][target], tableau[i][target + nQubits]);
+  }
+}
+
+void Tableau::applyS(const std::size_t target) {
+  assert(target < nQubits);
+  for (auto i = 0U; i < nQubits; ++i) {
+    tableau[i][2U * nQubits] ^=
+        tableau[i][target] & tableau[i][target + nQubits];
+    tableau[i][target + nQubits] ^= tableau[i][target];
+  }
+}
+
+// Sdag = S * S * S
+void Tableau::applySdag(const std::size_t target) {
+  assert(target < nQubits);
+  applyS(target);
+  applyS(target);
+  applyS(target);
+}
+
+// Sx = Sdag * H * Sdag
+void Tableau::applySx(const std::size_t target) {
+  assert(target < nQubits);
+  applySdag(target);
+  applyH(target);
+  applySdag(target);
+}
+
+// Sxdag = S * H * S
+void Tableau::applySxdag(const std::size_t target) {
+  assert(target < nQubits);
+  applyS(target);
+  applyH(target);
+  applyS(target);
+}
+
+// X = H * Z * H
+void Tableau::applyX(const std::size_t target) {
+  assert(target < nQubits);
+  applyH(target);
+  applyZ(target);
+  applyH(target);
+}
+
+// Y = X * Z
+void Tableau::applyY(const std::size_t target) {
+  assert(target < nQubits);
+  applyX(target);
+  applyZ(target);
+}
+
+// Z = S * S
+void Tableau::applyZ(const std::size_t target) {
+  assert(target < nQubits);
+  applyS(target);
+  applyS(target);
+}
+
+void Tableau::applyCX(const std::size_t control, const std::size_t target) {
+  assert(control < nQubits);
+  assert(target < nQubits);
+  assert(control != target);
+  for (auto i = 0U; i < nQubits; ++i) {
+    const auto xa = tableau[i][target];
+    const auto za = tableau[i][target + nQubits];
+    const auto xb = tableau[i][control];
+    const auto zb = tableau[i][control + nQubits];
+    tableau[i][2 * nQubits] ^= (xa & zb) & ((xb ^ za) ^ 1);
+    tableau[i][control + nQubits] = za ^ zb;
+    tableau[i][target]            = xb ^ xa;
+  }
+}
+
+void Tableau::applyCY(const std::size_t control, const std::size_t target) {
+  assert(control < nQubits);
+  assert(target < nQubits);
+  assert(control != target);
+  applySdag(target);
+  applyCX(control, target);
+  applyS(target);
+}
+
+void Tableau::applyCZ(const std::size_t control, const std::size_t target) {
+  assert(control < nQubits);
+  assert(target < nQubits);
+  assert(control != target);
+  applyH(target);
+  applyCX(control, target);
+  applyH(target);
+}
+
+void Tableau::applySwap(const std::size_t q1, const std::size_t q2) {
+  assert(q1 < nQubits);
+  assert(q2 < nQubits);
+  assert(q1 != q2);
+  applyCX(q1, q2);
+  applyCX(q2, q1);
+  applyCX(q1, q2);
+}
+
+Tableau::Tableau(const qc::QuantumComputation& qc, const std::size_t begin,
+                 const std::size_t end)
+    : Tableau(qc.getNqubits()) {
+  std::size_t currentG = 0;
+  for (const auto& gate : qc) {
+    if (gate->getType() == qc::OpType::Compound) {
+      const auto* const compOp =
+          dynamic_cast<const qc::CompoundOperation* const>(gate.get());
+      auto cit = compOp->begin();
+      while (cit != compOp->end()) {
+        if (currentG >= begin && (currentG < end)) {
+          applyGate((*cit).get());
+        }
+        ++cit;
+        ++currentG;
+      }
+    } else {
+      if (currentG >= begin && (currentG < end)) {
+        applyGate(gate.get());
+      }
+      ++currentG;
+    }
+    if (currentG >= end) {
+      break;
+    }
+  }
+}
+} // namespace cs


Test patch:
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 71841688f..ff7dc8143 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -74,3 +74,8 @@ if(TARGET ${PROJECT_NAME}_exact_lib)
   package_add_test(${PROJECT_NAME}_exact_test exact ${CMAKE_CURRENT_SOURCE_DIR}/test_exact.cpp
                    ${CMAKE_CURRENT_SOURCE_DIR}/test_encodings.cpp)
 endif()
+
+if(TARGET ${PROJECT_NAME}_cliffordsynthesis_lib)
+  package_add_test(${PROJECT_NAME}_cliffordsynthesis_test cliffordsynthesis
+                   ${CMAKE_CURRENT_SOURCE_DIR}/test_tableau.cpp)
+endif()
diff --git a/test/test_tableau.cpp b/test/test_tableau.cpp
new file mode 100644
index 000000000..4fe22d5ab
--- /dev/null
+++ b/test/test_tableau.cpp
@@ -0,0 +1,326 @@
+//
+// This file is part of the MQT QMAP library released under the MIT license.
+// See README.md or go to https://github.com/cda-tum/qmap for more information.
+//
+
+#include "QuantumComputation.hpp"
+#include "cliffordsynthesis/Tableau.hpp"
+
+#include "gtest/gtest.h"
+
+namespace cs {
+
+class TestTableau : public ::testing::Test {
+protected:
+  // 0 0 | 1 0 | 0
+  // 0 0 | 0 1 | 0
+  Tableau tableau{2};
+};
+
+TEST_F(TestTableau, InitialTableau) {
+  EXPECT_EQ(tableau.getQubitCount(), 2);
+  // X part is zero
+  EXPECT_EQ(tableau[0][0], 0);
+  EXPECT_EQ(tableau[0][1], 0);
+  EXPECT_EQ(tableau[1][0], 0);
+  EXPECT_EQ(tableau[1][1], 0);
+  // Z part is identity
+  EXPECT_EQ(tableau[0][2], 1);
+  EXPECT_EQ(tableau[0][3], 0);
+  EXPECT_EQ(tableau[1][2], 0);
+  EXPECT_EQ(tableau[1][3], 1);
+  // R part is zero
+  EXPECT_EQ(tableau[0][4], 0);
+  EXPECT_EQ(tableau[1][4], 0);
+
+  std::stringstream ss;
+  ss << tableau;
+
+  const std::string representation = "0;0;1;0;0;\n"
+                                     "0;0;0;1;0;\n";
+  EXPECT_EQ(ss.str(), representation);
+
+  const auto tFromStr = Tableau(representation);
+  EXPECT_EQ(tableau, tFromStr);
+
+  const std::string stabilizers      = "[+ZI, +IZ]";
+  const auto        tFromStabilizers = Tableau(stabilizers);
+  EXPECT_EQ(tableau, tFromStabilizers);
+}
+
+TEST_F(TestTableau, H) {
+  // H on |0> is |+>, which is stabilized by X
+  tableau.applyH(0);
+  std::string expected = "1;0;0;0;0;\n"
+                         "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+XI, +IZ]"));
+
+  tableau.applyH(1);
+  expected = "1;0;0;0;0;\n"
+             "0;1;0;0;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+XI, +IX]"));
+
+  tableau.applyH(1);
+  expected = "1;0;0;0;0;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+XI, +IZ]"));
+
+  tableau.applyH(0);
+  expected = "0;0;1;0;0;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +IZ]"));
+}
+
+TEST_F(TestTableau, X) {
+  // X on |0> is |1>, which is stabilized by -Z
+  tableau.applyX(0);
+  std::string expected = "0;0;1;0;1;\n"
+                         "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+
+  tableau.applyX(1);
+  expected = "0;0;1;0;1;\n"
+             "0;0;0;1;1;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, -IZ]"));
+
+  tableau.applyX(1);
+  expected = "0;0;1;0;1;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+
+  tableau.applyX(0);
+  expected = "0;0;1;0;0;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +IZ]"));
+}
+
+TEST_F(TestTableau, S) {
+  // S on |0> is |0>, which is stabilized by +Z
+  tableau.applyS(0);
+  std::string expected = "0;0;1;0;0;\n"
+                         "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +IZ]"));
+
+  // S on |1> is i|1>, which is stabilized by -Z
+  tableau.applyX(0);
+  tableau.applyS(0);
+  expected = "0;0;1;0;1;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+
+  // S on |+> is |R> = 1/sqrt(2) (|0> + i|1>), which is stabilized by Y
+  tableau.applyH(1);
+  tableau.applyS(1);
+  expected = "0;0;1;0;1;\n"
+             "0;1;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IY]"));
+
+  tableau.applySdag(1);
+  expected = "0;0;1;0;1;\n"
+             "0;1;0;0;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IX]"));
+}
+
+TEST_F(TestTableau, Z) {
+  // Z on |0> is |0>, which is stabilized by +Z
+  tableau.applyZ(0);
+  std::string expected = "0;0;1;0;0;\n"
+                         "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +IZ]"));
+
+  // Z on |1> is -|1>, which is stabilized by -Z
+  tableau.applyX(0);
+  tableau.applyZ(0);
+  expected = "0;0;1;0;1;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+
+  // Z on |+> is |->, which is stabilized by -X
+  tableau.applyH(1);
+  tableau.applyZ(1);
+  expected = "0;0;1;0;1;\n"
+             "0;1;0;0;1;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, -IX]"));
+}
+
+TEST_F(TestTableau, Sx) {
+  // Applying two Sx gates on |0> is equivalent to applying an X gate
+  tableau.applySx(0);
+  tableau.applySx(0);
+
+  std::string expected = "0;0;1;0;1;\n"
+                         "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+
+  tableau.applySxdag(0);
+  tableau.applySxdag(0);
+  expected = "0;0;1;0;0;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +IZ]"));
+}
+
+TEST_F(TestTableau, Y) {
+  // Y on |0> is i|1>, which is stabilized by -Z
+  tableau.applyY(0);
+  std::string expected = "0;0;1;0;1;\n"
+                         "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+
+  // Y on |1> is -i|0>, which is stabilized by +Z
+  tableau.applyX(1);
+  tableau.applyY(1);
+  expected = "0;0;1;0;1;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +IZ]"));
+}
+
+TEST_F(TestTableau, CX) {
+  // CX is stabilized by +ZI, +ZZ
+  tableau.applyCX(0, 1);
+  std::string expected = "0;0;1;0;0;\n"
+                         "0;0;1;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +ZZ]"));
+
+  // undo CX
+  tableau.applyCX(0, 1);
+  expected = "0;0;1;0;0;\n"
+             "0;0;0;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+ZI, +IZ]"));
+
+  // CX on |10> is |11>, which is stabilized by -ZI, +ZZ
+  tableau.applyX(0);
+  tableau.applyCX(0, 1);
+  expected = "0;0;1;0;1;\n"
+             "0;0;1;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[-ZI, +ZZ]"));
+}
+
+TEST_F(TestTableau, BellState) {
+  // |00> + |11> is stabilized by +XX, +ZZ
+  tableau.applyH(0);
+  tableau.applyCX(0, 1);
+  const std::string expected = "1;1;0;0;0;\n"
+                               "0;0;1;1;0;\n";
+  EXPECT_EQ(tableau, Tableau(expected));
+  EXPECT_EQ(tableau, Tableau("[+XX, +ZZ]"));
+}
+
+TEST_F(TestTableau, CircuitTranslation) {
+  using namespace dd::literals;
+
+  qc::QuantumComputation qc(2U);
+  qc.x(0);
+  qc.y(0);
+  qc.z(0);
+  qc.h(0);
+  qc.s(0);
+  qc.sdag(0);
+  qc.sx(0);
+  qc.sxdag(0);
+  qc.x(1, 0_pc);
+  qc.y(1, 0_pc);
+  qc.z(1, 0_pc);
+  qc.swap(0, 1);
+
+  auto compOP = std::make_unique<qc::CompoundOperation>(2U);
+  compOP->emplace_back<qc::StandardOperation>(2U, 0, qc::H);
+  compOP->emplace_back<qc::StandardOperation>(2U, 0_pc, 1, qc::X);
+  qc.emplace_back(compOP);
+
+  EXPECT_NO_THROW(tableau = cs::Tableau(qc););
+}
+
+TEST_F(TestTableau, UnsupportedOperations) {
+  using namespace dd::literals;
+
+  qc::QuantumComputation qc(3U);
+
+  // three-qubit operation not supported
+  qc.x(0, {1_pc, 2_pc});
+  EXPECT_THROW(tableau = cs::Tableau(qc), std::runtime_error);
+
+  // single-qubit gate not supported
+  qc.clear();
+  qc.t(0);
+  EXPECT_THROW(tableau = cs::Tableau(qc), std::runtime_error);
+
+  // controlled two-qubit gate not supported
+  qc.clear();
+  qc.s(0, 1_pc);
+  EXPECT_THROW(tableau = cs::Tableau(qc), std::runtime_error);
+}
+
+TEST_F(TestTableau, BVAccess) {
+  const auto bv0 = 0b01;
+  const auto bv1 = 0b10;
+  const auto bv2 = 0b00;
+
+  tableau.populateTableauFrom(bv0, 2, 0);
+  tableau.populateTableauFrom(bv1, 2, 1);
+  tableau.populateTableauFrom(bv2, 2, 2);
+  tableau.populateTableauFrom(bv2, 2, 3);
+  tableau.populateTableauFrom(bv2, 2, 4);
+
+  const auto col0 = tableau.getBVFrom(0);
+  const auto col1 = tableau.getBVFrom(1);
+  const auto col2 = tableau.getBVFrom(2);
+  const auto col3 = tableau.getBVFrom(3);
+  const auto col4 = tableau.getBVFrom(4);
+
+  EXPECT_EQ(col0, bv0);
+  EXPECT_EQ(col1, bv1);
+  EXPECT_EQ(col2, bv2);
+  EXPECT_EQ(col3, bv2);
+  EXPECT_EQ(col4, bv2);
+
+  const std::string expected = "1;0;0;0;0;\n"
+                               "0;1;0;0;0;\n";
+
+  EXPECT_EQ(tableau, Tableau(expected));
+}
+
+TEST_F(TestTableau, LargeBV) {
+  // Assert that a tableau for 128 qubits can be properly created
+  tableau = Tableau(128);
+  for (std::size_t i = 0U; i < 128U; ++i) {
+    EXPECT_EQ(tableau.getBVFrom<128>(128U + i), std::bitset<128>().set(i));
+  }
+
+  // Set the phase for all qubits to 1
+  tableau.populateTableauFrom(std::bitset<128>().flip(), 128, 256);
+  for (std::size_t i = 0U; i < 128U; ++i) {
+    EXPECT_EQ(tableau[i][256], 1U);
+  }
+}
+
+TEST_F(TestTableau, TableauIO) {
+  const auto filename = "tableau.txt";
+  tableau.dump(filename);
+  auto tableau2 = Tableau{};
+  tableau2.import(filename);
+  EXPECT_EQ(tableau, tableau2);
+}
+
+} // namespace cs