✨ Add native gate synthesis MLIR pass and decomposition stack

CHANGELOG.md

@@ -12,6 +12,9 @@ with the exception that minor releases may include breaking changes.
 
 ### Added
 
+- ✨ Add a `fuse-two-qubit-unitary-runs` pass
+  for fusing compile-time two-qubit unitary windows via Weyl/KAK resynthesis
+  ([#1655]) ([**@simon1hofmann**])
 - ✨ Add a `fuse-single-qubit-unitary-runs` pass
   for fusing compile-time single-qubit unitary runs via Euler resynthesis
   ([#1672]) ([**@simon1hofmann**], [**@burgholzer**])
@@ -631,6 +634,7 @@ changelogs._
 [#1664]: https://github.com/munich-quantum-toolkit/core/pull/1664
 [#1662]: https://github.com/munich-quantum-toolkit/core/pull/1662
 [#1660]: https://github.com/munich-quantum-toolkit/core/pull/1660
+[#1655]: https://github.com/munich-quantum-toolkit/core/pull/1655
 [#1652]: https://github.com/munich-quantum-toolkit/core/pull/1652
 [#1638]: https://github.com/munich-quantum-toolkit/core/pull/1638
 [#1637]: https://github.com/munich-quantum-toolkit/core/pull/1637

mlir/include/mlir/Compiler/CompilerPipeline.h

@@ -49,6 +49,18 @@ struct QuantumCompilerConfig {
 
   /// Enable Hadamard lifting
   bool enableHadamardLifting = false;
+
+  /// Comma-separated native gate menu. Recognised tokens: `u`, `x`, `sx`,
+  /// `rz` (or `p`), `rx`, `ry`, `r`, `cx`, `cz`, `rzz`.
+  /// Illustrative menus (use `cx` or `cz` as the entangler, or
+  /// both):
+  /// - `"x,sx,rz,cx"` / `"x,sx,rz,cz"` — IBM basic (no fractional 2q)
+  /// - `"x,sx,rz,rx,rzz,cx"` / `"...,cz"` — IBM fractional
+  /// - `"u,cx"` / `"u,cz"` — generic single-qubit U3 + CX/CZ
+  /// - `"r,cz"` — IQM-style default
+  /// - `"rx,rz,cx"`, `"rx,ry,cz"`, `"ry,rz,cx"` — supported RX/RY/RZ pairs plus
+  /// entangler
+  std::string nativeGates;
 };
 
 /**
@@ -84,7 +96,7 @@ struct CompilationRecord {
  * 2. QC cleanup pipeline
  * 3. QCO dialect (value semantics) - enables SSA-based optimizations
  * 4. QCO cleanup pipeline
- * 5. Quantum optimization passes
+ * 5. Optimization and native gate synthesis
  * 6. QCO cleanup pipeline
  * 7. QC dialect - converted back for backend lowering
  * 8. QC cleanup pipeline

mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.h

@@ -0,0 +1,231 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#pragma once
+
+#include "UnitaryMatrices.h"
+#include "WeylDecomposition.h"
+#include "mlir/Dialect/QCO/Utils/Matrix.h"
+
+#include <llvm/ADT/SmallVector.h>
+
+#include <array>
+#include <complex>
+#include <cstdint>
+#include <optional>
+#include <utility>
+
+namespace mlir::qco::decomposition {
+
+/// Intermediate single-qubit ``2×2`` unitaries produced while expanding a
+/// two-qubit basis decomposition.
+using TwoQubitLocalUnitaryList = llvm::SmallVector<Matrix2x2, 8>;
+
+/**
+ * Result of a two-qubit basis decomposition expressed as raw single-qubit
+ * factors interleaved with a fixed number of basis-gate (entangler) uses.
+ *
+ * The factors are stored in emission order. For `i` in `[0, numBasisUses)` the
+ * pair `(singleQubitFactors[2*i], singleQubitFactors[2*i + 1])` is applied to
+ * qubits `1` and `0` respectively, followed by one entangler. The final pair
+ * `(singleQubitFactors[2*numBasisUses], singleQubitFactors[2*numBasisUses+1])`
+ * is applied after the last entangler. The list therefore has length
+ * `2 * (numBasisUses + 1)`.
+ */
+struct TwoQubitNativeDecomposition {
+  /// Number of basis-gate (entangler) uses.
+  std::uint8_t numBasisUses = 0;
+  /// Single-qubit factors in emission order (see struct comment).
+  TwoQubitLocalUnitaryList singleQubitFactors;
+  /// Residual global phase (radians) not represented by factors/entanglers.
+  double globalPhase = 0.0;
+};
+
+/**
+ * Decomposer that must be initialized with a two-qubit basis gate that will
+ * be used to generate a circuit equivalent to a canonical gate (RXX+RYY+RZZ).
+ *
+ * @note Adapted from TwoQubitBasisDecomposer in the IBM Qiskit framework.
+ *       (C) Copyright IBM 2023
+ *
+ *       This code is licensed under the Apache License, Version 2.0. You may
+ *       obtain a copy of this license in the LICENSE.txt file in the root
+ *       directory of this source tree or at
+ *       https://www.apache.org/licenses/LICENSE-2.0.
+ *
+ *       Any modifications or derivative works of this code must retain this
+ *       copyright notice, and modified files need to carry a notice
+ *       indicating that they have been altered from the originals.
+ */
+class TwoQubitBasisDecomposer {
+public:
+  /**
+   * Create decomposer that allows two-qubit decompositions based on the
+   * specified entangler matrix.
+   * This entangler will appear between 0 and 3 times in each decomposition.
+   * The 4x4 matrix must be in MQT operand order (qubit 0 = MSB).
+   */
+  [[nodiscard]] static TwoQubitBasisDecomposer
+  create(const Matrix4x4& basisMatrix, double basisFidelity);
+
+  /**
+   * Perform decomposition using the basis gate of this decomposer.
+   *
+   * @param targetDecomposition Prepared Weyl decomposition of unitary matrix
+   *                            to be decomposed.
+   * @param numBasisGateUses Force use of given number of basis gates. When
+   *                         unset, the optimal count is selected from the
+   *                         Hilbert-Schmidt traces.
+   * @return The single-qubit factors and entangler count, or `std::nullopt`
+   *         when more than one basis gate would be required but the basis gate
+   *         is not super-controlled.
+   */
+  [[nodiscard]] std::optional<TwoQubitNativeDecomposition> twoQubitDecompose(
+      const decomposition::TwoQubitWeylDecomposition& targetDecomposition,
+      std::optional<std::uint8_t> numBasisGateUses) const;
+
+protected:
+  // NOLINTBEGIN(modernize-pass-by-value)
+  /**
+   * Constructs decomposer instance.
+   */
+  TwoQubitBasisDecomposer(
+      double basisFidelity,
+      const decomposition::TwoQubitWeylDecomposition& basisDecomposer,
+      bool isSuperControlled, const Matrix2x2& u0l, const Matrix2x2& u0r,
+      const Matrix2x2& u1l, const Matrix2x2& u1ra, const Matrix2x2& u1rb,
+      const Matrix2x2& u2la, const Matrix2x2& u2lb, const Matrix2x2& u2ra,
+      const Matrix2x2& u2rb, const Matrix2x2& u3l, const Matrix2x2& u3r,
+      const Matrix2x2& q0l, const Matrix2x2& q0r, const Matrix2x2& q1la,
+      const Matrix2x2& q1lb, const Matrix2x2& q1ra, const Matrix2x2& q1rb,
+      const Matrix2x2& q2l, const Matrix2x2& q2r)
+      : basisFidelity{basisFidelity}, basisDecomposer{basisDecomposer},
+        isSuperControlled{isSuperControlled}, u0l{u0l}, u0r{u0r}, u1l{u1l},
+        u1ra{u1ra}, u1rb{u1rb}, u2la{u2la}, u2lb{u2lb}, u2ra{u2ra}, u2rb{u2rb},
+        u3l{u3l}, u3r{u3r}, q0l{q0l}, q0r{q0r}, q1la{q1la}, q1lb{q1lb},
+        q1ra{q1ra}, q1rb{q1rb}, q2l{q2l}, q2r{q2r} {}
+  // NOLINTEND(modernize-pass-by-value)
+
+  /**
+   * Calculate decompositions when no basis gate is required.
+   *
+   * Decompose target :math:`\sim U_d(x, y, z)` with 0 uses of the
+   * basis gate. Result :math:`U_r` has trace:
+   *
+   * .. math::
+   *
+   *     \Big\vert\text{Tr}(U_r\cdot U_\text{target}^{\dag})\Big\vert =
+   *     4\Big\vert (\cos(x)\cos(y)\cos(z)+ j \sin(x)\sin(y)\sin(z)\Big\vert
+   *
+   * which is optimal for all targets and bases
+   */
+  [[nodiscard]] static TwoQubitLocalUnitaryList
+  decomp0(const decomposition::TwoQubitWeylDecomposition& target);
+
+  /**
+   * Calculate decompositions when one basis gate is required.
+   *
+   * Decompose target :math:`\sim U_d(x, y, z)` with 1 use of the
+   * basis gate :math:`\sim U_d(a, b, c)`. Result :math:`U_r` has trace:
+   *
+   * .. math::
+   *
+   *     \Big\vert\text{Tr}(U_r \cdot U_\text{target}^{\dag})\Big\vert =
+   *     4\Big\vert \cos(x-a)\cos(y-b)\cos(z-c) + j
+   *     \sin(x-a)\sin(y-b)\sin(z-c)\Big\vert
+   *
+   * which is optimal for all targets and bases with ``z==0`` or ``c==0``.
+   */
+  [[nodiscard]] TwoQubitLocalUnitaryList
+  decomp1(const decomposition::TwoQubitWeylDecomposition& target) const;
+
+  /**
+   * Calculate decompositions when two basis gates are required.
+   *
+   * Decompose target :math:`\sim U_d(x, y, z)` with 2 uses of the
+   * basis gate.
+   *
+   * For supercontrolled basis :math:`\sim U_d(\pi/4, b, 0)`, all b, result
+   * :math:`U_r` has trace
+   *
+   * .. math::
+   *
+   *     \Big\vert\text{Tr}(U_r \cdot U_\text{target}^\dag) \Big\vert =
+   * 4\cos(z)
+   *
+   * which is the optimal approximation for basis of CNOT-class
+   * :math:`\sim U_d(\pi/4, 0, 0)` or DCNOT-class
+   * :math:`\sim U_d(\pi/4, \pi/4, 0)` and any target. It may be sub-optimal
+   * for :math:`b \neq 0` (i.e. there exists an exact decomposition for any
+   * target using :math:`B \sim U_d(\pi/4, \pi/8, 0)`, but it may not be this
+   * decomposition). This is an exact decomposition for supercontrolled basis
+   * and target :math:`\sim U_d(x, y, 0)`. No guarantees for
+   * non-supercontrolled basis.
+   */
+  [[nodiscard]] TwoQubitLocalUnitaryList decomp2Supercontrolled(
+      const decomposition::TwoQubitWeylDecomposition& target) const;
+
+  /**
+   * Calculate decompositions when three basis gates are required.
+   *
+   * Decompose target with 3 uses of the basis.
+   *
+   * This is an exact decomposition for supercontrolled basis
+   * :math:`\sim U_d(\pi/4, b, 0)`, all b, and any target. No guarantees for
+   * non-supercontrolled basis.
+   */
+  [[nodiscard]] TwoQubitLocalUnitaryList decomp3Supercontrolled(
+      const decomposition::TwoQubitWeylDecomposition& target) const;
+
+  /**
+   * Calculate traces for a combination of the parameters of the canonical
+   * gates of the target and basis decompositions.
+   * This can be used to determine the smallest number of basis gates that are
+   * necessary to construct an equivalent to the canonical gate.
+   */
+  [[nodiscard]] std::array<std::complex<double>, 4>
+  traces(const decomposition::TwoQubitWeylDecomposition& target) const;
+
+  [[nodiscard]] static bool relativeEq(double lhs, double rhs, double epsilon,
+                                       double maxRelative);
+
+private:
+  // fidelity with which the basis gate decomposition has been calculated
+  double basisFidelity;
+  // cached decomposition for basis gate
+  decomposition::TwoQubitWeylDecomposition basisDecomposer;
+  // true if basis gate is super-controlled
+  bool isSuperControlled;
+
+  // pre-built components for decomposition with 3 basis gates
+  Matrix2x2 u0l;
+  Matrix2x2 u0r;
+  Matrix2x2 u1l;
+  Matrix2x2 u1ra;
+  Matrix2x2 u1rb;
+  Matrix2x2 u2la;
+  Matrix2x2 u2lb;
+  Matrix2x2 u2ra;
+  Matrix2x2 u2rb;
+  Matrix2x2 u3l;
+  Matrix2x2 u3r;
+
+  // pre-built components for decomposition with 2 basis gates
+  Matrix2x2 q0l;
+  Matrix2x2 q0r;
+  Matrix2x2 q1la;
+  Matrix2x2 q1lb;
+  Matrix2x2 q1ra;
+  Matrix2x2 q1rb;
+  Matrix2x2 q2l;
+  Matrix2x2 q2r;
+};
+
+} // namespace mlir::qco::decomposition

mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Euler.h

@@ -32,6 +32,7 @@ enum class EulerBasis : std::uint8_t {
   XYX = 3,  ///< `RX(phi) * RY(theta) * RX(lambda)`.
   U = 4,    ///< `U(theta, phi, lambda)`.
   ZSXX = 5, ///< `RZ` / `SX` / `X` synthesis via ZYZ decomposition.
+  R = 6,    ///< `R(.,0) * R(.,pi/2) * R(.,0)` (XYX with `Rx`/`Ry` as `R`).
 };
 
 /**
@@ -42,6 +43,31 @@ enum class EulerBasis : std::uint8_t {
  */
 [[nodiscard]] std::optional<EulerBasis> parseEulerBasis(StringRef basis);
 
+/**
+ * @brief Euler angles `(theta, phi, lambda)` and global phase for a 2x2
+ * unitary.
+ *
+ * The decomposition obeys `matrix == e^{i*phase} * K(phi) * A(theta) *
+ * K(lambda)` where `(K, A)` are the rotation axes of the chosen @ref
+ * EulerBasis.
+ */
+struct EulerAngles {
+  double theta = 0.0;  ///< Middle rotation angle.
+  double phi = 0.0;    ///< First outer rotation angle.
+  double lambda = 0.0; ///< Second outer rotation angle.
+  double phase = 0.0;  ///< Global phase in radians.
+};
+
+/**
+ * @brief Extracts `(theta, phi, lambda, phase)` of @p matrix in @p basis.
+ *
+ * @param matrix The single-qubit unitary to decompose.
+ * @param basis The target Euler basis.
+ * @return The extracted Euler angles and global phase.
+ */
+[[nodiscard]] EulerAngles anglesFromUnitary(const Matrix2x2& matrix,
+                                            EulerBasis basis);
+
 /**
  * @brief Synthesizes a composed single-qubit unitary as gates in @p basis.
  *

mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Helpers.h

@@ -0,0 +1,43 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#pragma once
+
+#include "mlir/Dialect/QCO/Utils/Matrix.h"
+
+#include <complex>
+
+/// Numeric helpers used by the decomposition passes.
+
+namespace mlir::qco::helpers {
+
+/// Check whether `matrix` is unitary within `tolerance` (i.e. `M^H M` is
+/// approximately the identity).
+[[nodiscard]] bool isUnitaryMatrix(const Matrix2x2& matrix,
+                                   double tolerance = 1e-12);
+
+/**
+ * Euclidean remainder of a modulo b.
+ * The returned value is never negative.
+ */
+[[nodiscard]] double remEuclid(double a, double b);
+
+/**
+ * Convert a two-qubit trace overlap into the average gate fidelity metric used
+ * by the decomposition cost code.
+ */
+[[nodiscard]] double traceToFidelity(const std::complex<double>& x);
+
+/**
+ * Return the scalar `e^(i * globalPhase)` factor for a stored global phase.
+ */
+[[nodiscard]] std::complex<double> globalPhaseFactor(double globalPhase);
+
+} // namespace mlir::qco::helpers