diff --git a/include/mqt-core/zx/FunctionalityConstruction.hpp b/include/mqt-core/zx/FunctionalityConstruction.hpp
index 08411a13af..cabec855ec 100644
--- a/include/mqt-core/zx/FunctionalityConstruction.hpp
+++ b/include/mqt-core/zx/FunctionalityConstruction.hpp
@@ -85,8 +85,14 @@ class FunctionalityConstruction {
                       EdgeType type = EdgeType::Simple);
   static void addCphase(ZXDiagram& diag, const PiExpression& phase, Qubit ctrl,
                         Qubit target, std::vector<Vertex>& qubits);
+  static void addMcphase(ZXDiagram& diag, const PiExpression& phase,
+                         const std::vector<Qubit>& controls, Qubit target,
+                         std::vector<Vertex>& qubits);
   static void addSwap(ZXDiagram& diag, Qubit target, Qubit target2,
                       std::vector<Vertex>& qubits);
+  static void addMcswap(ZXDiagram& diag, const std::vector<Qubit>& controls,
+                        Qubit target, Qubit target2,
+                        std::vector<Vertex>& qubits);
   static void
   addRzz(ZXDiagram& diag, const PiExpression& phase, Qubit target,
          Qubit target2, std::vector<Vertex>& qubits,
@@ -99,6 +105,11 @@ class FunctionalityConstruction {
   addRzx(ZXDiagram& diag, const PiExpression& phase, Qubit target,
          Qubit target2, std::vector<Vertex>& qubits,
          const std::optional<double>& unconvertedPhase = std::nullopt);
+  static void
+  addMcrzz(ZXDiagram& diag, const PiExpression& phase,
+           const std::vector<Qubit>& controls, const Qubit target,
+           const Qubit target2, std::vector<Vertex>& qubits,
+           const std::optional<double>& unconvertedPhase = std::nullopt);
   static void addDcx(ZXDiagram& diag, Qubit qubit1, Qubit qubit2,
                      std::vector<Vertex>& qubits);
   static void
@@ -113,6 +124,24 @@ class FunctionalityConstruction {
                const std::optional<double>& unconvertedBeta = std::nullopt);
   static void addCcx(ZXDiagram& diag, Qubit ctrl0, Qubit ctrl1, Qubit target,
                      std::vector<Vertex>& qubits);
+  static void addCcz(ZXDiagram& diag, Qubit ctrl0, Qubit ctrl1, Qubit target,
+                     std::vector<Vertex>& qubits);
+  static void addCrx(ZXDiagram& diag, const PiExpression& phase,
+                     const Qubit control, const Qubit target,
+                     std::vector<Vertex>& qubits);
+  static void addMcrx(ZXDiagram& diag, const PiExpression& phase,
+                      std::vector<Qubit> controls, const Qubit target,
+                      std::vector<Vertex>& qubits);
+  static void addCrz(ZXDiagram& diag, const PiExpression& phase,
+                     const Qubit control, const Qubit target,
+                     std::vector<Vertex>& qubits);
+  static void addMcrz(ZXDiagram& diag, const PiExpression& phase,
+                      std::vector<Qubit> controls, const Qubit target,
+                      std::vector<Vertex>& qubits);
+  static void addMcx(ZXDiagram& diag, std::vector<Qubit> controls,
+                     const Qubit target, std::vector<Vertex>& qubits);
+  static void addMcz(ZXDiagram& diag, std::vector<Qubit> controls,
+                     const Qubit target, std::vector<Vertex>& qubits);
   static op_it parseOp(ZXDiagram& diag, op_it it, op_it end,
                        std::vector<Vertex>& qubits, const qc::Permutation& p);
   static op_it parseCompoundOp(ZXDiagram& diag, op_it it, op_it end,
diff --git a/src/zx/FunctionalityConstruction.cpp b/src/zx/FunctionalityConstruction.cpp
index 23697afaef..8efea0e1bb 100644
--- a/src/zx/FunctionalityConstruction.cpp
+++ b/src/zx/FunctionalityConstruction.cpp
@@ -140,6 +140,32 @@ void FunctionalityConstruction::addCphase(ZXDiagram& diag,
   addZSpider(diag, target, qubits, newPhase);
 }
 
+void FunctionalityConstruction::addMcphase(ZXDiagram& diag,
+                                           const PiExpression& phase,
+                                           const std::vector<Qubit>& controls,
+                                           const Qubit target,
+                                           std::vector<Vertex>& qubits) {
+  auto newConst = phase.getConst() / 2;
+  auto newPhase = phase / 2.0;
+  newPhase.setConst(newConst);
+  addZSpider(diag, target, qubits, newPhase);
+  addMcx(diag, controls, target, qubits);
+  addZSpider(diag, target, qubits, -newPhase);
+  addMcx(diag, controls, target, qubits);
+  switch (controls.size()) {
+  case 1:
+    addZSpider(diag, controls[0], qubits, newPhase);
+    return;
+  case 2:
+    addCphase(diag, newPhase, controls[0], controls[1], qubits);
+    return;
+  default:
+    addMcphase(diag, newPhase,
+               std::vector<Qubit>(controls.begin(), controls.end() - 1),
+               controls.back(), qubits);
+  }
+}
+
 void FunctionalityConstruction::addRzz(
     ZXDiagram& diag, const PiExpression& phase, const Qubit target,
     const Qubit target2, std::vector<Vertex>& qubits,
@@ -207,6 +233,20 @@ void FunctionalityConstruction::addRzx(
   }
 }
 
+void FunctionalityConstruction::addMcrzz(
+    ZXDiagram& diag, const PiExpression& phase,
+    const std::vector<Qubit>& controls, const Qubit target, const Qubit target2,
+    std::vector<Vertex>& qubits,
+    const std::optional<double>& unconvertedPhase) {
+  addRzz(diag, phase, target, target2, qubits, unconvertedPhase);
+  addMcx(diag, controls, target, qubits);
+  addRzz(diag, -phase, target, target2, qubits,
+         unconvertedPhase.has_value()
+             ? std::optional<double>(-unconvertedPhase.value())
+             : std::nullopt);
+  addMcx(diag, controls, target, qubits);
+}
+
 void FunctionalityConstruction::addDcx(ZXDiagram& diag, const Qubit qubit1,
                                        const Qubit qubit2,
                                        std::vector<Vertex>& qubits) {
@@ -286,6 +326,22 @@ void FunctionalityConstruction::addSwap(ZXDiagram& diag, const Qubit target,
   qubits[c] = t1;
 }
 
+void FunctionalityConstruction::addMcswap(ZXDiagram& diag,
+                                          const std::vector<Qubit>& controls,
+                                          const Qubit target,
+                                          const Qubit target2,
+                                          std::vector<Vertex>& qubits) {
+  std::vector controls1 = controls;
+  std::vector controls2 = controls;
+
+  controls1.emplace_back(target);
+  controls2.emplace_back(target2);
+
+  addMcx(diag, controls1, target2, qubits);
+  addMcx(diag, controls2, target, qubits);
+  addMcx(diag, controls1, target2, qubits);
+}
+
 void FunctionalityConstruction::addCcx(ZXDiagram& diag, const Qubit ctrl0,
                                        const Qubit ctrl1, const Qubit target,
                                        std::vector<Vertex>& qubits) {
@@ -307,6 +363,125 @@ void FunctionalityConstruction::addCcx(ZXDiagram& diag, const Qubit ctrl0,
   addCnot(diag, ctrl0, ctrl1, qubits);
 }
 
+void FunctionalityConstruction::addCcz(ZXDiagram& diag, const Qubit ctrl0,
+                                       const Qubit ctrl1, const Qubit target,
+                                       std::vector<Vertex>& qubits) {
+
+  addCnot(diag, ctrl1, target, qubits);
+  addZSpider(diag, target, qubits, PiExpression(PiRational(-1, 4)));
+  addCnot(diag, ctrl0, target, qubits);
+  addZSpider(diag, target, qubits, PiExpression(PiRational(1, 4)));
+  addCnot(diag, ctrl1, target, qubits);
+  addZSpider(diag, ctrl1, qubits, PiExpression(PiRational(1, 4)));
+  addZSpider(diag, target, qubits, PiExpression(PiRational(-1, 4)));
+  addCnot(diag, ctrl0, target, qubits);
+  addZSpider(diag, target, qubits, PiExpression(PiRational(1, 4)));
+  addCnot(diag, ctrl0, ctrl1, qubits);
+  addZSpider(diag, ctrl0, qubits, PiExpression(PiRational(1, 4)));
+  addZSpider(diag, ctrl1, qubits, PiExpression(PiRational(-1, 4)));
+  addZSpider(diag, target, qubits, PiExpression(PiRational(0, 1)),
+             EdgeType::Hadamard);
+  addCnot(diag, ctrl0, ctrl1, qubits);
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+}
+
+void FunctionalityConstruction::addCrx(ZXDiagram& diag,
+                                       const PiExpression& phase,
+                                       const Qubit control, const Qubit target,
+                                       std::vector<Vertex>& qubits) {
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+  addCrz(diag, phase, control, target, qubits);
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+}
+
+void FunctionalityConstruction::addMcrx(ZXDiagram& diag,
+                                        const PiExpression& phase,
+                                        std::vector<Qubit> controls,
+                                        const Qubit target,
+                                        std::vector<Vertex>& qubits) {
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+  addMcrz(diag, phase, std::move(controls), target, qubits);
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+}
+
+void FunctionalityConstruction::addCrz(ZXDiagram& diag,
+                                       const PiExpression& phase,
+                                       const Qubit q0, const Qubit q1,
+                                       std::vector<Vertex>& qubits) {
+  // CRZ decomposition uses reversed CNOT direction
+  addCnot(diag, q1, q0, qubits);
+  addZSpider(diag, q0, qubits, -phase / 2);
+  addZSpider(diag, q1, qubits, phase / 2);
+  addCnot(diag, q1, q0, qubits);
+}
+
+void FunctionalityConstruction::addMcrz(ZXDiagram& diag,
+                                        const PiExpression& phase,
+                                        std::vector<Qubit> controls,
+                                        const Qubit target,
+                                        std::vector<Vertex>& qubits) {
+
+  const Qubit nextControl = controls.back();
+  controls.pop_back();
+
+  addCrz(diag, phase / 2, nextControl, target, qubits);
+  addMcx(diag, controls, target, qubits);
+  addCrz(diag, -phase / 2, nextControl, target, qubits);
+  addMcx(diag, controls, target, qubits);
+}
+
+void FunctionalityConstruction::addMcx(ZXDiagram& diag,
+                                       std::vector<Qubit> controls,
+                                       const Qubit target,
+                                       std::vector<Vertex>& qubits) {
+
+  switch (controls.size()) {
+  case 0:
+    addXSpider(diag, target, qubits, PiExpression(PiRational(1, 1)));
+    return;
+  case 1:
+    addCnot(diag, controls.front(), target, qubits);
+    return;
+  case 2:
+    addCcx(diag, controls.front(), controls.back(), target, qubits);
+    return;
+  default:
+    const auto half = static_cast<std::ptrdiff_t>((controls.size() + 1) / 2);
+    const std::vector<Qubit> first(controls.begin(), controls.begin() + half);
+    const std::vector<Qubit> second(controls.begin() + half, controls.end());
+
+    addRx(diag, PiExpression(PiRational(1, 4)), target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addMcx(diag, first, target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addRx(diag, PiExpression(-PiRational(1, 4)), target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addMcx(diag, second, target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addRx(diag, PiExpression(PiRational(1, 4)), target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addMcx(diag, first, target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addRx(diag, PiExpression(-PiRational(1, 4)), target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    addMcx(diag, second, target, qubits);
+    addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+    const Qubit lastControl = controls.back();
+    controls.pop_back();
+    addMcphase(diag, PiExpression(PiRational(1, 2)), controls, lastControl,
+               qubits);
+  }
+}
+
+void FunctionalityConstruction::addMcz(ZXDiagram& diag,
+                                       std::vector<Qubit> controls,
+                                       const Qubit target,
+                                       std::vector<Vertex>& qubits) {
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+  addMcx(diag, std::move(controls), target, qubits);
+  addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+}
+
 FunctionalityConstruction::op_it
 FunctionalityConstruction::parseOp(ZXDiagram& diag, op_it it, op_it end,
                                    std::vector<Vertex>& qubits,
@@ -429,21 +604,27 @@ FunctionalityConstruction::parseOp(ZXDiagram& diag, op_it it, op_it end,
     case qc::OpType::RXX: {
       const auto target2 = static_cast<Qubit>(p.at(op->getTargets()[1]));
       const auto& phase = parseParam(op.get(), 0);
+      addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
       if (phase.isConstant()) {
         addRxx(diag, phase, target, target2, qubits, op->getParameter().at(0));
       } else {
         addRxx(diag, phase, target, target2, qubits);
       }
+      addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
       break;
     }
     case qc::OpType::RZX: {
       const auto target2 = static_cast<Qubit>(p.at(op->getTargets()[1]));
       const auto& phase = parseParam(op.get(), 0);
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
       if (phase.isConstant()) {
         addRzx(diag, phase, target, target2, qubits, op->getParameter().at(0));
       } else {
         addRzx(diag, phase, target, target2, qubits);
       }
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
       break;
     }
     case qc::OpType::RYY: {
@@ -538,10 +719,15 @@ FunctionalityConstruction::parseOp(ZXDiagram& diag, op_it it, op_it end,
                    qubits[static_cast<std::size_t>(target)],
                    EdgeType::Hadamard);
       break;
+    case qc::OpType::RZ:
+      addCrz(diag, parseParam(op.get(), 0), ctrl, target, qubits);
+      break;
 
+    case qc::OpType::RX:
+      addCrx(diag, parseParam(op.get(), 0), ctrl, target, qubits);
+      break;
     case qc::OpType::I:
       break;
-
     case qc::OpType::P:
       addCphase(diag, parseParam(op.get(), 0), ctrl, target, qubits);
       break;
@@ -565,7 +751,7 @@ FunctionalityConstruction::parseOp(ZXDiagram& diag, op_it it, op_it end,
       throw ZXException("Unsupported Controlled Operation: " +
                         qc::toString(op->getType()));
     }
-  } else if (op->getNcontrols() == 2) {
+  } else if (op->getNcontrols() == 2 && op->getNtargets() == 1) {
     // three-qubit controlled gates (ccx or ccz)
     Qubit ctrl0 = 0;
     Qubit ctrl1 = 0;
@@ -582,14 +768,119 @@ FunctionalityConstruction::parseOp(ZXDiagram& diag, op_it it, op_it end,
     case qc::OpType::X:
       addCcx(diag, ctrl0, ctrl1, target, qubits);
       break;
-
     case qc::OpType::Z:
+      addCcz(diag, ctrl0, ctrl1, target, qubits);
+      break;
+    case qc::OpType::P:
+      addMcphase(diag, parseParam(op.get(), 0), {ctrl0, ctrl1}, target, qubits);
+      break;
+    case qc::OpType::T:
+      addMcphase(diag, PiExpression{PiRational(1, 4)}, {ctrl0, ctrl1}, target,
+                 qubits);
+      break;
+    case qc::OpType::S:
+      addMcphase(diag, PiExpression{PiRational(1, 2)}, {ctrl0, ctrl1}, target,
+                 qubits);
+      break;
+    case qc::OpType::RZ:
+      addMcrz(diag, parseParam(op.get(), 0), {ctrl0, ctrl1}, target, qubits);
+      break;
+    case qc::OpType::RX:
+      addMcrx(diag, parseParam(op.get(), 0), {ctrl0, ctrl1}, target, qubits);
+      break;
+    default:
+      throw ZXException("Unsupported Multi-control operation: " +
+                        qc::toString(op->getType()));
+    }
+  } else if (op->getNtargets() == 1) {
+    const auto target = static_cast<Qubit>(p.at(op->getTargets().front()));
+    std::vector<Qubit> controls;
+    controls.reserve(op->getNcontrols());
+    for (const auto& ctrl : op->getControls()) {
+      controls.emplace_back(static_cast<Qubit>(p.at(ctrl.qubit)));
+    }
+    switch (op->getType()) {
+    case qc::OpType::X:
+      addMcx(diag, controls, target, qubits);
+      break;
+    case qc::OpType::Z:
+      addMcz(diag, controls, target, qubits);
+      break;
+    case qc::OpType::P:
+      addMcphase(diag, parseParam(op.get(), 0), controls, target, qubits);
+      break;
+    case qc::OpType::T:
+      addMcphase(diag, PiExpression{PiRational(1, 4)}, controls, target,
+                 qubits);
+      break;
+    case qc::OpType::S:
+      addMcphase(diag, PiExpression{PiRational(1, 2)}, controls, target,
+                 qubits);
+      break;
+    case qc::OpType::RZ:
+      addMcrz(diag, parseParam(op.get(), 0), controls, target, qubits);
+      break;
+    case qc::OpType::RX:
+      addMcrx(diag, parseParam(op.get(), 0), controls, target, qubits);
+      break;
+    default:
+      throw ZXException("Unsupported Multi-control operation (" +
+                        std::to_string(op->getNcontrols()) + " ctrls)" +
+                        qc::toString(op->getType()));
+    }
+  } else if (op->getNtargets() == 2) {
+    // At this point, op must have getNtargets() >= 2 (all 1-target cases
+    // handled above)
+    const auto target = static_cast<Qubit>(p.at(op->getTargets().front()));
+    const auto target2 = static_cast<Qubit>(p.at(op->getTargets()[1]));
+    std::vector<Qubit> controls;
+    controls.reserve(op->getNcontrols());
+    for (const auto& ctrl : op->getControls()) {
+      controls.emplace_back(static_cast<Qubit>(p.at(ctrl.qubit)));
+    }
+    switch (op->getType()) {
+    case qc::OpType::SWAP:
+      addMcswap(diag, controls, target, target2, qubits);
+      break;
+    case qc::OpType::RZZ: {
+      const auto& phase = parseParam(op.get(), 0);
+      if (phase.isConstant()) {
+        addMcrzz(diag, phase, controls, target, target2, qubits,
+                 op->getParameter().at(0));
+      } else {
+        addMcrzz(diag, phase, controls, target, target2, qubits);
+      }
+      break;
+    }
+    case qc::OpType::RXX: {
+      const auto& phase = parseParam(op.get(), 0);
       addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
-      addCcx(diag, ctrl0, ctrl1, target, qubits);
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
+      if (phase.isConstant()) {
+        addMcrzz(diag, phase, controls, target, target2, qubits,
+                 op->getParameter().at(0));
+      } else {
+        addMcrzz(diag, phase, controls, target, target2, qubits);
+      }
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
       addZSpider(diag, target, qubits, PiExpression(), EdgeType::Hadamard);
       break;
+    }
+    case qc::OpType::RZX: {
+      const auto& phase = parseParam(op.get(), 0);
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
+      if (phase.isConstant()) {
+        addMcrzz(diag, phase, controls, target, target2, qubits,
+                 op->getParameter().at(0));
+      } else {
+        addMcrzz(diag, phase, controls, target, target2, qubits);
+      }
+      addZSpider(diag, target2, qubits, PiExpression(), EdgeType::Hadamard);
+      break;
+    }
     default:
-      throw ZXException("Unsupported Multi-control operation: " +
+      throw ZXException("Unsupported Multi-control operation (" +
+                        std::to_string(op->getNcontrols()) + " ctrls)" +
                         qc::toString(op->getType()));
     }
   } else {
@@ -701,15 +992,31 @@ bool FunctionalityConstruction::transformableToZX(const qc::Operation* op) {
     case qc::OpType::S:
     case qc::OpType::Tdg:
     case qc::OpType::Sdg:
+    case qc::OpType::RX:
+    case qc::OpType::RZ:
       return true;
-
     default:
       return false;
     }
-  } else if (op->getNcontrols() == 2) {
+  } else if (op->getNtargets() == 1) {
     switch (op->getType()) {
     case qc::OpType::X:
     case qc::OpType::Z:
+    case qc::OpType::P:
+    case qc::OpType::T:
+    case qc::OpType::S:
+    case qc::OpType::RZ:
+    case qc::OpType::RX:
+      return true;
+    default:
+      return false;
+    }
+  } else if (op->getNtargets() == 2) {
+    switch (op->getType()) {
+    case qc::OpType::SWAP:
+    case qc::OpType::RZZ:
+    case qc::OpType::RXX:
+    case qc::OpType::RZX:
       return true;
     default:
       return false;
diff --git a/test/zx/test_zx_functionality.cpp b/test/zx/test_zx_functionality.cpp
index 6ef29b69e8..035f5c0961 100644
--- a/test/zx/test_zx_functionality.cpp
+++ b/test/zx/test_zx_functionality.cpp
@@ -28,12 +28,33 @@
 #include <sstream>
 #include <string>
 #include <utility>
+#include <vector>
 
 namespace zx {
 class ZXFunctionalityTest : public ::testing::Test {
 public:
   qc::QuantumComputation qc;
 };
+namespace {
+void checkEquivalence(const qc::QuantumComputation& qc1,
+                      const qc::QuantumComputation& qc2,
+                      const std::vector<qc::Qubit>& qubits) {
+  EXPECT_TRUE(FunctionalityConstruction::transformableToZX(&qc1));
+  EXPECT_TRUE(FunctionalityConstruction::transformableToZX(&qc2));
+  EXPECT_EQ(qc1.getNqubits(), qc2.getNqubits());
+
+  auto d1 = FunctionalityConstruction::buildFunctionality(&qc1);
+  auto d2 = FunctionalityConstruction::buildFunctionality(&qc2);
+  d1.concat(d2.invert());
+  fullReduce(d1);
+  EXPECT_TRUE(d1.isIdentity());
+  EXPECT_TRUE(d1.globalPhaseIsZero());
+  for (const auto q : qubits) {
+    ASSERT_LT(q, qc1.getNqubits());
+    EXPECT_TRUE(d1.connected(d1.getInput(q), d1.getOutput(q)));
+  }
+}
+} // namespace
 
 TEST_F(ZXFunctionalityTest, parseQasm) {
   const std::string testfile = "OPENQASM 2.0;"
@@ -91,70 +112,29 @@ TEST_F(ZXFunctionalityTest, complexCircuit) {
   std::stringstream ss{};
   ss << "// i 1 0 2\n"
      << "// o 0 1 2\n"
-     << "OPENQASM 2.0;"
-     << "include \"qelib1.inc\";"
-     << "qreg q[3];"
-     << "sx q[0];"
-     << "sxdg q[0];"
-     << "h q[0];"
-     << "cx q[0],q[1];"
-     << "z q[1];"
-     << "x q[2];"
-     << "y q[0];"
-     << "rx(pi/4) q[0];"
-     << "rz(0.1) q[1];"
-     << "p(0.1) q[1];"
-     << "ry(pi/4) q[2];"
-     << "t q[0];"
-     << "s q[2];"
-     << "u2(pi/4, pi/4) q[1];"
-     << "u3(pi/4, pi/4, pi/4) q[2];"
-     << "barrier q[0],q[1],q[2];"
-     << "swap q[0],q[1];"
-     << "cz q[1],q[2];"
-     << "cp(pi/4) q[0],q[1];"
-     << "ctrl(2) @ x q[0],q[1],q[2];"
-     << "ctrl(2) @ z q[1],q[2],q[0];"
-     << "cp(pi/2) q[0], q[1];"
-     << "cp(pi/4) q[0], q[1];"
-     << "cp(pi/8) q[0], q[1];"
-     << "rzz(pi/4) q[0], q[1];"
-     << "rxx(pi/4) q[0], q[1];"
-     << "ryy(pi/4) q[0], q[1];"
-     << "rzx(pi/4) q[0], q[1];"
-     << "ecr q[0], q[1];"
-     << "dcx q[0], q[1];"
-     << "r(pi/8, pi/4) q[2];"
-     << "r(-pi/8, pi/4) q[2];"
-     << "dcx q[1], q[0];"
-     << "ecr q[0], q[1];"
-     << "rzx(-pi/4) q[0], q[1];"
-     << "ryy(-pi/4) q[0], q[1];"
-     << "rxx(-pi/4) q[0], q[1];"
-     << "rzz(-pi/4) q[0], q[1];"
-     << "cp(-pi/8) q[0], q[1];"
-     << "cp(-pi/4) q[0], q[1];"
-     << "cp(-pi/2) q[0], q[1];"
-     << "ctrl(2) @ z q[1],q[2],q[0];"
-     << "ctrl(2) @ x q[0],q[1],q[2];"
-     << "cp(-pi/4) q[0],q[1];"
-     << "cz q[1],q[2];"
-     << "cx q[1],q[0];"
-     << "cx q[0],q[1];"
-     << "cx q[1],q[0];"
-     << "u3(-pi/4,-pi/4,-pi/4) q[2];"
-     << "u2(-5*pi/4,3*pi/4) q[1];"
-     << "sdg q[2];"
-     << "tdg q[0];"
-     << "ry(-pi/4) q[2];"
-     << "p(-0.1) q[1];"
-     << "rz(-0.1) q[1];"
-     << "rx(-pi/4) q[0];"
-     << "y q[0];"
-     << "x q[2];"
-     << "z q[1];"
-     << "cx q[0],q[1];"
-     << "h q[0];\n";
+     << "OPENQASM 2.0;" << "include \"qelib1.inc\";" << "qreg q[3];"
+     << "sx q[0];" << "sxdg q[0];" << "h q[0];" << "cx q[0],q[1];" << "z q[1];"
+     << "x q[2];" << "y q[0];" << "rx(pi/4) q[0];" << "rz(0.1) q[1];"
+     << "p(0.1) q[1];" << "ry(pi/4) q[2];" << "t q[0];" << "s q[2];"
+     << "u2(pi/4, pi/4) q[1];" << "u3(pi/4, pi/4, pi/4) q[2];"
+     << "barrier q[0],q[1],q[2];" << "swap q[0],q[1];" << "cz q[1],q[2];"
+     << "cp(pi/4) q[0],q[1];" << "ctrl(2) @ x q[0],q[1],q[2];"
+     << "ctrl(2) @ z q[1],q[2],q[0];" << "cp(pi/2) q[0], q[1];"
+     << "cp(pi/4) q[0], q[1];" << "cp(pi/8) q[0], q[1];"
+     << "rzz(pi/4) q[0], q[1];" << "rxx(pi/4) q[0], q[1];"
+     << "ryy(pi/4) q[0], q[1];" << "rzx(pi/4) q[0], q[1];" << "ecr q[0], q[1];"
+     << "dcx q[0], q[1];" << "r(pi/8, pi/4) q[2];" << "r(-pi/8, pi/4) q[2];"
+     << "dcx q[1], q[0];" << "ecr q[0], q[1];" << "rzx(-pi/4) q[0], q[1];"
+     << "ryy(-pi/4) q[0], q[1];" << "rxx(-pi/4) q[0], q[1];"
+     << "rzz(-pi/4) q[0], q[1];" << "cp(-pi/8) q[0], q[1];"
+     << "cp(-pi/4) q[0], q[1];" << "cp(-pi/2) q[0], q[1];"
+     << "ctrl(2) @ z q[1],q[2],q[0];" << "ctrl(2) @ x q[0],q[1],q[2];"
+     << "cp(-pi/4) q[0],q[1];" << "cz q[1],q[2];" << "cx q[1],q[0];"
+     << "cx q[0],q[1];" << "cx q[1],q[0];" << "u3(-pi/4,-pi/4,-pi/4) q[2];"
+     << "u2(-5*pi/4,3*pi/4) q[1];" << "sdg q[2];" << "tdg q[0];"
+     << "ry(-pi/4) q[2];" << "p(-0.1) q[1];" << "rz(-0.1) q[1];"
+     << "rx(-pi/4) q[0];" << "y q[0];" << "x q[2];" << "z q[1];"
+     << "cx q[0],q[1];" << "h q[0];\n";
   qc = qasm3::Importer::import(ss);
 
   EXPECT_TRUE(FunctionalityConstruction::transformableToZX(&qc));
@@ -187,7 +167,7 @@ TEST_F(ZXFunctionalityTest, nestedCompoundGate) {
 }
 
 TEST_F(ZXFunctionalityTest, Phase) {
-  using namespace qc::literals;
+
   qc = qc::QuantumComputation(2);
   qc.p(PI / 4, 0);
   qc.cp(PI / 4, 1, 0);
@@ -219,18 +199,149 @@ TEST_F(ZXFunctionalityTest, Compound) {
   EXPECT_TRUE(diag.isIdentity());
 }
 
-TEST_F(ZXFunctionalityTest, UnsupportedMultiControl) {
-  using namespace qc::literals;
+TEST_F(ZXFunctionalityTest, CRZ) {
+  qc = qc::QuantumComputation(2);
+  qc.crz(PI / 2, 0, 1);
+
+  auto qcPrime = qc::QuantumComputation(2);
+
+  qcPrime.cx(0, 1);
+  qcPrime.rz(-PI / 4, 1);
+  qcPrime.cx(0, 1);
+  qcPrime.rz(PI / 4, 1);
+
+  checkEquivalence(qc, qcPrime, {0, 1});
+}
+
+TEST_F(ZXFunctionalityTest, CCZ) {
+
+  qc = qc::QuantumComputation(3);
+  qc.mcz({1, 2}, 0);
+
+  auto qcPrime = qc::QuantumComputation(3);
+  qcPrime.h(0);
+  qcPrime.mcx({1, 2}, 0);
+  qcPrime.h(0);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2});
+}
+
+TEST_F(ZXFunctionalityTest, MCX) {
+
   qc = qc::QuantumComputation(4);
   qc.mcx({1, 2, 3}, 0);
-  EXPECT_FALSE(FunctionalityConstruction::transformableToZX(&qc));
-  EXPECT_THROW(const ZXDiagram diag =
-                   FunctionalityConstruction::buildFunctionality(&qc),
-               ZXException);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcx({3, 2, 1}, 0);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCX0) {
+
+  qc = qc::QuantumComputation(1);
+  qc.mcx({}, 0);
+
+  auto qcPrime = qc::QuantumComputation(1);
+
+  qcPrime.x(0);
+  checkEquivalence(qc, qcPrime, {0});
+}
+
+TEST_F(ZXFunctionalityTest, MCX1) {
+
+  qc = qc::QuantumComputation(2);
+  qc.mcx({1}, 0);
+
+  auto qcPrime = qc::QuantumComputation(2);
+
+  qcPrime.cx(1, 0);
+
+  checkEquivalence(qc, qcPrime, {0, 1});
+}
+
+TEST_F(ZXFunctionalityTest, MCZ) {
+
+  qc = qc::QuantumComputation(4);
+  qc.mcz({1, 2, 3}, 0);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcz({1, 2, 3}, 0);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCZ0) {
+
+  qc = qc::QuantumComputation(1);
+  qc.mcz({}, 0);
+
+  auto qcPrime = qc::QuantumComputation(1);
+  qcPrime.z(0);
+
+  checkEquivalence(qc, qcPrime, {0});
+}
+
+TEST_F(ZXFunctionalityTest, MCZ1) {
+
+  qc = qc::QuantumComputation(2);
+  qc.mcz({1}, 0);
+
+  auto qcPrime = qc::QuantumComputation(2);
+  qcPrime.cz(1, 0);
+
+  checkEquivalence(qc, qcPrime, {0, 1});
+}
+
+TEST_F(ZXFunctionalityTest, MCZ2) {
+
+  qc = qc::QuantumComputation(4);
+  qc.mcz({1, 2}, 0);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.h(0);
+  qcPrime.mcx({1, 2}, 0);
+  qcPrime.h(0);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCRZ) {
+
+  qc = qc::QuantumComputation(4);
+  qc.mcrz(PI / 4, {1, 2, 3}, 0);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.h(0);
+  qcPrime.mcrx(PI / 4, {1, 2, 3}, 0);
+  qcPrime.h(0);
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCRZ0) {
+
+  qc = qc::QuantumComputation(1);
+  qc.mcrz(PI / 4, {}, 0);
+
+  auto qcPrime = qc::QuantumComputation(1);
+  qcPrime.rz(PI / 4, 0);
+
+  checkEquivalence(qc, qcPrime, {0});
+}
+
+TEST_F(ZXFunctionalityTest, MCRZ1) {
+
+  qc = qc::QuantumComputation(2);
+  qc.mcrz(PI / 4, {1}, 0);
+
+  auto qcPrime = qc::QuantumComputation(2);
+  qcPrime.crz(PI / 4, 1, 0);
+
+  checkEquivalence(qc, qcPrime, {0, 1});
 }
 
 TEST_F(ZXFunctionalityTest, UnsupportedControl) {
-  using namespace qc::literals;
+
   qc = qc::QuantumComputation(2);
   qc.cy(1, 0);
   EXPECT_FALSE(FunctionalityConstruction::transformableToZX(&qc));
@@ -240,7 +351,7 @@ TEST_F(ZXFunctionalityTest, UnsupportedControl) {
 }
 
 TEST_F(ZXFunctionalityTest, UnsupportedControl2) {
-  using namespace qc::literals;
+
   qc = qc::QuantumComputation(3);
   qc.mcy({1, 2}, 0);
   EXPECT_FALSE(FunctionalityConstruction::transformableToZX(&qc));
@@ -262,13 +373,7 @@ TEST_F(ZXFunctionalityTest, InitialLayout) {
   qcPrime.x(1);
   qcPrime.z(0);
 
-  auto d = FunctionalityConstruction::buildFunctionality(&qc);
-  auto dPrime = FunctionalityConstruction::buildFunctionality(&qcPrime);
-
-  d.concat(dPrime);
-
-  fullReduce(d);
-  EXPECT_TRUE(d.isIdentity());
+  checkEquivalence(qc, qcPrime, {0, 1});
 }
 
 TEST_F(ZXFunctionalityTest, FromSymbolic) {
@@ -289,21 +394,13 @@ TEST_F(ZXFunctionalityTest, RZ) {
   qc.rz(PI / 8, 0);
 
   auto qcPrime = qc::QuantumComputation(1);
-  qcPrime.p(PI / 8, 0);
-
-  auto d = FunctionalityConstruction::buildFunctionality(&qc);
-  auto dPrime = FunctionalityConstruction::buildFunctionality(&qcPrime);
+  qcPrime.rz(PI / 8, 0);
 
-  d.concat(dPrime.invert());
-
-  fullReduce(d);
-  EXPECT_FALSE(d.isIdentity());
-  EXPECT_FALSE(d.globalPhaseIsZero());
-  EXPECT_TRUE(d.connected(d.getInput(0), d.getOutput(0)));
+  checkEquivalence(qc, qcPrime, {0});
 }
 
 TEST_F(ZXFunctionalityTest, ISWAP) {
-  using namespace qc::literals;
+
   qc = qc::QuantumComputation(2);
   qc.iswap(0, 1);
 
@@ -315,15 +412,7 @@ TEST_F(ZXFunctionalityTest, ISWAP) {
   qcPrime.cx(1, 0);
   qc.h(1);
 
-  auto d = FunctionalityConstruction::buildFunctionality(&qc);
-  auto dPrime = FunctionalityConstruction::buildFunctionality(&qcPrime);
-
-  d.concat(dPrime.invert());
-
-  fullReduce(d);
-  EXPECT_TRUE(d.isIdentity());
-  EXPECT_TRUE(d.globalPhaseIsZero());
-  EXPECT_TRUE(d.connected(d.getInput(0), d.getOutput(0)));
+  checkEquivalence(qc, qcPrime, {0, 1});
 }
 
 TEST_F(ZXFunctionalityTest, XXplusYY) {
@@ -349,17 +438,7 @@ TEST_F(ZXFunctionalityTest, XXplusYY) {
   qcPrime.rz(qc::PI_2, 0);
   qcPrime.rz(-beta, 1);
 
-  auto d = FunctionalityConstruction::buildFunctionality(&qc);
-
-  auto dPrime = FunctionalityConstruction::buildFunctionality(&qcPrime);
-
-  d.concat(dPrime.invert());
-
-  fullReduce(d);
-
-  EXPECT_TRUE(d.isIdentity());
-  EXPECT_TRUE(d.globalPhaseIsZero());
-  EXPECT_TRUE(d.connected(d.getInput(0), d.getOutput(0)));
+  checkEquivalence(qc, qcPrime, {0, 1});
 }
 
 TEST_F(ZXFunctionalityTest, XXminusYY) {
@@ -385,16 +464,171 @@ TEST_F(ZXFunctionalityTest, XXminusYY) {
   qcPrime.rz(qc::PI_2, 0);
   qcPrime.rz(beta, 1);
 
-  auto d = FunctionalityConstruction::buildFunctionality(&qc);
+  checkEquivalence(qc, qcPrime, {0, 1});
+}
+
+TEST_F(ZXFunctionalityTest, SWAP) {
+  qc = qc::QuantumComputation(2);
+  qc.swap(0, 1);
+
+  auto qcPrime = qc::QuantumComputation(2);
+  qcPrime.cx(1, 0);
+  qcPrime.cx(0, 1);
+  qcPrime.cx(1, 0);
+
+  checkEquivalence(qc, qcPrime, {0, 1});
+}
+
+TEST_F(ZXFunctionalityTest, CSWAP) {
+  qc = qc::QuantumComputation(3);
+  qc.mcswap({0}, 1, 2);
+
+  auto qcPrime = qc::QuantumComputation(3);
+  qcPrime.mcx({0, 1}, 2);
+  qcPrime.mcx({0, 2}, 1);
+  qcPrime.mcx({0, 1}, 2);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2});
+}
+
+TEST_F(ZXFunctionalityTest, MCSWAP) {
+  qc = qc::QuantumComputation(4);
+  qc.mcswap({0, 1}, 2, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcx({0, 1, 2}, 3);
+  qcPrime.mcx({0, 1, 3}, 2);
+  qcPrime.mcx({0, 1, 2}, 3);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCRzz) {
+  qc = qc::QuantumComputation(4);
+  qc.mcrzz(qc::PI_2 / 3, {0, 1}, 2, 3);
+  qc.mcrzz(2 * qc::PI, {0, 1}, 2, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcrzz(qc::PI_2 / 3, {0, 1}, 2, 3);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCRxx) {
+  qc = qc::QuantumComputation(4);
+  qc.mcrxx(qc::PI_2, {0, 1}, 2, 3);
+  qc.mcrzx(2 * qc::PI, {0, 1}, 2, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcrxx(qc::PI_2, {0, 1}, 2, 3);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCRzx) {
+  qc = qc::QuantumComputation(4);
+  qc.mcrzx(qc::PI_2, {0, 1}, 2, 3);
+  qc.mcrzx(2 * qc::PI, {0, 1}, 2, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcrzx(qc::PI_2, {0, 1}, 2, 3);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCRx0) {
+  qc = qc::QuantumComputation(1);
+  qc.mcrx(PI / 4, {}, 0);
+
+  auto qcPrime = qc::QuantumComputation(1);
+  qcPrime.mcrx(PI / 4, {}, 0);
+  checkEquivalence(qc, qcPrime, {0});
+}
+
+TEST_F(ZXFunctionalityTest, MCRx1) {
+  qc = qc::QuantumComputation(2);
+  qc.mcrx(PI / 4, {0}, 1);
+
+  auto qcPrime = qc::QuantumComputation(2);
+  qcPrime.mcrx(PI / 4, {0}, 1);
+  checkEquivalence(qc, qcPrime, {0, 1});
+}
+
+TEST_F(ZXFunctionalityTest, MCRx2) {
+  qc = qc::QuantumComputation(3);
+  qc.mcrx(PI / 4, {0, 1}, 2);
+
+  auto qcPrime = qc::QuantumComputation(3);
+  qcPrime.mcrx(PI / 4, {0, 1}, 2);
+  checkEquivalence(qc, qcPrime, {0, 1, 2});
+}
+
+TEST_F(ZXFunctionalityTest, MCRx3) {
+  qc = qc::QuantumComputation(4);
+  qc.mcrx(PI / 4, {0, 1, 2}, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcrx(PI / 4, {0, 1, 2}, 3);
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCPhase) {
+  qc = qc::QuantumComputation(4);
+  qc.mcp(PI / 2, {0, 1, 2}, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcp(PI / 2, {0, 1, 2}, 3);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCS) {
+  qc = qc::QuantumComputation(4);
+  qc.mcs({0, 1, 2}, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mcs({0, 1, 2}, 3);
 
-  auto dPrime = FunctionalityConstruction::buildFunctionality(&qcPrime);
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
+
+TEST_F(ZXFunctionalityTest, MCT) {
+  qc = qc::QuantumComputation(4);
+  qc.mct({0, 1, 2}, 3);
+
+  auto qcPrime = qc::QuantumComputation(4);
+  qcPrime.mct({0, 1, 2}, 3);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2, 3});
+}
 
-  d.concat(dPrime.invert());
+TEST_F(ZXFunctionalityTest, MCPhase2) {
+  qc = qc::QuantumComputation(3);
+  qc.mcp(PI / 2, {0, 1}, 2);
+
+  auto qcPrime = qc::QuantumComputation(3);
+  qcPrime.mcp(PI / 2, {0, 1}, 2);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2});
+}
+
+TEST_F(ZXFunctionalityTest, MCS2) {
+  qc = qc::QuantumComputation(3);
+  qc.mcs({0, 1}, 2);
+
+  auto qcPrime = qc::QuantumComputation(3);
+  qcPrime.mcs({0, 1}, 2);
+
+  checkEquivalence(qc, qcPrime, {0, 1, 2});
+}
+
+TEST_F(ZXFunctionalityTest, MCT2) {
+  qc = qc::QuantumComputation(3);
+  qc.mct({0, 1}, 2);
 
-  fullReduce(d);
+  auto qcPrime = qc::QuantumComputation(3);
+  qcPrime.mct({0, 1}, 2);
 
-  EXPECT_TRUE(d.isIdentity());
-  EXPECT_TRUE(d.globalPhaseIsZero());
-  EXPECT_TRUE(d.connected(d.getInput(0), d.getOutput(0)));
+  checkEquivalence(qc, qcPrime, {0, 1, 2});
 }
 } // namespace zx