diff --git a/CHANGELOG.md b/CHANGELOG.md index 614fe3c174..2a2888002c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,7 +42,7 @@ with the exception that minor releases may include breaking changes. [#1569], [#1570], [#1572], [#1573], [#1580], [#1602], [#1620], [#1623], [#1624], [#1626], [#1627], [#1635], [#1638], [#1673], [#1675], [#1700], [#1710], [#1717], [#1728], [#1730], [#1749], [#1751], [#1762], [#1765], - [#1774], [#1780], [#1781], [#1782], [#1787], [#1802]) + [#1774], [#1780], [#1781], [#1782], [#1787], [#1802], [#1803]) ([**@burgholzer**], [**@denialhaag**], [**@taminob**], [**@DRovara**], [**@li-mingbao**], [**@Ectras**], [**@MatthiasReumann**], [**@simon1hofmann**]) @@ -598,6 +598,7 @@ changelogs._ +[#1803]: https://github.com/munich-quantum-toolkit/core/pull/1803 [#1802]: https://github.com/munich-quantum-toolkit/core/pull/1802 [#1787]: https://github.com/munich-quantum-toolkit/core/pull/1787 [#1782]: https://github.com/munich-quantum-toolkit/core/pull/1782 diff --git a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Euler.h b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Euler.h index 8fb0018b28..7cea37b3a6 100644 --- a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Euler.h +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Euler.h @@ -42,6 +42,31 @@ enum class EulerBasis : std::uint8_t { */ [[nodiscard]] std::optional 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. * diff --git a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h new file mode 100644 index 0000000000..d5db42398d --- /dev/null +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h @@ -0,0 +1,214 @@ +/* + * 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 + +#include +#include +#include +#include + +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 2026 + * + * 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 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& requestedFidelity); + + void finalizeSpecializationPhase(bool flippedFromOriginal, + double preSpecializationA, + double preSpecializationB, + double preSpecializationC, + const std::optional& fidelity); + + 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 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 2026 + * + * 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: + /** + * @brief Precomputes basis-gate data for repeated target decompositions. + * + * Creation performs a full Weyl decomposition of @p basisMatrix and + * precomputes the single-qubit templates used by @ref twoQubitDecompose. + * Reuse one instance for many targets that share the same entangler + * (e.g. CX) via @ref decomposeTarget. + */ + [[nodiscard]] static TwoQubitBasisDecomposer + create(const Matrix4x4& basisMatrix, double basisFidelity); + + /** + * @brief Decomposes a target Weyl decomposition into single-qubit factors and + * basis-gate uses. + * + * @param targetDecomposition Weyl decomposition of the target unitary. + * @param numBasisGateUses Requested number of basis-gate applications in + * `{0, 1, 2, 3}`. Pass `std::nullopt` to pick the count that + * maximizes `traceToFidelity(trace[i]) * basisFidelity^i` over + * `i ∈ {0, 1, 2, 3}`. + * @return A native decomposition on success, or `std::nullopt` when the + * requested count is unsupported (e.g. more than one basis gate for a + * non-super-controlled basis, or a value outside `{0, 1, 2, 3}`). + */ + [[nodiscard]] std::optional + twoQubitDecompose(const TwoQubitWeylDecomposition& targetDecomposition, + std::optional numBasisGateUses) const; + + /** + * @brief Decomposes @p targetUnitary using this cached basis decomposer. + * + * Only the target undergoes Weyl decomposition; basis precomputation from + * @ref create is reused. + */ + [[nodiscard]] std::optional decomposeTarget( + const Matrix4x4& targetUnitary, + std::optional numBasisGateUses = std::nullopt) 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 + decomp0(const TwoQubitWeylDecomposition& target); + [[nodiscard]] SmallVector + decomp1(const TwoQubitWeylDecomposition& target) const; + [[nodiscard]] SmallVector + decomp2Supercontrolled(const TwoQubitWeylDecomposition& target) const; + [[nodiscard]] SmallVector + decomp3Supercontrolled(const TwoQubitWeylDecomposition& target) const; + [[nodiscard]] std::array, 4> + traces(const TwoQubitWeylDecomposition& target) const; + + TwoQubitBasisDecomposer() = default; + + double basisFidelity{}; + TwoQubitWeylDecomposition basisWeyl; + bool isSuperControlled{}; + SmbPrecomputed smb{}; +}; + +/** + * @brief Convenience wrapper that builds a fresh basis decomposer per call. + * + * For a fixed basis gate decomposed many times, prefer caching + * `TwoQubitBasisDecomposer::create(basisMatrix, basisFidelity)` and calling + * `TwoQubitBasisDecomposer::decomposeTarget` for each target. + */ +[[nodiscard]] std::optional +decomposeTwoQubitWithBasis( + const Matrix4x4& target, const Matrix4x4& basisMatrix, + double basisFidelity = 1.0, + std::optional numBasisUses = std::nullopt); + +} // namespace mlir::qco::decomposition diff --git a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h index cdea703ed4..d9d46a714e 100644 --- a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h +++ b/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 diff --git a/mlir/lib/Dialect/QCO/Transforms/CMakeLists.txt b/mlir/lib/Dialect/QCO/Transforms/CMakeLists.txt index 3564b55dc5..77f594db17 100644 --- a/mlir/lib/Dialect/QCO/Transforms/CMakeLists.txt +++ b/mlir/lib/Dialect/QCO/Transforms/CMakeLists.txt @@ -22,6 +22,8 @@ add_mlir_library( DEPENDS MLIRQCOTransformsIncGen) +mqt_mlir_target_use_project_options(MLIRQCOTransforms) + # collect header files file(GLOB_RECURSE PASSES_HEADERS_SOURCE ${MQT_MLIR_SOURCE_INCLUDE_DIR}/mlir/Dialect/QCO/Transforms/*.h) diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp new file mode 100644 index 0000000000..814702544d --- /dev/null +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp @@ -0,0 +1,312 @@ +/* + * 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 + */ + +#include "mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h" +#include "mlir/Dialect/QCO/Utils/Matrix.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mlir::qco::decomposition { + +using namespace std::complex_literals; + +static constexpr double DEFAULT_WEYL_FIDELITY = 1.0 - 1e-12; +static constexpr double PI = std::numbers::pi; +static constexpr double PI_OVER_4 = PI / 4.0; +static constexpr double INV_SQRT2 = 1.0 / std::numbers::sqrt2; + +static const Matrix2x2 K12_R_ARR = Matrix2x2::fromElements( + 1i * INV_SQRT2, INV_SQRT2, -INV_SQRT2, -1i * INV_SQRT2); +static const Matrix2x2 K12_L_ARR = + Matrix2x2::fromElements(Complex{0.5, 0.5}, Complex{0.5, 0.5}, + Complex{-0.5, 0.5}, Complex{0.5, -0.5}); +static const Matrix2x2 K22_L = + Matrix2x2::fromElements(INV_SQRT2, -INV_SQRT2, INV_SQRT2, INV_SQRT2); +static const Matrix2x2 K22_R = Matrix2x2::fromElements(0, 1, -1, 0); + +static double remEuclid(const double a, const double b) { + if (b == 0.0) { + llvm::reportFatalInternalError("remEuclid expects non-zero divisor"); + } + const auto r = std::fmod(a, b); + return (r < 0.0) ? r + std::abs(b) : r; +} + +static bool relativeEq(double lhs, double rhs, double epsilon, + double maxRelative) { + if (lhs == rhs) { + return true; + } + if (std::isinf(lhs) || std::isinf(rhs)) { + return false; + } + const auto absDiff = std::abs(lhs - rhs); + if (absDiff <= epsilon) { + return true; + } + const auto absLhs = std::abs(lhs); + const auto absRhs = std::abs(rhs); + if (absRhs > absLhs) { + return absDiff <= absRhs * maxRelative; + } + return absDiff <= absLhs * maxRelative; +} + +static double traceToFidelity(const Complex& trace) { + const auto traceAbs = std::abs(trace); + return (4.0 + (traceAbs * traceAbs)) / 20.0; +} + +//===----------------------------------------------------------------------===// +// TwoQubitBasisDecomposer +//===----------------------------------------------------------------------===// + +TwoQubitBasisDecomposer +TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, + double basisFidelity) { + if (!std::isfinite(basisFidelity) || basisFidelity < 0.0 || + basisFidelity > 1.0) { + llvm::reportFatalInternalError(llvm::formatv( + "TwoQubitBasisDecomposer: basisFidelity must be finite and in [0, 1] " + "(got {0})", + basisFidelity)); + } + + const auto basisWeyl = + TwoQubitWeylDecomposition::create(basisMatrix, DEFAULT_WEYL_FIDELITY); + const auto isSuperControlled = + relativeEq(basisWeyl.a(), PI_OVER_4, 1e-13, 1e-09) && + relativeEq(basisWeyl.c(), 0.0, 1e-13, 1e-09); + + const auto b = basisWeyl.b(); + Complex temp{0.5, -0.5}; + const Matrix2x2 k11l = Matrix2x2::fromElements( + temp * (-1i * std::exp(-1i * b)), temp * std::exp(-1i * b), + temp * (-1i * std::exp(1i * b)), temp * -std::exp(1i * b)); + const Matrix2x2 k11r = Matrix2x2::fromElements( + INV_SQRT2 * (1i * std::exp(-1i * b)), INV_SQRT2 * -std::exp(-1i * b), + INV_SQRT2 * std::exp(1i * b), INV_SQRT2 * (-1i * std::exp(1i * b))); + const Matrix2x2 k32lK21l = Matrix2x2::fromElements( + INV_SQRT2 * Complex{1., std::cos(2. * b)}, + INV_SQRT2 * (1i * std::sin(2. * b)), INV_SQRT2 * (1i * std::sin(2. * b)), + INV_SQRT2 * Complex{1., -std::cos(2. * b)}); + temp = Complex{0.5, 0.5}; + const Matrix2x2 k21r = Matrix2x2::fromElements( + temp * (-1i * std::exp(-2i * b)), temp * std::exp(-2i * b), + temp * (1i * std::exp(2i * b)), temp * std::exp(2i * b)); + const Matrix2x2 k31l = Matrix2x2::fromElements( + INV_SQRT2 * std::exp(-1i * b), INV_SQRT2 * std::exp(-1i * b), + INV_SQRT2 * -std::exp(1i * b), INV_SQRT2 * std::exp(1i * b)); + const Matrix2x2 k31r = Matrix2x2::fromElements(1i * std::exp(1i * b), 0, 0, + -1i * std::exp(-1i * b)); + const Matrix2x2 k32r = Matrix2x2::fromElements( + temp * std::exp(1i * b), temp * -std::exp(-1i * b), + temp * (-1i * std::exp(1i * b)), temp * (-1i * std::exp(-1i * b))); + + const auto k1lDagger = basisWeyl.k1l().adjoint(); + const auto k1rDagger = basisWeyl.k1r().adjoint(); + const auto k2lDagger = basisWeyl.k2l().adjoint(); + const auto k2rDagger = basisWeyl.k2r().adjoint(); + + TwoQubitBasisDecomposer decomposer; + decomposer.basisFidelity = basisFidelity; + decomposer.basisWeyl = basisWeyl; + decomposer.isSuperControlled = isSuperControlled; + decomposer.smb = TwoQubitBasisDecomposer::SmbPrecomputed{ + .u0l = k31l * k1lDagger, + .u0r = k31r * k1rDagger, + .u1l = k2lDagger * k32lK21l * k1lDagger, + .u1ra = k2rDagger * k32r, + .u1rb = k21r * k1rDagger, + .u2la = k2lDagger * K22_L, + .u2lb = k11l * k1lDagger, + .u2ra = k2rDagger * K22_R, + .u2rb = k11r * k1rDagger, + .u3l = k2lDagger * K12_L_ARR, + .u3r = k2rDagger * K12_R_ARR, + .q0l = K12_L_ARR.adjoint() * k1lDagger, + .q0r = K12_R_ARR.adjoint() * iPauliZ() * k1rDagger, + .q1la = k2lDagger * k11l.adjoint(), + .q1lb = k11l * k1lDagger, + .q1ra = k2rDagger * iPauliZ() * k11r.adjoint(), + .q1rb = k11r * k1rDagger, + .q2l = k2lDagger * K12_L_ARR, + .q2r = k2rDagger * K12_R_ARR, + }; + return decomposer; +} + +std::optional +TwoQubitBasisDecomposer::decomposeTarget( + const Matrix4x4& targetUnitary, + const std::optional numBasisGateUses) const { + const auto targetWeyl = + TwoQubitWeylDecomposition::create(targetUnitary, DEFAULT_WEYL_FIDELITY); + return twoQubitDecompose(targetWeyl, numBasisGateUses); +} + +std::optional +TwoQubitBasisDecomposer::twoQubitDecompose( + const TwoQubitWeylDecomposition& targetDecomposition, + std::optional numBasisGateUses) const { + const auto traceValues = traces(targetDecomposition); + auto getDefaultNbasis = [&]() -> std::uint8_t { + auto bestValue = std::numeric_limits::lowest(); + auto bestIndex = -1; + double fidelityPower = 1.0; + for (int i = 0; std::cmp_less(i, traceValues.size()); ++i) { + const auto value = traceToFidelity(traceValues[i]) * fidelityPower; + fidelityPower *= basisFidelity; + if (std::isnan(value)) { + continue; + } + if (value > bestValue) { + bestIndex = i; + bestValue = value; + } + } + if (bestIndex < 0) { + llvm::reportFatalInternalError("Unable to select basis-gate count: all " + "candidate fidelities are NaN"); + } + return static_cast(bestIndex); + }; + + const auto bestNbasis = numBasisGateUses.value_or(getDefaultNbasis()); + if (bestNbasis > 1 && !isSuperControlled) { + return std::nullopt; + } + + SmallVector factors; + switch (bestNbasis) { + case 0: + factors = decomp0(targetDecomposition); + break; + case 1: + factors = decomp1(targetDecomposition); + break; + case 2: + factors = decomp2Supercontrolled(targetDecomposition); + break; + case 3: + factors = decomp3Supercontrolled(targetDecomposition); + break; + default: + return std::nullopt; + } + + double globalPhase = targetDecomposition.globalPhase(); + globalPhase -= bestNbasis * basisWeyl.globalPhase(); + if (bestNbasis == 2) { + globalPhase += PI; + } + globalPhase = remEuclid(globalPhase, 2.0 * PI); + + return TwoQubitNativeDecomposition{ + .numBasisUses = bestNbasis, + .singleQubitFactors = std::move(factors), + .globalPhase = globalPhase, + }; +} + +SmallVector +TwoQubitBasisDecomposer::decomp0(const TwoQubitWeylDecomposition& target) { + return SmallVector{ + target.k1r() * target.k2r(), + target.k1l() * target.k2l(), + }; +} + +SmallVector TwoQubitBasisDecomposer::decomp1( + const TwoQubitWeylDecomposition& target) const { + return SmallVector{ + basisWeyl.k2r().adjoint() * target.k2r(), + basisWeyl.k2l().adjoint() * target.k2l(), + target.k1r() * basisWeyl.k1r().adjoint(), + target.k1l() * basisWeyl.k1l().adjoint(), + }; +} + +SmallVector TwoQubitBasisDecomposer::decomp2Supercontrolled( + const TwoQubitWeylDecomposition& target) const { + if (!isSuperControlled) { + llvm::reportFatalInternalError( + "Basis gate of TwoQubitBasisDecomposer is not super-controlled " + "- no guarantee for exact decomposition with two basis gates"); + } + return SmallVector{ + smb.q2r * target.k2r(), + smb.q2l * target.k2l(), + smb.q1ra * rzMatrix(2. * target.b()) * smb.q1rb, + smb.q1la * rzMatrix(-2. * target.a()) * smb.q1lb, + target.k1r() * smb.q0r, + target.k1l() * smb.q0l, + }; +} + +SmallVector TwoQubitBasisDecomposer::decomp3Supercontrolled( + const TwoQubitWeylDecomposition& target) const { + if (!isSuperControlled) { + llvm::reportFatalInternalError( + "Basis gate of TwoQubitBasisDecomposer is not super-controlled " + "- no guarantee for exact decomposition with three basis gates"); + } + return SmallVector{ + smb.u3r * target.k2r(), + smb.u3l * target.k2l(), + smb.u2ra * rzMatrix(2. * target.b()) * smb.u2rb, + smb.u2la * rzMatrix(-2. * target.a()) * smb.u2lb, + smb.u1ra * rzMatrix(-2. * target.c()) * smb.u1rb, + smb.u1l, + target.k1r() * smb.u0r, + target.k1l() * smb.u0l, + }; +} + +std::array, 4> +TwoQubitBasisDecomposer::traces(const TwoQubitWeylDecomposition& target) const { + return { + 4. * std::complex{std::cos(target.a()) * std::cos(target.b()) * + std::cos(target.c()), + std::sin(target.a()) * std::sin(target.b()) * + std::sin(target.c())}, + 4. * std::complex{std::cos(PI_OVER_4 - target.a()) * + std::cos(basisWeyl.b() - target.b()) * + std::cos(target.c()), + std::sin(PI_OVER_4 - target.a()) * + std::sin(basisWeyl.b() - target.b()) * + std::sin(target.c())}, + std::complex{4. * std::cos(target.c()), 0.}, + std::complex{4., 0.}, + }; +} + +std::optional +decomposeTwoQubitWithBasis(const Matrix4x4& target, + const Matrix4x4& basisMatrix, + const double basisFidelity, + const std::optional numBasisUses) { + const auto decomposer = + TwoQubitBasisDecomposer::create(basisMatrix, basisFidelity); + return decomposer.decomposeTarget(target, numBasisUses); +} + +} // namespace mlir::qco::decomposition diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Euler.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Euler.cpp index 085d493abf..077cdde0e1 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Euler.cpp +++ b/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: diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp new file mode 100644 index 0000000000..a857fdc9cc --- /dev/null +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -0,0 +1,635 @@ +/* + * 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 + */ + +#include "mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h" + +#include "mlir/Dialect/QCO/Transforms/Decomposition/Euler.h" +#include "mlir/Dialect/QCO/Utils/Matrix.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mlir::qco::decomposition { + +using namespace std::complex_literals; + +namespace { + +enum class Specialization : std::uint8_t { + General, + IdEquiv, + SWAPEquiv, + PartialSWAPEquiv, + PartialSWAPFlipEquiv, + ControlledEquiv, + MirrorControlledEquiv, + FSimaabEquiv, + FSimabbEquiv, + FSimabmbEquiv, +}; + +struct ChamberState { + std::array cs{}; + Matrix2x2 k1l; + Matrix2x2 k1r; + Matrix2x2 k2l; + Matrix2x2 k2r; + double globalPhase{}; + double a{}; + double b{}; + double c{}; +}; + +} // namespace + +static constexpr double DIAGONALIZATION_PRECISION = 1e-13; +static constexpr double PI = std::numbers::pi; +static constexpr double PI_OVER_4 = PI / 4.0; + +static double remEuclid(const double a, const double b) { + if (b == 0.0) { + llvm::reportFatalInternalError("remEuclid expects non-zero divisor"); + } + const auto r = std::fmod(a, b); + return (r < 0.0) ? r + std::abs(b) : r; +} + +/** Maps `|tr(U^dag V)|` to average two-qubit gate fidelity. */ +static double traceToFidelity(const Complex& trace) { + const auto traceAbs = std::abs(trace); + return (4.0 + (traceAbs * traceAbs)) / 20.0; +} + +static const Matrix4x4 MAGIC_BASIS_NON_NORMALIZED = Matrix4x4::fromElements( // + 1, 1i, 0, 0, // + 0, 0, 1i, 1, // + 0, 0, 1i, -1, // + 1, -1i, 0, 0); +static const Matrix4x4 MAGIC_BASIS_NON_NORMALIZED_DAGGER = + Matrix4x4::fromElements( // + 0.5, 0, 0, 0.5, // + Complex{0.0, -0.5}, 0, 0, Complex{0.0, 0.5}, // + 0, Complex{0.0, -0.5}, Complex{0.0, -0.5}, 0, // + 0, 0.5, -0.5, 0); + +static Matrix4x4 magicBasisTransform(const Matrix4x4& unitary, + bool outOfMagicBasis) { + if (outOfMagicBasis) { + return MAGIC_BASIS_NON_NORMALIZED_DAGGER * unitary * + MAGIC_BASIS_NON_NORMALIZED; + } + return MAGIC_BASIS_NON_NORMALIZED * unitary * + MAGIC_BASIS_NON_NORMALIZED_DAGGER; +} + +static double closestPartialSwap(double a, double b, double c) { + const auto m = (a + b + c) / 3.; + const auto [am, bm, cm] = std::array{a - m, b - m, c - m}; + const auto [ab, bc, ca] = std::array{a - b, b - c, c - a}; + return m + (am * bm * cm * (6. + (ab * ab) + (bc * bc) + (ca * ca)) / 18.); +} + +static std::pair> +diagonalizeComplexSymmetric(const Matrix4x4& m, double precision) { + auto state = std::mt19937{2023}; + std::normal_distribution dist; + + const auto mReal = m.realPart(); + const auto mImag = m.imagPart(); + + double bestErr = 1e300; + constexpr auto maxDiagonalizationAttempts = 100; + for (int i = 0; i < maxDiagonalizationAttempts; ++i) { + double randA{}; + double randB{}; + // Fixed perturbation coefficients for the first diagonalization attempt, + // carried over from Qiskit's two-qubit Weyl decomposition (legacy Python + // RNG values). The loop usually succeeds on this trial; fixing randA/randB + // keeps behavior deterministic while later attempts sample the + // distribution. + if (i == 0) { + randA = 1.2602066112249388; + randB = 0.22317849046722027; + } else { + randA = dist(state); + randB = dist(state); + } + std::array m2Real{}; + for (std::size_t k = 0; k < m2Real.size(); ++k) { + m2Real[k] = (randA * mReal[k]) + (randB * mImag[k]); + } + const Matrix4x4 p = Matrix4x4::symmetricEigen4(m2Real).eigenvectors; + const std::array d = (p.transpose() * m * p).diagonal(); + + const auto compare = p * Matrix4x4::fromDiagonal(d) * p.transpose(); + { + double err = 0.0; + for (std::size_t r = 0; r < 4; ++r) { + for (std::size_t cc = 0; cc < 4; ++cc) { + err = std::max(err, std::abs(compare(r, cc) - m(r, cc))); + } + } + bestErr = std::min(bestErr, err); + } + if (compare.isApprox(m, precision)) { + assert((p.transpose() * p).isIdentity(WEYL_TOLERANCE)); + assert(std::abs(Matrix4x4::fromDiagonal(d).determinant() - 1.0) < + WEYL_TOLERANCE); + return std::make_pair(p, d); + } + } + llvm::reportFatalInternalError(llvm::formatv( + "TwoQubitWeylDecomposition: failed to diagonalize M2 ({0} iterations). " + "best error = {1:e}, precision = {2:e}", + maxDiagonalizationAttempts, bestErr, precision)); +} + +static std::tuple +decomposeTwoQubitProductGate(const Matrix4x4& specialUnitary) { + Matrix2x2 r = + Matrix2x2::fromElements(specialUnitary(0, 0), specialUnitary(0, 1), + specialUnitary(1, 0), specialUnitary(1, 1)); + auto detR = r.determinant(); + if (std::abs(detR) < 0.1) { + r = Matrix2x2::fromElements(specialUnitary(2, 0), specialUnitary(2, 1), + specialUnitary(3, 0), specialUnitary(3, 1)); + detR = r.determinant(); + } + if (std::abs(detR) < 0.1) { + llvm::reportFatalInternalError( + "decomposeTwoQubitProductGate: unable to decompose: det_r < 0.1"); + } + r *= (1.0 / std::sqrt(detR)); + const Matrix2x2 rTConj = r.adjoint(); + + const Matrix4x4 temp = + specialUnitary * Matrix4x4::kron(Matrix2x2::identity(), rTConj); + + Matrix2x2 l = + Matrix2x2::fromElements(temp(0, 0), temp(0, 2), temp(2, 0), temp(2, 2)); + auto detL = l.determinant(); + if (std::abs(detL) < 0.9) { + llvm::reportFatalInternalError( + "decomposeTwoQubitProductGate: unable to decompose: detL < 0.9"); + } + l *= (1.0 / std::sqrt(detL)); + const auto phase = std::arg(detL) / 2.; + + return {l, r, phase}; +} + +static std::complex getTrace(double a, double b, double c, double ap, + double bp, double cp) { + const auto da = a - ap; + const auto db = b - bp; + const auto dc = c - cp; + return 4. * std::complex{std::cos(da) * std::cos(db) * std::cos(dc), + std::sin(da) * std::sin(db) * std::sin(dc)}; +} + +static Specialization +bestSpecialization(const TwoQubitWeylDecomposition& decomposition, + const std::optional& requestedFidelity) { + auto isClose = [&](double ap, double bp, double cp) -> bool { + const auto tr = getTrace(decomposition.a(), decomposition.b(), + decomposition.c(), ap, bp, cp); + if (requestedFidelity) { + return traceToFidelity(tr) >= *requestedFidelity; + } + return false; + }; + + const auto closestAbc = closestPartialSwap( + decomposition.a(), decomposition.b(), decomposition.c()); + const auto closestAbMinusC = closestPartialSwap( + decomposition.a(), decomposition.b(), -decomposition.c()); + + if (isClose(0., 0., 0.)) { + return Specialization::IdEquiv; + } + if (isClose(PI_OVER_4, PI_OVER_4, PI_OVER_4) || + isClose(PI_OVER_4, PI_OVER_4, -PI_OVER_4)) { + return Specialization::SWAPEquiv; + } + if (isClose(closestAbc, closestAbc, closestAbc)) { + return Specialization::PartialSWAPEquiv; + } + if (isClose(closestAbMinusC, closestAbMinusC, -closestAbMinusC)) { + return Specialization::PartialSWAPFlipEquiv; + } + if (isClose(decomposition.a(), 0., 0.)) { + return Specialization::ControlledEquiv; + } + if (isClose(PI_OVER_4, PI_OVER_4, decomposition.c())) { + return Specialization::MirrorControlledEquiv; + } + if (isClose((decomposition.a() + decomposition.b()) / 2., + (decomposition.a() + decomposition.b()) / 2., + decomposition.c())) { + return Specialization::FSimaabEquiv; + } + if (isClose(decomposition.a(), (decomposition.b() + decomposition.c()) / 2., + (decomposition.b() + decomposition.c()) / 2.)) { + return Specialization::FSimabbEquiv; + } + if (isClose(decomposition.a(), (decomposition.b() - decomposition.c()) / 2., + (decomposition.c() - decomposition.b()) / 2.)) { + return Specialization::FSimabmbEquiv; + } + return Specialization::General; +} + +static std::pair projectToSU4(const Matrix4x4& unitary) { + auto u = unitary; + const auto detU = u.determinant(); + u *= std::pow(detU, -0.25); + return {u, std::arg(detU) / 4.0}; +} + +static std::tuple, + std::array> +computeOrderedWeylCoordinates(const Matrix4x4& u) { + const auto uP = magicBasisTransform(u, /*outOfMagicBasis=*/true); + const Matrix4x4 m2 = uP.transpose() * uP; + auto [p, d] = diagonalizeComplexSymmetric(m2, DIAGONALIZATION_PRECISION); + + std::array dReal{}; + for (std::size_t i = 0; i < d.size(); ++i) { + dReal[i] = -std::arg(d[i]) / 2.0; + } + dReal[3] = -dReal[0] - dReal[1] - dReal[2]; + + std::array cs{}; + for (std::size_t i = 0; i < cs.size(); ++i) { + cs[i] = remEuclid((dReal[i] + dReal[3]) / 2.0, 2.0 * PI); + } + + // Sort coordinates by min(x mod pi/2, pi/2 - x mod pi/2). + std::array cstemp{}; + for (std::size_t i = 0; i < cs.size(); ++i) { + const auto tmp = remEuclid(cs[i], PI / 2.0); + cstemp[i] = std::min(tmp, (PI / 2.0) - tmp); + } + std::array order{0, 1, 2}; + std::ranges::stable_sort( + order, [&](auto a, auto b) { return cstemp[a] < cstemp[b]; }); + order = {order[1], order[2], order[0]}; + cs = {cs[order[0]], cs[order[1]], cs[order[2]]}; + { + const std::array reordered{dReal[order[0]], dReal[order[1]], + dReal[order[2]]}; + dReal[0] = reordered[0]; + dReal[1] = reordered[1]; + dReal[2] = reordered[2]; + } + + const Matrix4x4 pOrig = p; + for (std::size_t i = 0; i < order.size(); ++i) { + p.setColumn(i, pOrig.column(order[i])); + } + if (p.determinant().real() < 0.0) { + auto lastColumn = p.column(3); + for (auto& entry : lastColumn) { + entry = -entry; + } + p.setColumn(3, lastColumn); + } + assert(std::abs(p.determinant() - 1.0) < WEYL_TOLERANCE); + + return {uP, p, cs, dReal}; +} + +static ChamberState buildChamberState(const Matrix4x4& u, const Matrix4x4& uP, + Matrix4x4 p, std::array cs, + const std::array& dReal, + double globalPhase) { + std::array tempDiag{}; + for (std::size_t k = 0; k < tempDiag.size(); ++k) { + tempDiag[k] = std::exp(1i * dReal[k]); + } + const Matrix4x4 temp = Matrix4x4::fromDiagonal(tempDiag); + + Matrix4x4 k1 = uP * p * temp; + // Orthogonality checks use WEYL_TOLERANCE (diagonalization residual scale). + assert((k1.transpose() * k1).isIdentity(WEYL_TOLERANCE)); + assert(k1.determinant().real() > 0.0); + k1 = magicBasisTransform(k1, /*outOfMagicBasis=*/false); + + Matrix4x4 k2 = p.adjoint(); + assert((k2.transpose() * k2).isIdentity(WEYL_TOLERANCE)); + assert(k2.determinant().real() > 0.0); + k2 = magicBasisTransform(k2, /*outOfMagicBasis=*/false); + + std::array tempConjDiag{}; + for (std::size_t k = 0; k < tempConjDiag.size(); ++k) { + tempConjDiag[k] = std::conj(tempDiag[k]); + } + assert((k1 * + magicBasisTransform(Matrix4x4::fromDiagonal(tempConjDiag), + /*outOfMagicBasis=*/false) * + k2) + .isApprox(u, WEYL_TOLERANCE)); + + auto [k1l, k1r, phaseL] = decomposeTwoQubitProductGate(k1); + auto [k2l, k2r, phaseR] = decomposeTwoQubitProductGate(k2); + assert(Matrix4x4::kron(k1l, k1r).isApprox(k1, WEYL_TOLERANCE)); + assert(Matrix4x4::kron(k2l, k2r).isApprox(k2, WEYL_TOLERANCE)); + globalPhase += phaseL + phaseR; + + if (cs[0] > (PI / 2.0)) { + cs[0] -= 3.0 * (PI / 2.0); + k1l = k1l * iPauliY(); + k1r = k1r * iPauliY(); + globalPhase += (PI / 2.0); + } + if (cs[1] > (PI / 2.0)) { + cs[1] -= 3.0 * (PI / 2.0); + k1l = k1l * iPauliX(); + k1r = k1r * iPauliX(); + globalPhase += (PI / 2.0); + } + auto conjs = 0; + if (cs[0] > PI_OVER_4) { + cs[0] = (PI / 2.0) - cs[0]; + k1l = k1l * iPauliY(); + k2r = iPauliY() * k2r; + conjs += 1; + globalPhase -= (PI / 2.0); + } + if (cs[1] > PI_OVER_4) { + cs[1] = (PI / 2.0) - cs[1]; + k1l = k1l * iPauliX(); + k2r = iPauliX() * k2r; + conjs += 1; + globalPhase += (PI / 2.0); + if (conjs == 1) { + globalPhase -= PI; + } + } + if (cs[2] > (PI / 2.0)) { + cs[2] -= 3.0 * (PI / 2.0); + k1l = k1l * iPauliZ(); + k1r = k1r * iPauliZ(); + globalPhase += (PI / 2.0); + if (conjs == 1) { + globalPhase -= PI; + } + } + if (conjs == 1) { + cs[2] = (PI / 2.0) - cs[2]; + k1l = k1l * iPauliZ(); + k2r = iPauliZ() * k2r; + globalPhase += (PI / 2.0); + } + if (cs[2] > PI_OVER_4) { + cs[2] -= (PI / 2.0); + k1l = k1l * iPauliZ(); + k1r = k1r * iPauliZ(); + globalPhase -= (PI / 2.0); + } + + ChamberState chamber; + chamber.cs = cs; + chamber.k1l = k1l; + chamber.k1r = k1r; + chamber.k2l = k2l; + chamber.k2r = k2r; + chamber.globalPhase = globalPhase; + chamber.a = cs[1]; + chamber.b = cs[0]; + chamber.c = cs[2]; + return chamber; +} + +//===----------------------------------------------------------------------===// +// TwoQubitWeylDecomposition +//===----------------------------------------------------------------------===// + +void TwoQubitWeylDecomposition::finalizeSpecializationPhase( + bool flippedFromOriginal, double preSpecializationA, + double preSpecializationB, double preSpecializationC, + const std::optional& fidelity) { + const auto trace = + flippedFromOriginal + ? getTrace((PI / 2.0) - preSpecializationA, preSpecializationB, + -preSpecializationC, a_, b_, c_) + : getTrace(preSpecializationA, preSpecializationB, preSpecializationC, + a_, b_, c_); + const double calculatedFidelity = traceToFidelity(trace); + if (fidelity && calculatedFidelity + 1.0e-13 < *fidelity) { + llvm::reportFatalInternalError(llvm::formatv( + "TwoQubitWeylDecomposition: Calculated fidelity of " + "specialization is worse than requested fidelity ({0:F4} vs {1:F4})!", + calculatedFidelity, *fidelity)); + } + globalPhase_ += std::arg(trace); +} + +TwoQubitWeylDecomposition +TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, + std::optional fidelity) { + if (fidelity && + (!std::isfinite(*fidelity) || *fidelity < 0.0 || *fidelity > 1.0)) { + llvm::reportFatalInternalError(llvm::formatv( + "TwoQubitWeylDecomposition: fidelity must be finite and in [0, 1] " + "(got {0})", + *fidelity)); + } + + const auto [u, globalPhase0] = projectToSU4(unitaryMatrix); + auto [uP, p, cs, dReal] = computeOrderedWeylCoordinates(u); + const auto chamber = buildChamberState(u, uP, p, cs, dReal, globalPhase0); + TwoQubitWeylDecomposition decomposition; + decomposition.a_ = chamber.a; + decomposition.b_ = chamber.b; + decomposition.c_ = chamber.c; + decomposition.globalPhase_ = chamber.globalPhase; + decomposition.k1l_ = chamber.k1l; + decomposition.k2l_ = chamber.k2l; + decomposition.k1r_ = chamber.k1r; + decomposition.k2r_ = chamber.k2r; + + assert((Matrix4x4::kron(decomposition.k1l_, decomposition.k1r_) * + decomposition.getCanonicalMatrix() * + Matrix4x4::kron(decomposition.k2l_, decomposition.k2r_) * + std::exp(Complex{0.0, 1.0} * decomposition.globalPhase_)) + .isApprox(unitaryMatrix, WEYL_TOLERANCE)); + + const bool flippedFromOriginal = decomposition.applySpecialization(fidelity); + decomposition.finalizeSpecializationPhase(flippedFromOriginal, chamber.a, + chamber.b, chamber.c, fidelity); + + const auto reconstructed = + Matrix4x4::kron(decomposition.k1l_, decomposition.k1r_) * + decomposition.getCanonicalMatrix() * + Matrix4x4::kron(decomposition.k2l_, decomposition.k2r_) * + std::exp(Complex{0.0, 1.0} * decomposition.globalPhase_); + if (!reconstructed.isApprox(unitaryMatrix, WEYL_TOLERANCE)) { + llvm::reportFatalInternalError( + "TwoQubitWeylDecomposition: failed to reconstruct unitary after " + "specialization"); + } + + return decomposition; +} + +Matrix4x4 TwoQubitWeylDecomposition::getCanonicalMatrix(double a, double b, + double c) { + return rzzMatrix(-2.0 * c) * ryyMatrix(-2.0 * b) * rxxMatrix(-2.0 * a); +} + +bool TwoQubitWeylDecomposition::applySpecialization( + const std::optional& requestedFidelity) { + bool flippedFromOriginal = false; + const auto newSpecialization = bestSpecialization(*this, requestedFidelity); + if (newSpecialization == Specialization::General) { + return flippedFromOriginal; + } + + switch (newSpecialization) { + case Specialization::IdEquiv: + a_ = 0.; + b_ = 0.; + c_ = 0.; + k1l_ = k1l_ * k2l_; + k2l_ = Matrix2x2::identity(); + k1r_ = k1r_ * k2r_; + k2r_ = Matrix2x2::identity(); + break; + case Specialization::SWAPEquiv: + if (c_ > 0.) { + k1l_ = k1l_ * k2r_; + k1r_ = k1r_ * k2l_; + k2l_ = Matrix2x2::identity(); + k2r_ = Matrix2x2::identity(); + } else { + flippedFromOriginal = true; + globalPhase_ += (PI / 2.0); + k1l_ = k1l_ * iPauliZ() * k2r_; + k1r_ = k1r_ * iPauliZ() * k2l_; + k2l_ = Matrix2x2::identity(); + k2r_ = Matrix2x2::identity(); + } + a_ = PI_OVER_4; + b_ = PI_OVER_4; + c_ = PI_OVER_4; + break; + case Specialization::PartialSWAPEquiv: { + const auto closest = closestPartialSwap(a_, b_, c_); + const auto k2lDagger = k2l_.adjoint(); + a_ = closest; + b_ = closest; + c_ = closest; + k1l_ = k1l_ * k2l_; + k1r_ = k1r_ * k2l_; + k2r_ = k2lDagger * k2r_; + k2l_ = Matrix2x2::identity(); + break; + } + case Specialization::PartialSWAPFlipEquiv: { + const auto closest = closestPartialSwap(a_, b_, -c_); + const auto k2lDagger = k2l_.adjoint(); + a_ = closest; + b_ = closest; + c_ = -closest; + k1l_ = k1l_ * k2l_; + k1r_ = k1r_ * iPauliZ() * k2l_ * iPauliZ(); + k2r_ = iPauliZ() * k2lDagger * iPauliZ() * k2r_; + k2l_ = Matrix2x2::identity(); + break; + } + case Specialization::ControlledEquiv: { + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::XYX); + const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = + anglesFromUnitary(k2r_, EulerBasis::XYX); + b_ = 0.; + c_ = 0.; + globalPhase_ = globalPhase_ + k2lphase + k2rphase; + k1l_ = k1l_ * rxMatrix(k2lphi); + k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); + k1r_ = k1r_ * rxMatrix(k2rphi); + k2r_ = ryMatrix(k2rtheta) * rxMatrix(k2rlambda); + break; + } + case Specialization::MirrorControlledEquiv: { + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::ZYZ); + const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = + anglesFromUnitary(k2r_, EulerBasis::ZYZ); + a_ = PI_OVER_4; + b_ = PI_OVER_4; + globalPhase_ = globalPhase_ + k2lphase + k2rphase; + k1l_ = k1l_ * rzMatrix(k2rphi); + k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); + k1r_ = k1r_ * rzMatrix(k2lphi); + k2r_ = ryMatrix(k2rtheta) * rzMatrix(k2rlambda); + break; + } + case Specialization::FSimaabEquiv: { + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::ZYZ); + const auto ab = (a_ + b_) / 2.; + a_ = ab; + b_ = ab; + globalPhase_ = globalPhase_ + k2lphase; + k1l_ = k1l_ * rzMatrix(k2lphi); + k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); + k1r_ = k1r_ * rzMatrix(k2lphi); + k2r_ = rzMatrix(-k2lphi) * k2r_; + break; + } + case Specialization::FSimabbEquiv: { + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::XYX); + const auto bc = (b_ + c_) / 2.; + b_ = bc; + c_ = bc; + globalPhase_ = globalPhase_ + k2lphase; + k1l_ = k1l_ * rxMatrix(k2lphi); + k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); + k1r_ = k1r_ * rxMatrix(k2lphi); + k2r_ = rxMatrix(-k2lphi) * k2r_; + break; + } + case Specialization::FSimabmbEquiv: { + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::XYX); + const auto bc = (b_ - c_) / 2.; + b_ = bc; + c_ = -bc; + globalPhase_ = globalPhase_ + k2lphase; + k1l_ = k1l_ * rxMatrix(k2lphi); + k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); + k1r_ = k1r_ * iPauliZ() * rxMatrix(k2lphi) * iPauliZ(); + k2r_ = iPauliZ() * rxMatrix(-k2lphi) * iPauliZ() * k2r_; + break; + } + case Specialization::General: + llvm_unreachable("unreachable specialization"); + } + return flippedFromOriginal; +} + +} // namespace mlir::qco::decomposition diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 50d524daa1..274b3417e3 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -1015,4 +1015,71 @@ DynamicMatrix operator*(const Complex& scalar, const DynamicMatrix& matrix) { return matrix * scalar; } +Matrix2x2 rxMatrix(const double theta) { + const auto halfTheta = theta / 2.0; + const Complex cos{std::cos(halfTheta), 0.0}; + const Complex isin{0.0, -std::sin(halfTheta)}; + return Matrix2x2::fromElements(cos, isin, isin, cos); +} + +Matrix2x2 ryMatrix(const double theta) { + const auto halfTheta = theta / 2.0; + const Complex cos{std::cos(halfTheta), 0.0}; + const Complex sin{std::sin(halfTheta), 0.0}; + return Matrix2x2::fromElements(cos, -sin, sin, cos); +} + +Matrix2x2 rzMatrix(const double theta) { + return Matrix2x2::fromElements( + Complex{std::cos(theta / 2.0), -std::sin(theta / 2.0)}, 0.0, 0.0, + Complex{std::cos(theta / 2.0), std::sin(theta / 2.0)}); +} + +const Matrix2x2& iPauliX() { + static const Matrix2x2 MATRIX = + Matrix2x2::fromElements(0.0, Complex{0.0, 1.0}, Complex{0.0, 1.0}, 0.0); + return MATRIX; +} + +const Matrix2x2& iPauliY() { + static const Matrix2x2 MATRIX = Matrix2x2::fromElements(0.0, 1.0, -1.0, 0.0); + return MATRIX; +} + +const Matrix2x2& iPauliZ() { + static const Matrix2x2 MATRIX = + Matrix2x2::fromElements(Complex{0.0, 1.0}, 0.0, 0.0, Complex{0.0, -1.0}); + return MATRIX; +} + +Matrix4x4 rxxMatrix(const double theta) { + const auto cosTheta = std::cos(theta / 2.0); + const Complex misin{0.0, -std::sin(theta / 2.0)}; + return Matrix4x4::fromElements(cosTheta, 0.0, 0.0, misin, // + 0.0, cosTheta, misin, 0.0, // + 0.0, misin, cosTheta, 0.0, // + misin, 0.0, 0.0, cosTheta); +} + +Matrix4x4 ryyMatrix(const double theta) { + const auto cosTheta = std::cos(theta / 2.0); + const Complex isin{0.0, std::sin(theta / 2.0)}; + const Complex misin{0.0, -std::sin(theta / 2.0)}; + return Matrix4x4::fromElements(cosTheta, 0.0, 0.0, isin, // + 0.0, cosTheta, misin, 0.0, // + 0.0, misin, cosTheta, 0.0, // + isin, 0.0, 0.0, cosTheta); +} + +Matrix4x4 rzzMatrix(const double theta) { + const auto cosTheta = std::cos(theta / 2.0); + const auto sinTheta = std::sin(theta / 2.0); + const Complex em{cosTheta, -sinTheta}; + const Complex ep{cosTheta, sinTheta}; + return Matrix4x4::fromElements(em, 0.0, 0.0, 0.0, // + 0.0, ep, 0.0, 0.0, // + 0.0, 0.0, ep, 0.0, // + 0.0, 0.0, 0.0, em); +} + } // namespace mlir::qco diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/CMakeLists.txt b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/CMakeLists.txt index f493bb9e4d..c6160821cd 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/CMakeLists.txt +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/CMakeLists.txt @@ -7,12 +7,20 @@ # Licensed under the MIT License set(target_name mqt-core-mlir-unittest-decomposition) -add_executable(${target_name} test_euler_decomposition.cpp) +add_executable(${target_name} test_euler_decomposition.cpp test_weyl_decomposition.cpp) -target_link_libraries(${target_name} PRIVATE GTest::gtest_main MLIRQCOProgramBuilder - MLIRQCOTransforms) -target_link_libraries(${target_name} PRIVATE MLIRPass MLIRFuncDialect MLIRArithDialect MLIRIR - MLIRSupport MLIRQTensorDialect) +target_link_libraries( + ${target_name} + PRIVATE GTest::gtest_main + MLIRQCOProgramBuilder + MLIRQCOTransforms + MLIRQCOUtils + MLIRPass + MLIRFuncDialect + MLIRArithDialect + MLIRIR + MLIRSupport + MLIRQTensorDialect) mqt_mlir_configure_unittest_target(${target_name}) diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp index fa9ba1d4ea..aedc27d80f 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -104,18 +105,6 @@ class EulerSynthesisExactTest // Euler synthesis support //===----------------------------------------------------------------------===// -[[nodiscard]] static Matrix2x2 rzMatrix(const double theta) { - const auto m00 = std::polar(1.0, -theta / 2.0); - const auto m11 = std::polar(1.0, theta / 2.0); - return Matrix2x2::fromElements(m00, 0, 0, m11); -} - -[[nodiscard]] static Matrix2x2 ryMatrix(const double theta) { - const auto m00 = std::cos(theta / 2.0); - const auto m01 = -std::sin(theta / 2.0); - return Matrix2x2::fromElements(m00, m01, -m01, m00); -} - [[nodiscard]] static Matrix2x2 randomUnitaryMatrix(std::mt19937& rng) { std::uniform_real_distribution dist(-std::numbers::pi, std::numbers::pi); const Matrix2x2 su2 = @@ -453,6 +442,74 @@ TEST(EulerSynthesisTest, RandomReconstructionAllBases) { } } +TEST(EulerAnglesCoverageTest, ParamsZYZUsesOffDiagonal01When10IsNearZero) { + Matrix2x2 matrix = rxMatrix(0.4); + matrix(1, 0) = Complex{0.0, 0.0}; + ASSERT_GT(std::abs(matrix(0, 1)), mlir::utils::TOLERANCE); + const EulerAngles angles = anglesFromUnitary(matrix, ZYZ); + EXPECT_TRUE(std::isfinite(angles.theta)); + EXPECT_TRUE(std::isfinite(angles.phi)); + EXPECT_TRUE(std::isfinite(angles.lambda)); +} + +TEST(EulerAnglesCoverageTest, PhaseOnlyDecompositionSkipsRotationGates) { + TestFixture fx; + fx.setUp(); + constexpr double scale = 1.0 + 1e-10; + const Matrix2x2 matrix = Matrix2x2::fromElements(scale, 0, 0, scale); + ASSERT_FALSE(matrix.isApprox(Matrix2x2::identity())); + const EulerAngles angles = anglesFromUnitary(matrix, ZYZ); + EXPECT_LE(std::abs(angles.theta), mlir::utils::TOLERANCE); + EXPECT_LE(std::abs(angles.phi), mlir::utils::TOLERANCE); + EXPECT_LE(std::abs(angles.lambda), mlir::utils::TOLERANCE); + const auto circuit = synthesizeMatrix(fx.ctx(), matrix, ZYZ); + ASSERT_TRUE(succeeded(verify(*circuit.mlirModule))); + EXPECT_EQ(countZYZGates(circuit.func), 0U); +} + +TEST(EulerAnglesCoverageTest, UBasisZeroThetaEmitsSingleUGate) { + TestFixture fx; + fx.setUp(); + const Matrix2x2 matrix = rzMatrix(0.7); + expectSynthesizedMatrix(fx.ctx(), matrix, U, + [](func::FuncOp funcOp, const Matrix2x2& /*matrix*/) { + EXPECT_EQ(countOps(funcOp), 1U); + EXPECT_EQ(countZYZGates(funcOp), 0U); + }); +} + +TEST(EulerAnglesCoverageTest, UBasisNonzeroThetaEmitsSingleUGate) { + TestFixture fx; + fx.setUp(); + const Matrix2x2 matrix = ryMatrix(1.2); + expectSynthesizedMatrix(fx.ctx(), matrix, U, + [](func::FuncOp funcOp, const Matrix2x2& /*matrix*/) { + EXPECT_EQ(countOps(funcOp), 1U); + EXPECT_EQ(countZYZGates(funcOp), 0U); + }); +} + +TEST(EulerAnglesCoverageTest, Mod2PiMapsPiBoundaryThroughSynthesis) { + TestFixture fx; + fx.setUp(); + constexpr double eps = 0.5 * mlir::utils::TOLERANCE; + const Complex global = std::polar(1.0, std::numbers::pi - eps); + const Matrix2x2 matrix = Matrix2x2::fromElements(global, 0, 0, global); + expectSynthesizedMatrix(fx.ctx(), matrix, U, + [](func::FuncOp funcOp, const Matrix2x2& /*matrix*/) { + EXPECT_EQ(countOps(funcOp), 1U); + EXPECT_EQ(countOps(funcOp), 1U); + }); +} + +TEST(EulerAnglesCoverageTest, Mod2PiPreservesNonFinitePhase) { + TestFixture fx; + fx.setUp(); + const Matrix2x2 matrix = Matrix2x2::fromElements( + Complex{std::numeric_limits::quiet_NaN(), 0}, 0, 0, 1); + EXPECT_NO_FATAL_FAILURE(synthesizeMatrix(fx.ctx(), matrix, ZYZ)); +} + //===----------------------------------------------------------------------===// // FuseSingleQubitUnitaryRuns support //===----------------------------------------------------------------------===// diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp new file mode 100644 index 0000000000..bb11c6fb2b --- /dev/null +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -0,0 +1,409 @@ +/* + * 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 + */ + +#include "mlir/Dialect/QCO/IR/QCOOps.h" +#include "mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h" +#include "mlir/Dialect/QCO/Utils/Matrix.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace mlir::qco; +using namespace mlir::qco::decomposition; + +static Complex globalPhaseFactor(const double phase) { + return std::exp(Complex{0.0, 1.0} * phase); +} + +static const Matrix4x4& twoQubitControlledX01() { + static const Matrix4x4 MATRIX = Matrix4x4::fromElements(1.0, 0.0, 0.0, 0.0, // + 0.0, 1.0, 0.0, 0.0, // + 0.0, 0.0, 0.0, 1.0, // + 0.0, 0.0, 1.0, 0.0); + return MATRIX; +} + +static const Matrix4x4& twoQubitControlledX10() { + static const Matrix4x4 MATRIX = Matrix4x4::fromElements(1.0, 0.0, 0.0, 0.0, // + 0.0, 0.0, 0.0, 1.0, // + 0.0, 0.0, 1.0, 0.0, // + 0.0, 1.0, 0.0, 0.0); + return MATRIX; +} + +static bool isUnitaryMatrix(const Matrix2x2& matrix, + const double tolerance = MATRIX_TOLERANCE) { + return (matrix.adjoint() * matrix).isIdentity(tolerance); +} + +static Matrix4x4 randomUnitary4x4(std::mt19937& rng) { + std::normal_distribution normalDist(0.0, 1.0); + std::vector columns(4, std::vector(4, std::complex{0.0, 0.0})); + for (auto& column : columns) { + for (auto& entry : column) { + entry = std::complex(normalDist(rng), normalDist(rng)); + } + } + for (std::size_t j = 0; j < 4; ++j) { + for (std::size_t k = 0; k < j; ++k) { + std::complex projection{0.0, 0.0}; + for (std::size_t i = 0; i < 4; ++i) { + projection += std::conj(columns[k][i]) * columns[j][i]; + } + for (std::size_t i = 0; i < 4; ++i) { + columns[j][i] -= projection * columns[k][i]; + } + } + double norm = 0.0; + for (std::size_t i = 0; i < 4; ++i) { + norm += std::norm(columns[j][i]); + } + norm = std::sqrt(norm); + for (std::size_t i = 0; i < 4; ++i) { + columns[j][i] /= norm; + } + } + const Matrix4x4 unitary = Matrix4x4::fromElements( + columns[0][0], columns[1][0], columns[2][0], columns[3][0], columns[0][1], + columns[1][1], columns[2][1], columns[3][1], columns[0][2], columns[1][2], + columns[2][2], columns[3][2], columns[0][3], columns[1][3], columns[2][3], + columns[3][3]); + assert((unitary.adjoint() * unitary).isIdentity(WEYL_TOLERANCE)); + return unitary; +} + +static Matrix4x4 restoreWeyl(const TwoQubitWeylDecomposition& decomposition) { + return Matrix4x4::kron(decomposition.k1l(), decomposition.k1r()) * + decomposition.getCanonicalMatrix() * + Matrix4x4::kron(decomposition.k2l(), decomposition.k2r()) * + globalPhaseFactor(decomposition.globalPhase()); +} + +static Matrix4x4 restoreBasis(const TwoQubitNativeDecomposition& decomposition, + const Matrix4x4& entangler) { + const auto& factors = decomposition.singleQubitFactors; + const auto layer = [&](std::size_t i) { + return Matrix4x4::kron(factors[(2 * i) + 1], factors[2 * i]); + }; + Matrix4x4 matrix = layer(0); + for (std::uint8_t i = 0; i < decomposition.numBasisUses; ++i) { + matrix = entangler * matrix; + matrix = layer(static_cast(i) + 1) * matrix; + } + return matrix * globalPhaseFactor(decomposition.globalPhase); +} + +static auto productMatrixCases() { + return ::testing::Values( + []() { return Matrix4x4::identity(); }, + []() { return Matrix4x4::kron(rzMatrix(1.0), ryMatrix(3.1)); }, + []() { return Matrix4x4::kron(Matrix2x2::identity(), rxMatrix(0.1)); }); +} + +static auto entangledMatrixCases() { + return ::testing::Values( + []() { return rzzMatrix(2.0); }, + []() { return ryyMatrix(1.0) * rzzMatrix(3.0) * rxxMatrix(2.0); }, + []() { + return TwoQubitWeylDecomposition::getCanonicalMatrix(1.5, -0.2, 0.0) * + Matrix4x4::kron(rxMatrix(1.0), Matrix2x2::identity()); + }, + []() { + return Matrix4x4::kron(rxMatrix(1.0), ryMatrix(1.0)) * + TwoQubitWeylDecomposition::getCanonicalMatrix(1.1, 0.2, 3.0) * + Matrix4x4::kron(rxMatrix(1.0), Matrix2x2::identity()); + }, + []() { + return Matrix4x4::kron(HOp::getUnitaryMatrix(), iPauliZ()) * + twoQubitControlledX01() * Matrix4x4::kron(iPauliX(), iPauliY()); + }); +} + +static auto cxBasisCases() { + return ::testing::Values([]() { return twoQubitControlledX01(); }, + []() { return twoQubitControlledX10(); }); +} + +static auto specializedMatrixCases() { + return ::testing::Values( + []() { + return twoQubitControlledX01() * twoQubitControlledX10() * + twoQubitControlledX01(); + }, + []() { + return TwoQubitWeylDecomposition::getCanonicalMatrix(0.5, 0.5, 0.5); + }, + []() { + return TwoQubitWeylDecomposition::getCanonicalMatrix(0.5, 0.5, -0.5); + }, + []() { return twoQubitControlledX01() * twoQubitControlledX10(); }, + []() { + return TwoQubitWeylDecomposition::getCanonicalMatrix(0.5, 0.5, 0.1); + }, + []() { + return TwoQubitWeylDecomposition::getCanonicalMatrix(0.5, 0.1, 0.1); + }, + []() { + return TwoQubitWeylDecomposition::getCanonicalMatrix(0.5, 0.1, -0.1); + }); +} + +TEST(DecompositionHelpersTest, MatrixUtilitySanity) { + EXPECT_NEAR(std::abs(globalPhaseFactor(1.25)), 1.0, 1e-14); + EXPECT_FALSE(isUnitaryMatrix(Matrix2x2::fromElements(2.0, 0.0, 0.0, 2.0))); + EXPECT_TRUE(isUnitaryMatrix(Matrix2x2::identity())); +} + +TEST(DecompositionHelpersTest, GateMatrixFactoriesMatchCanonicalForm) { + for (const double theta : {0.0, 0.25, 1.0, 2.5, -1.3}) { + EXPECT_TRUE(rxxMatrix(theta).isApprox( + TwoQubitWeylDecomposition::getCanonicalMatrix(-theta / 2.0, 0.0, 0.0), + WEYL_TOLERANCE)); + EXPECT_TRUE(ryyMatrix(theta).isApprox( + TwoQubitWeylDecomposition::getCanonicalMatrix(0.0, -theta / 2.0, 0.0), + WEYL_TOLERANCE)); + EXPECT_TRUE(rzzMatrix(theta).isApprox( + TwoQubitWeylDecomposition::getCanonicalMatrix(0.0, 0.0, -theta / 2.0), + WEYL_TOLERANCE)); + } +} + +namespace { + +class WeylDecompositionTest : public testing::TestWithParam {}; + +class BasisDecomposerTest : public testing::TestWithParam< + std::tuple> { +protected: + void SetUp() override { + basisMatrix = std::get<0>(GetParam())(); + target = std::get<1>(GetParam())(); + targetDecomposition = std::make_unique( + TwoQubitWeylDecomposition::create(target, 1.0)); + } + + Matrix4x4 target; + Matrix4x4 basisMatrix; + std::unique_ptr targetDecomposition; +}; + +} // namespace + +TEST_P(WeylDecompositionTest, ReconstructsWithinRequestedFidelity) { + const Matrix4x4 originalMatrix = GetParam()(); + for (const double fidelity : {1.0, 1.0 - 1e-12}) { + const auto decomposition = + TwoQubitWeylDecomposition::create(originalMatrix, fidelity); + EXPECT_TRUE( + restoreWeyl(decomposition).isApprox(originalMatrix, WEYL_TOLERANCE)); + } +} + +TEST(WeylDecompositionStandalone, + CnotProducesValidWeylParametersAndUnitaryLocals) { + const Matrix4x4 cnot = + Matrix4x4::fromElements(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0); + const auto decomp = TwoQubitWeylDecomposition::create(cnot, std::nullopt); + constexpr double piOver4 = 0.7853981633974483; + for (const double angle : {decomp.a(), decomp.b(), decomp.c()}) { + EXPECT_GE(angle, -1e-10); + EXPECT_LE(angle, piOver4 + 1e-10); + } + EXPECT_TRUE(isUnitaryMatrix(decomp.k1l())); + EXPECT_TRUE(isUnitaryMatrix(decomp.k2l())); + EXPECT_TRUE(isUnitaryMatrix(decomp.k1r())); + EXPECT_TRUE(isUnitaryMatrix(decomp.k2r())); +} + +TEST(WeylDecompositionStandalone, Random) { + std::mt19937 rng{1234567UL}; + for (int i = 0; i < 5000; ++i) { + const Matrix4x4 originalMatrix = randomUnitary4x4(rng); + const auto decomposition = TwoQubitWeylDecomposition::create( + originalMatrix, std::optional{1.0 - 1e-12}); + EXPECT_TRUE( + restoreWeyl(decomposition).isApprox(originalMatrix, WEYL_TOLERANCE)); + } +} + +INSTANTIATE_TEST_SUITE_P(ProductTwoQubitMatrices, WeylDecompositionTest, + productMatrixCases()); +INSTANTIATE_TEST_SUITE_P(TwoQubitMatrices, WeylDecompositionTest, + entangledMatrixCases()); +INSTANTIATE_TEST_SUITE_P(SpecializedMatrices, WeylDecompositionTest, + specializedMatrixCases()); + +TEST_P(BasisDecomposerTest, ReconstructsWithinRequestedFidelity) { + for (const double fidelity : {1.0, 1.0 - 1e-12}) { + const auto decomposer = + TwoQubitBasisDecomposer::create(basisMatrix, fidelity); + const auto decomposed = + decomposer.twoQubitDecompose(*targetDecomposition, std::nullopt); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_TRUE(restoreBasis(*decomposed, basisMatrix) + .isApprox(target, WEYL_TOLERANCE)); + } +} + +TEST(BasisDecomposerTest, Random) { + std::mt19937 rng{123456UL}; + const mlir::SmallVector basisMatrices{twoQubitControlledX01(), + twoQubitControlledX10()}; + std::uniform_int_distribution distBasisGate{0, 1}; + + for (int i = 0; i < 2000; ++i) { + const Matrix4x4 originalMatrix = randomUnitary4x4(rng); + const auto targetDecomposition = TwoQubitWeylDecomposition::create( + originalMatrix, std::optional{1.0}); + const Matrix4x4 basisMatrix = basisMatrices[distBasisGate(rng)]; + const auto decomposer = TwoQubitBasisDecomposer::create(basisMatrix, 1.0); + const auto decomposed = + decomposer.twoQubitDecompose(targetDecomposition, std::nullopt); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_TRUE(restoreBasis(*decomposed, basisMatrix) + .isApprox(originalMatrix, WEYL_TOLERANCE)); + } +} + +TEST(BasisDecomposerNumBasisTest, ForcesZeroBasisUsesForIdentityTarget) { + const Matrix4x4 basis = twoQubitControlledX01(); + const auto decomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const Matrix4x4 target = Matrix4x4::identity(); + const auto weyl = TwoQubitWeylDecomposition::create(target, 1.0); + const auto decomposed = decomposer.twoQubitDecompose(weyl, std::uint8_t{0}); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_EQ(decomposed->numBasisUses, 0); + EXPECT_TRUE( + restoreBasis(*decomposed, basis).isApprox(target, WEYL_TOLERANCE)); +} + +TEST(BasisDecomposerTest, DecomposeTwoQubitWithBasisReconstructsTarget) { + const Matrix4x4 basis = twoQubitControlledX01(); + const Matrix4x4 target = + Matrix4x4::kron(rxMatrix(0.4), ryMatrix(0.6)) * + TwoQubitWeylDecomposition::getCanonicalMatrix(0.3, 0.2, 0.1) * + Matrix4x4::kron(rzMatrix(0.2), Matrix2x2::identity()); + const auto decomposed = decomposeTwoQubitWithBasis(target, basis); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_TRUE( + restoreBasis(*decomposed, basis).isApprox(target, WEYL_TOLERANCE)); +} + +TEST(BasisDecomposerTest, CachedDecomposerMatchesOneShotAcrossTargets) { + const Matrix4x4 basis = twoQubitControlledX01(); + const auto cachedDecomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const mlir::SmallVector targets{ + Matrix4x4::identity(), + twoQubitControlledX01(), + Matrix4x4::kron(rxMatrix(0.2), ryMatrix(0.3)) * + TwoQubitWeylDecomposition::getCanonicalMatrix(0.1, 0.2, 0.3) * + Matrix4x4::kron(rzMatrix(0.1), Matrix2x2::identity()), + }; + for (const Matrix4x4& target : targets) { + const auto oneShot = decomposeTwoQubitWithBasis(target, basis); + const auto cached = cachedDecomposer.decomposeTarget(target); + ASSERT_TRUE(oneShot.has_value()); + ASSERT_TRUE(cached.has_value()); + EXPECT_TRUE(restoreBasis(*oneShot, basis).isApprox(target, WEYL_TOLERANCE)); + EXPECT_TRUE(restoreBasis(*cached, basis).isApprox(target, WEYL_TOLERANCE)); + EXPECT_EQ(cached->numBasisUses, oneShot->numBasisUses); + EXPECT_EQ(cached->singleQubitFactors.size(), + oneShot->singleQubitFactors.size()); + } +} + +TEST(BasisDecomposerTest, RejectsMultipleBasisUsesForNonSuperControlledBasis) { + const Matrix4x4 basis = rzzMatrix(1.0); + const auto decomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const auto weyl = + TwoQubitWeylDecomposition::create(Matrix4x4::identity(), 1.0); + EXPECT_FALSE(decomposer.twoQubitDecompose(weyl, std::uint8_t{2}).has_value()); +} + +TEST(BasisDecomposerTest, RejectsInvalidBasisGateUseCount) { + const Matrix4x4 basis = twoQubitControlledX01(); + const auto decomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const auto weyl = + TwoQubitWeylDecomposition::create(twoQubitControlledX01(), 1.0); + EXPECT_FALSE(decomposer.twoQubitDecompose(weyl, std::uint8_t{4}).has_value()); +} + +TEST(BasisDecomposerForcedCountTest, OneBasisUseProducesFactors) { + const Matrix4x4 basis = twoQubitControlledX01(); + const auto decomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const auto weyl = + TwoQubitWeylDecomposition::create(twoQubitControlledX01(), 1.0); + const auto decomposed = decomposer.twoQubitDecompose(weyl, std::uint8_t{1}); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_EQ(decomposed->numBasisUses, 1); + EXPECT_EQ(decomposed->singleQubitFactors.size(), 4U); +} + +TEST(BasisDecomposerForcedCountTest, TwoBasisUsesProducesFactors) { + const Matrix4x4 basis = twoQubitControlledX01(); + const auto decomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const auto weyl = + TwoQubitWeylDecomposition::create(twoQubitControlledX01(), 1.0); + const auto decomposed = decomposer.twoQubitDecompose(weyl, std::uint8_t{2}); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_EQ(decomposed->numBasisUses, 2); + EXPECT_EQ(decomposed->singleQubitFactors.size(), 6U); +} + +TEST(BasisDecomposerForcedCountTest, ThreeBasisUsesProducesFactors) { + const Matrix4x4 basis = twoQubitControlledX01(); + const auto decomposer = TwoQubitBasisDecomposer::create(basis, 1.0); + const auto weyl = + TwoQubitWeylDecomposition::create(twoQubitControlledX01(), 1.0); + const auto decomposed = decomposer.twoQubitDecompose(weyl, std::uint8_t{3}); + ASSERT_TRUE(decomposed.has_value()); + EXPECT_EQ(decomposed->numBasisUses, 3); + EXPECT_EQ(decomposed->singleQubitFactors.size(), 8U); +} + +TEST(WeylDecompositionStandalone, SwapNegativeCSpecializationReconstructs) { + constexpr double piOver4 = std::numbers::pi / 4.0; + const Matrix4x4 swapNegativeC = + TwoQubitWeylDecomposition::getCanonicalMatrix(piOver4, piOver4, -piOver4); + const auto decomposition = + TwoQubitWeylDecomposition::create(swapNegativeC, 1.0); + EXPECT_TRUE( + restoreWeyl(decomposition).isApprox(swapNegativeC, WEYL_TOLERANCE)); +} + +TEST(WeylDecompositionStandalone, ControlledSpecializationReconstructs) { + const Matrix4x4 controlledLike = + Matrix4x4::kron(rxMatrix(0.3), ryMatrix(0.4)) * + TwoQubitWeylDecomposition::getCanonicalMatrix(0.6, 0.0, 0.0) * + Matrix4x4::kron(Matrix2x2::identity(), rzMatrix(0.2)); + const auto decomposition = + TwoQubitWeylDecomposition::create(controlledLike, 1.0); + EXPECT_TRUE( + restoreWeyl(decomposition).isApprox(controlledLike, WEYL_TOLERANCE)); +} + +INSTANTIATE_TEST_SUITE_P(ProductTwoQubitMatrices, BasisDecomposerTest, + testing::Combine(cxBasisCases(), + productMatrixCases())); +INSTANTIATE_TEST_SUITE_P(TwoQubitMatrices, BasisDecomposerTest, + testing::Combine(cxBasisCases(), + entangledMatrixCases()));