✨ Add two-qubit Weyl (KAK) decomposition library

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

@@ -42,6 +42,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/Weyl.h

@@ -0,0 +1,169 @@
+/*
+ * 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 <mlir/Support/LLVM.h>
+
+#include <array>
+#include <complex>
+#include <cstdint>
+#include <optional>
+
+namespace mlir::qco::decomposition {
+
+/** Numeric tolerance for Weyl internal matrix checks. */
+inline constexpr double WEYL_TOLERANCE = 100 * MATRIX_TOLERANCE;
+
+/**
+ * @brief Weyl decomposition of a 2-qubit unitary.
+ *
+ * A 4x4 unitary is factored as
+ * `(K1l ⊗ K1r) · U_canon(a,b,c) · (K2l ⊗ K2r)` up to global phase, where
+ * `U_canon(a,b,c) = RXX(-2a) · RYY(-2b) · RZZ(-2c)`.
+ *
+ * @note Adapted from TwoQubitWeylDecomposition 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 TwoQubitWeylDecomposition {
+public:
+  [[nodiscard]] static TwoQubitWeylDecomposition
+  create(const Matrix4x4& unitaryMatrix, std::optional<double> fidelity);
+
+  [[nodiscard]] Matrix4x4 getCanonicalMatrix() const {
+    return getCanonicalMatrix(a_, b_, c_);
+  }
+
+  [[nodiscard]] double a() const { return a_; }
+  [[nodiscard]] double b() const { return b_; }
+  [[nodiscard]] double c() const { return c_; }
+  [[nodiscard]] double globalPhase() const { return globalPhase_; }
+
+  /** @brief Single-qubit factor on qubit 0 after the canonical gate. */
+  [[nodiscard]] const Matrix2x2& k1l() const { return k1l_; }
+  /** @brief Single-qubit factor on qubit 0 before the canonical gate. */
+  [[nodiscard]] const Matrix2x2& k2l() const { return k2l_; }
+  /** @brief Single-qubit factor on qubit 1 after the canonical gate. */
+  [[nodiscard]] const Matrix2x2& k1r() const { return k1r_; }
+  /** @brief Single-qubit factor on qubit 1 before the canonical gate. */
+  [[nodiscard]] const Matrix2x2& k2r() const { return k2r_; }
+
+  [[nodiscard]] static Matrix4x4 getCanonicalMatrix(double a, double b,
+                                                    double c);
+
+private:
+  bool applySpecialization(const std::optional<double>& requestedFidelity);
+
+  double a_{};
+  double b_{};
+  double c_{};
+  double globalPhase_{};
+  Matrix2x2 k1l_;
+  Matrix2x2 k2l_;
+  Matrix2x2 k1r_;
+  Matrix2x2 k2r_;
+};
+
+/**
+ * @brief Native two-qubit synthesis result.
+ *
+ * Emission order: for each `i`, apply `singleQubitFactors[2*i+1]` on q0 and
+ * `singleQubitFactors[2*i]` on q1, then the basis gate (except after the last
+ * pair).
+ */
+struct TwoQubitNativeDecomposition {
+  std::uint8_t numBasisUses = 0;
+  SmallVector<Matrix2x2, 8> singleQubitFactors;
+  double globalPhase = 0.0;
+};
+
+/**
+ * @brief Decomposer for a fixed two-qubit basis gate (e.g. CX/CZ).
+ *
+ * @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:
+  [[nodiscard]] static TwoQubitBasisDecomposer
+  create(const Matrix4x4& basisMatrix, double basisFidelity);
+
+  [[nodiscard]] std::optional<TwoQubitNativeDecomposition>
+  twoQubitDecompose(const TwoQubitWeylDecomposition& targetDecomposition,
+                    std::optional<std::uint8_t> numBasisGateUses) const;
+
+private:
+  struct SmbPrecomputed {
+    Matrix2x2 u0l;
+    Matrix2x2 u0r;
+    Matrix2x2 u1l;
+    Matrix2x2 u1ra;
+    Matrix2x2 u1rb;
+    Matrix2x2 u2la;
+    Matrix2x2 u2lb;
+    Matrix2x2 u2ra;
+    Matrix2x2 u2rb;
+    Matrix2x2 u3l;
+    Matrix2x2 u3r;
+    Matrix2x2 q0l;
+    Matrix2x2 q0r;
+    Matrix2x2 q1la;
+    Matrix2x2 q1lb;
+    Matrix2x2 q1ra;
+    Matrix2x2 q1rb;
+    Matrix2x2 q2l;
+    Matrix2x2 q2r;
+  };
+
+  [[nodiscard]] static SmallVector<Matrix2x2, 8>
+  decomp0(const TwoQubitWeylDecomposition& target);
+  [[nodiscard]] SmallVector<Matrix2x2, 8>
+  decomp1(const TwoQubitWeylDecomposition& target) const;
+  [[nodiscard]] SmallVector<Matrix2x2, 8>
+  decomp2Supercontrolled(const TwoQubitWeylDecomposition& target) const;
+  [[nodiscard]] SmallVector<Matrix2x2, 8>
+  decomp3Supercontrolled(const TwoQubitWeylDecomposition& target) const;
+  [[nodiscard]] std::array<std::complex<double>, 4>
+  traces(const TwoQubitWeylDecomposition& target) const;
+
+  double basisFidelity{};
+  TwoQubitWeylDecomposition basisDecomposer;
+  bool isSuperControlled{};
+  SmbPrecomputed smb{};
+};
+
+/** @brief Weyl-decomposes @p target using @p basisMatrix as entangler. */
+[[nodiscard]] std::optional<TwoQubitNativeDecomposition>
+decomposeTwoQubitWithBasis(
+    const Matrix4x4& target, const Matrix4x4& basisMatrix,
+    double basisFidelity = 1.0,
+    std::optional<std::uint8_t> numBasisUses = std::nullopt);
+
+} // namespace mlir::qco::decomposition

mlir/include/mlir/Dialect/QCO/Utils/Matrix.h

@@ -795,4 +795,31 @@ struct SymmetricEigen4 {
   Matrix4x4 eigenvectors{};
 };
 
+/** @brief `RX(theta)` unitary matrix. */
+[[nodiscard]] Matrix2x2 rxMatrix(double theta);
+
+/** @brief `RY(theta)` unitary matrix. */
+[[nodiscard]] Matrix2x2 ryMatrix(double theta);
+
+/** @brief `RZ(theta)` unitary matrix. */
+[[nodiscard]] Matrix2x2 rzMatrix(double theta);
+
+/** @brief `i` times Pauli-X. */
+[[nodiscard]] const Matrix2x2& iPauliX();
+
+/** @brief `i` times Pauli-Y. */
+[[nodiscard]] const Matrix2x2& iPauliY();
+
+/** @brief `i` times Pauli-Z. */
+[[nodiscard]] const Matrix2x2& iPauliZ();
+
+/** @brief `RXX(theta)` unitary matrix. */
+[[nodiscard]] Matrix4x4 rxxMatrix(double theta);
+
+/** @brief `RYY(theta)` unitary matrix. */
+[[nodiscard]] Matrix4x4 ryyMatrix(double theta);
+
+/** @brief `RZZ(theta)` unitary matrix. */
+[[nodiscard]] Matrix4x4 rzzMatrix(double theta);
+
 } // namespace mlir::qco

mlir/lib/Dialect/QCO/Transforms/Decomposition/Euler.cpp

@@ -102,21 +102,6 @@ static void emitGPhaseIfNeeded(OpBuilder& builder, Location loc, double phase) {
 // Euler decomposition (angles)
 //===----------------------------------------------------------------------===//
 
-/**
- * @brief Euler angles `(theta, phi, lambda)` and global phase for a 2x2
- * unitary.
- */
-namespace {
-
-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.
-};
-
-} // namespace
-
 /**
  * @brief Z-Y-Z Euler angles and global phase for a 2x2 unitary.
  *
@@ -197,15 +182,7 @@ struct EulerAngles {
           .phase = phase - (0.5 * (phi + lambda))};
 }
 
-/**
- * @brief Extracts `(theta, phi, lambda, phase)` for all Euler bases.
- *
- * @param matrix The single-qubit unitary to decompose.
- * @param basis The target Euler basis.
- * @return The extracted Euler angles and global phase.
- */
-[[nodiscard]] static EulerAngles anglesFromUnitary(const Matrix2x2& matrix,
-                                                   const EulerBasis basis) {
+EulerAngles anglesFromUnitary(const Matrix2x2& matrix, const EulerBasis basis) {
   switch (basis) {
   case EulerBasis::ZYZ:
   case EulerBasis::ZSXX: