From 3ca7f74c40500f3ef9bc50243e8958a2e0446757 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 18:22:28 +0200 Subject: [PATCH 01/16] =?UTF-8?q?=E2=9C=A8=20Weyl=20decomposition?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Euler.h | 25 + .../QCO/Transforms/Decomposition/Weyl.h | 187 ++++ mlir/include/mlir/Dialect/QCO/Utils/Matrix.h | 27 + .../QCO/Transforms/Decomposition/Euler.cpp | 25 +- .../QCO/Transforms/Decomposition/Weyl.cpp | 977 ++++++++++++++++++ mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 67 ++ .../Transforms/Decomposition/CMakeLists.txt | 18 +- .../test_euler_decomposition.cpp | 12 - .../Decomposition/test_weyl_decomposition.cpp | 256 +++++ 9 files changed, 1553 insertions(+), 41 deletions(-) create mode 100644 mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h create mode 100644 mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp create mode 100644 mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp 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..ac7ef42cee --- /dev/null +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h @@ -0,0 +1,187 @@ +/* + * 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 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 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 Left single-qubit factor after the canonical gate. + * + * ``` + * q1 - k2r - C - k1r - + * A + * q0 - k2l - N - k1l - + * ``` + */ + [[nodiscard]] const Matrix2x2& k1l() const { return k1l_; } + /** + * @brief Left single-qubit factor before the canonical gate. + */ + [[nodiscard]] const Matrix2x2& k2l() const { return k2l_; } + /** + * @brief Right single-qubit factor after the canonical gate. + */ + [[nodiscard]] const Matrix2x2& k1r() const { return k1r_; } + /** + * @brief Right single-qubit factor 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); + + double a_{}; + double b_{}; + double c_{}; + double globalPhase_{}; + Matrix2x2 k1l_; + Matrix2x2 k2l_; + Matrix2x2 k1r_; + Matrix2x2 k2r_; +}; + +/** + * @brief Single-qubit factors and entangler count for basis-gate synthesis. + * + * 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. + */ +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 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 + twoQubitDecompose(const TwoQubitWeylDecomposition& targetDecomposition, + std::optional 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 + 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; + + double basisFidelity{}; + TwoQubitWeylDecomposition basisDecomposer; + bool isSuperControlled{}; + SmbPrecomputed smb{}; +}; + +/** + * @brief Decomposes @p target with respect to a two-qubit basis gate matrix. + */ +[[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..60324d3873 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 Constant-folded `RX(theta)` unitary. */ +[[nodiscard]] Matrix2x2 rxMatrix(double theta); + +/** @brief Constant-folded `RY(theta)` unitary. */ +[[nodiscard]] Matrix2x2 ryMatrix(double theta); + +/** @brief Constant-folded `RZ(theta)` unitary. */ +[[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 Constant-folded two-qubit `RXX(theta)` unitary. */ +[[nodiscard]] Matrix4x4 rxxMatrix(double theta); + +/** @brief Constant-folded two-qubit `RYY(theta)` unitary. */ +[[nodiscard]] Matrix4x4 ryyMatrix(double theta); + +/** @brief Constant-folded two-qubit `RZZ(theta)` unitary. */ +[[nodiscard]] Matrix4x4 rzzMatrix(double theta); + } // namespace mlir::qco 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..8b47c02c87 --- /dev/null +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -0,0 +1,977 @@ +/* + * 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 +#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, +}; + +} // namespace + +static constexpr auto DIAGONALIZATION_PRECISION = WEYL_TOLERANCE / 10; + +/** Non-negative remainder of @p a modulo @p b (always in `[0, |b|)`). */ +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 constexpr double INV_SQRT2 = 1.0 / std::numbers::sqrt2; + +static Matrix4x4 magicBasisTransform(const Matrix4x4& unitary, + bool outOfMagicBasis) { + const Matrix4x4 bNonNormalized = Matrix4x4::fromElements( // + 1, 1i, 0, 0, // + 0, 0, 1i, 1, // + 0, 0, 1i, -1, // + 1, -1i, 0, 0); + const Matrix4x4 bNonNormalizedDagger = 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); + if (outOfMagicBasis) { + return bNonNormalizedDagger * unitary * bNonNormalized; + } + return bNonNormalized * unitary * bNonNormalizedDagger; +} + +static double closestPartialSwap(double a, double b, double c) { + auto m = (a + b + c) / 3.; + auto [am, bm, cm] = std::array{a - m, b - m, c - m}; + 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{}; + 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(); + + 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)); + 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) { + auto da = a - ap; + auto db = b - bp; + 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 { + auto tr = getTrace(decomposition.a(), decomposition.b(), decomposition.c(), + ap, bp, cp); + if (requestedFidelity) { + return traceToFidelity(tr) >= *requestedFidelity; + } + return false; + }; + + auto closestAbc = closestPartialSwap(decomposition.a(), decomposition.b(), + decomposition.c()); + auto closestAbMinusC = closestPartialSwap( + decomposition.a(), decomposition.b(), -decomposition.c()); + + if (isClose(0., 0., 0.)) { + return Specialization::IdEquiv; + } + if (isClose((std::numbers::pi / 4.0), (std::numbers::pi / 4.0), + (std::numbers::pi / 4.0)) || + isClose((std::numbers::pi / 4.0), (std::numbers::pi / 4.0), + -(std::numbers::pi / 4.0))) { + 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((std::numbers::pi / 4.0), (std::numbers::pi / 4.0), + 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 bool relativeEq(double lhs, double rhs, double epsilon, + double maxRelative) { + if (lhs == rhs) { + return true; + } + if (std::isinf(lhs) || std::isinf(rhs)) { + return false; + } + auto absDiff = std::abs(lhs - rhs); + if (absDiff <= epsilon) { + return true; + } + auto absLhs = std::abs(lhs); + auto absRhs = std::abs(rhs); + if (absRhs > absLhs) { + return absDiff <= absRhs * maxRelative; + } + return absDiff <= absLhs * maxRelative; +} + +TwoQubitWeylDecomposition +TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, + std::optional fidelity) { + auto u = unitaryMatrix; + auto detU = u.determinant(); + // Project into SU(4) by dividing out the fourth root of det(U): for a 4x4 + // unitary, |det(U)| == 1 so `det^{-1/4}` both enforces det == 1 and removes + // the global phase. The extracted phase is tracked separately in + // `globalPhase` (quarter of arg(det) to match the fourth-root choice) so the + // caller can reconstruct the original matrix exactly if needed. + auto detPow = std::pow(detU, -0.25); + u *= detPow; // remove global phase from unitary matrix + auto globalPhase = std::arg(detU) / 4.; + + // Numerical drift can still leave tiny determinant errors after root + // normalization. Re-normalize once more instead of aborting. + auto detNormalized = u.determinant(); + if (std::abs(detNormalized - Complex{1.0, 0.0}) > WEYL_TOLERANCE && + std::abs(detNormalized) > WEYL_TOLERANCE) { + u *= std::pow(detNormalized, -0.25); + } + + // transform unitary matrix to magic basis; this enables two properties: + // 1. if uP ∈ SO(4), V = A ⊗ B (SO(4) → SU(2) ⊗ SU(2)) + // 2. magic basis diagonalizes canonical gate, allowing calculation of + // canonical gate parameters later on + auto uP = magicBasisTransform(u, /*outOfMagicBasis=*/true); + const Matrix4x4 m2 = uP.transpose() * uP; + + // diagonalization yields eigenvectors (p) and eigenvalues (d); + // p is used to calculate K1/K2 (and thus the single-qubit gates + // surrounding the canonical gate); d is used to determine the Weyl + // coordinates and thus the parameters of the canonical gate + auto [p, d] = diagonalizeComplexSymmetric(m2, DIAGONALIZATION_PRECISION); + + // extract Weyl coordinates from eigenvalues, map to [0, 2*pi) + constexpr double pi = std::numbers::pi; + 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); + } + + // Reorder coordinates according to min(a, pi/2 - a) with + // a = x mod pi/2 for each Weyl coordinate x + 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]; + } + + // update eigenvectors (columns of p) according to new order of + // weyl coordinates + const Matrix4x4 pOrig = p; + for (std::size_t i = 0; i < order.size(); ++i) { + p.setColumn(i, pOrig.column(order[i])); + } + // apply correction for determinant if necessary + 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); + + // re-create complex eigenvalue matrix; this matrix contains the + // parameters of the canonical gate which is later used in the + // verification. Since the matrix is diagonal, the matrix exponential is + // equivalent to the element-wise exponential function. + 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); + + // combined matrix k1 of 1q gates after canonical gate + Matrix4x4 k1 = uP * p * temp; + // k1 must be orthogonal; the tolerance matches the iterative diagonalization + // residual rather than the (much tighter) default matrix tolerance. + assert((k1.transpose() * k1).isIdentity(WEYL_TOLERANCE)); + assert(k1.determinant().real() > 0.0); + k1 = magicBasisTransform(k1, /*outOfMagicBasis=*/false); + + // combined matrix k2 of 1q gates before canonical gate + Matrix4x4 k2 = p.adjoint(); + // k2 must be orthogonal; see the tolerance note on the k1 check above. + assert((k2.transpose() * k2).isIdentity(WEYL_TOLERANCE)); + assert(k2.determinant().real() > 0.0); + k2 = magicBasisTransform(k2, /*outOfMagicBasis=*/false); + + // ensure k1 and k2 are correct (when combined with the canonical gate + // parameters in-between, they are equivalent to u) + 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)); + + // calculate k1 = K1l ⊗ K1r + auto [K1l, K1r, phaseL] = decomposeTwoQubitProductGate(k1); + // decompose k2 = K2l ⊗ K2r + auto [K2l, K2r, phaseR] = decomposeTwoQubitProductGate(k2); + assert(Matrix4x4::kron(K1l, K1r).isApprox(k1, WEYL_TOLERANCE)); + assert(Matrix4x4::kron(K2l, K2r).isApprox(k2, WEYL_TOLERANCE)); + // accumulate global phase + globalPhase += phaseL + phaseR; + + // Flip into Weyl chamber + 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 / 4.0)) { + cs[0] = (pi / 2.0) - cs[0]; + K1l = K1l * iPauliY(); + K2r = iPauliY() * K2r; + conjs += 1; + globalPhase -= (pi / 2.0); + } + if (cs[1] > (pi / 4.0)) { + 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 / 4.0)) { + cs[2] -= (pi / 2.0); + K1l = K1l * iPauliZ(); + K1r = K1r * iPauliZ(); + globalPhase -= (pi / 2.0); + } + + // bind weyl coordinates as parameters of canonical gate + auto [a, b, c] = std::tie(cs[1], cs[0], cs[2]); + + TwoQubitWeylDecomposition decomposition; + decomposition.a_ = a; + decomposition.b_ = b; + decomposition.c_ = c; + decomposition.globalPhase_ = globalPhase; + decomposition.k1l_ = K1l; + decomposition.k2l_ = K2l; + decomposition.k1r_ = K1r; + decomposition.k2r_ = K2r; + + // make sure decomposition is equal to input + assert((Matrix4x4::kron(K1l, K1r) * decomposition.getCanonicalMatrix() * + Matrix4x4::kron(K2l, K2r) * std::exp(Complex{0.0, 1.0} * globalPhase)) + .isApprox(unitaryMatrix, WEYL_TOLERANCE)); + + // determine actual specialization of canonical gate so that the 1q + // matrices can potentially be simplified + auto flippedFromOriginal = decomposition.applySpecialization(fidelity); + + auto getTraceValue = [&]() { + if (flippedFromOriginal) { + return getTrace((pi / 2.0) - a, b, -c, decomposition.a_, decomposition.b_, + decomposition.c_); + } + return getTrace(a, b, c, decomposition.a_, decomposition.b_, + decomposition.c_); + }; + // use trace to calculate fidelity of applied specialization and + // adjust global phase + auto trace = getTraceValue(); + const double calculatedFidelity = traceToFidelity(trace); + // final check if specialization is close enough to the original matrix to + // satisfy the requested fidelity; since no forced specialization is + // allowed, this should never fail + 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)); + } + decomposition.globalPhase_ += std::arg(trace); + + // final check if decomposition is still valid after specialization + 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)); + + return decomposition; +} + +Matrix4x4 TwoQubitWeylDecomposition::getCanonicalMatrix(double a, double b, + double c) { + // Canonical gate `U_d(a, b, c) = exp(i * (a*XX + b*YY + c*ZZ))`. XX/YY/ZZ + // commute pairwise, so any product order is equivalent; the order below is + // chosen to match common Qiskit/QuantumFlow references. The negated rotation + // angles (`-2 * a`, ...) compensate for the `RXX/RYY/RZZ` convention + // `exp(-i * theta/2 * XX)`, so that the factored angles sum back to the + // intended `+a`, `+b`, `+c`. + const auto xx = rxxMatrix(-2.0 * a); + const auto yy = ryyMatrix(-2.0 * b); + const auto zz = rzzMatrix(-2.0 * c); + return zz * yy * xx; +} + +bool TwoQubitWeylDecomposition::applySpecialization( + const std::optional& requestedFidelity) { + bool flippedFromOriginal = false; + const auto newSpecialization = bestSpecialization(*this, requestedFidelity); + if (newSpecialization == Specialization::General) { + // U has no special symmetry. + // + // This gate binds all 6 possible parameters, so there is no need to + // make the single-qubit pre-/post-gates canonical. + return flippedFromOriginal; + } + + if (newSpecialization == Specialization::IdEquiv) { + // :math:`U \sim U_d(0,0,0)` + // Thus, :math:`\sim Id` + // + // This gate binds 0 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id` , :math:`K2_r = Id`. + a_ = 0.; + b_ = 0.; + c_ = 0.; + // unmodified global phase + k1l_ = k1l_ * k2l_; + k2l_ = Matrix2x2::identity(); + k1r_ = k1r_ * k2r_; + k2r_ = Matrix2x2::identity(); + } else if (newSpecialization == Specialization::SWAPEquiv) { + // :math:`U \sim U_d(\pi/4, \pi/4, \pi/4) \sim U(\pi/4, \pi/4, -\pi/4)` + // Thus, :math:`U \sim \text{SWAP}` + // + // This gate binds 0 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id` , :math:`K2_r = Id`. + if (c_ > 0.) { + // unmodified global phase + k1l_ = k1l_ * k2r_; + k1r_ = k1r_ * k2l_; + k2l_ = Matrix2x2::identity(); + k2r_ = Matrix2x2::identity(); + } else { + flippedFromOriginal = true; + + globalPhase_ += (std::numbers::pi / 2.0); + k1l_ = k1l_ * iPauliZ() * k2r_; + k1r_ = k1r_ * iPauliZ() * k2l_; + k2l_ = Matrix2x2::identity(); + k2r_ = Matrix2x2::identity(); + } + a_ = (std::numbers::pi / 4.0); + b_ = (std::numbers::pi / 4.0); + c_ = (std::numbers::pi / 4.0); + } else if (newSpecialization == Specialization::PartialSWAPEquiv) { + // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, \alpha\pi/4)` + // Thus, :math:`U \sim \text{SWAP}^\alpha` + // + // This gate binds 3 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id`. + auto closest = closestPartialSwap(a_, b_, c_); + auto k2lDagger = k2l_.adjoint(); + + a_ = closest; + b_ = closest; + c_ = closest; + // unmodified global phase + k1l_ = k1l_ * k2l_; + k1r_ = k1r_ * k2l_; + k2r_ = k2lDagger * k2r_; + k2l_ = Matrix2x2::identity(); + } else if (newSpecialization == Specialization::PartialSWAPFlipEquiv) { + // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, -\alpha\pi/4)` + // Thus, :math:`U \sim \text{SWAP}^\alpha` + // + // (a non-equivalent root of SWAP from the TwoQubitWeylPartialSWAPEquiv + // similar to how :math:`x = (\pm \sqrt(x))^2`) + // + // This gate binds 3 parameters, we make it canonical by setting: + // + // :math:`K2_l = Id` + auto closest = closestPartialSwap(a_, b_, -c_); + auto k2lDagger = k2l_.adjoint(); + + a_ = closest; + b_ = closest; + c_ = -closest; + // unmodified global phase + k1l_ = k1l_ * k2l_; + k1r_ = k1r_ * iPauliZ() * k2l_ * iPauliZ(); + k2r_ = iPauliZ() * k2lDagger * iPauliZ() * k2r_; + k2l_ = Matrix2x2::identity(); + } else if (newSpecialization == Specialization::ControlledEquiv) { + // :math:`U \sim U_d(\alpha, 0, 0)` + // Thus, :math:`U \sim \text{Ctrl-U}` + // + // This gate binds 4 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l) Rx(\lambda_l)` + // :math:`K2_r = Ry(\theta_r) Rx(\lambda_r)` + const EulerBasis eulerBasis = EulerBasis::XYX; + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, eulerBasis); + const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = + anglesFromUnitary(k2r_, EulerBasis::XYX); + // unmodified parameter a + 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); + } else if (newSpecialization == Specialization::MirrorControlledEquiv) { + // :math:`U \sim U_d(\pi/4, \pi/4, \alpha)` + // Thus, :math:`U \sim \text{SWAP} \cdot \text{Ctrl-U}` + // + // This gate binds 4 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)` + // :math:`K2_r = Ry(\theta_r)\cdot Rz(\lambda_r)` + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::ZYZ); + const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = + anglesFromUnitary(k2r_, EulerBasis::ZYZ); + a_ = (std::numbers::pi / 4.0); + b_ = (std::numbers::pi / 4.0); + // unmodified parameter c + globalPhase_ = globalPhase_ + k2lphase + k2rphase; + k1l_ = k1l_ * rzMatrix(k2rphi); + k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); + k1r_ = k1r_ * rzMatrix(k2lphi); + k2r_ = ryMatrix(k2rtheta) * rzMatrix(k2rlambda); + } else if (newSpecialization == Specialization::FSimaabEquiv) { + // :math:`U \sim U_d(\alpha, \alpha, \beta), \alpha \geq |\beta|` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l) \cdot Rz(\lambda_l)`. + auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, EulerBasis::ZYZ); + auto ab = (a_ + b_) / 2.; + + a_ = ab; + b_ = ab; + // unmodified parameter c + globalPhase_ = globalPhase_ + k2lphase; + k1l_ = k1l_ * rzMatrix(k2lphi); + k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); + k1r_ = k1r_ * rzMatrix(k2lphi); + k2r_ = rzMatrix(-k2lphi) * k2r_; + } else if (newSpecialization == Specialization::FSimabbEquiv) { + // :math:`U \sim U_d(\alpha, \beta, \beta), \alpha \geq \beta \geq 0` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l) \cdot Rx(\lambda_l)` + auto eulerBasis = EulerBasis::XYX; + auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, eulerBasis); + auto bc = (b_ + c_) / 2.; + + // unmodified parameter a + b_ = bc; + c_ = bc; + globalPhase_ = globalPhase_ + k2lphase; + k1l_ = k1l_ * rxMatrix(k2lphi); + k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); + k1r_ = k1r_ * rxMatrix(k2lphi); + k2r_ = rxMatrix(-k2lphi) * k2r_; + } else if (newSpecialization == Specialization::FSimabmbEquiv) { + // :math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` + // + // This gate binds 5 parameters, we make it canonical by setting: + // + // :math:`K2_l = Ry(\theta_l) \cdot Rx(\lambda_l)` + auto eulerBasis = EulerBasis::XYX; + auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + anglesFromUnitary(k2l_, eulerBasis); + auto bc = (b_ - c_) / 2.; + + // unmodified parameter a + 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_; + } else { + llvm::reportFatalInternalError( + "Unknown specialization for Weyl decomposition!"); + } + return flippedFromOriginal; +} + +TwoQubitBasisDecomposer +TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, + double basisFidelity) { + const Matrix2x2 k12RArr = Matrix2x2::fromElements( + 1i * INV_SQRT2, INV_SQRT2, -INV_SQRT2, -1i * INV_SQRT2); + const Matrix2x2 k12LArr = + Matrix2x2::fromElements(Complex{0.5, 0.5}, Complex{0.5, 0.5}, + Complex{-0.5, 0.5}, Complex{0.5, -0.5}); + + const auto basisDecomposer = + TwoQubitWeylDecomposition::create(basisMatrix, basisFidelity); + const auto isSuperControlled = + relativeEq(basisDecomposer.a(), std::numbers::pi / 4.0, 1e-13, 1e-09) && + relativeEq(basisDecomposer.c(), 0.0, 1e-13, 1e-09); + + // Create some useful matrices U1, U2, U3 are equivalent to the basis, + // expand as Ui = Ki1.Ubasis.Ki2 + auto b = basisDecomposer.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 k22l = + Matrix2x2::fromElements(INV_SQRT2, -INV_SQRT2, INV_SQRT2, INV_SQRT2); + const Matrix2x2 k22r = Matrix2x2::fromElements(0, 1, -1, 0); + 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))); + auto k1lDagger = basisDecomposer.k1l().adjoint(); + auto k1rDagger = basisDecomposer.k1r().adjoint(); + auto k2lDagger = basisDecomposer.k2l().adjoint(); + auto k2rDagger = basisDecomposer.k2r().adjoint(); + // Pre-build the fixed parts of the matrices used in 3-part decomposition + auto u0l = k31l * k1lDagger; + auto u0r = k31r * k1rDagger; + auto u1l = k2lDagger * k32lK21l * k1lDagger; + auto u1ra = k2rDagger * k32r; + auto u1rb = k21r * k1rDagger; + auto u2la = k2lDagger * k22l; + auto u2lb = k11l * k1lDagger; + auto u2ra = k2rDagger * k22r; + auto u2rb = k11r * k1rDagger; + auto u3l = k2lDagger * k12LArr; + auto u3r = k2rDagger * k12RArr; + // Pre-build the fixed parts of the matrices used in the 2-part decomposition + auto q0l = k12LArr.adjoint() * k1lDagger; + auto q0r = k12RArr.adjoint() * iPauliZ() * k1rDagger; + auto q1la = k2lDagger * k11l.adjoint(); + auto q1lb = k11l * k1lDagger; + auto q1ra = k2rDagger * iPauliZ() * k11r.adjoint(); + auto q1rb = k11r * k1rDagger; + auto q2l = k2lDagger * k12LArr; + auto q2r = k2rDagger * k12RArr; + + TwoQubitBasisDecomposer decomposer; + decomposer.basisFidelity = basisFidelity; + decomposer.basisDecomposer = basisDecomposer; + decomposer.isSuperControlled = isSuperControlled; + decomposer.smb.u0l = u0l; + decomposer.smb.u0r = u0r; + decomposer.smb.u1l = u1l; + decomposer.smb.u1ra = u1ra; + decomposer.smb.u1rb = u1rb; + decomposer.smb.u2la = u2la; + decomposer.smb.u2lb = u2lb; + decomposer.smb.u2ra = u2ra; + decomposer.smb.u2rb = u2rb; + decomposer.smb.u3l = u3l; + decomposer.smb.u3r = u3r; + decomposer.smb.q0l = q0l; + decomposer.smb.q0r = q0r; + decomposer.smb.q1la = q1la; + decomposer.smb.q1lb = q1lb; + decomposer.smb.q1ra = q1ra; + decomposer.smb.q1rb = q1rb; + decomposer.smb.q2l = q2l; + decomposer.smb.q2r = q2r; + return decomposer; +} + +std::optional +TwoQubitBasisDecomposer::twoQubitDecompose( + const TwoQubitWeylDecomposition& targetDecomposition, + std::optional numBasisGateUses) const { + auto traces = this->traces(targetDecomposition); + auto getDefaultNbasis = [&]() -> std::uint8_t { + // Pick the number of basis gate uses `i ∈ {0, 1, 2, 3}` that maximizes + // expected_fidelity(i) = traceToFidelity(traces[i]) * basisFidelity^i. + auto bestValue = std::numeric_limits::lowest(); + auto bestIndex = -1; + for (int i = 0; std::cmp_less(i, traces.size()); ++i) { + auto value = traceToFidelity(traces[i]) * std::pow(basisFidelity, i); + 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); + }; + // number of basis gates that need to be used in the decomposition + auto bestNbasis = numBasisGateUses.value_or(getDefaultNbasis()); + if (bestNbasis > 1 && !isSuperControlled) { + // cannot reliably decompose with more than one basis gate and a + // non-super-controlled basis gate + return std::nullopt; + } + auto chooseDecomposition = [&]() { + if (bestNbasis == 0) { + return decomp0(targetDecomposition); + } + if (bestNbasis == 1) { + return decomp1(targetDecomposition); + } + if (bestNbasis == 2) { + return decomp2Supercontrolled(targetDecomposition); + } + if (bestNbasis == 3) { + return decomp3Supercontrolled(targetDecomposition); + } + llvm::reportFatalInternalError( + "Invalid number of basis gates to use in basis decomposition (" + + llvm::Twine(bestNbasis) + ")!"); + llvm_unreachable(""); + }; + SmallVector factors = chooseDecomposition(); + + double globalPhase = targetDecomposition.globalPhase(); + globalPhase -= bestNbasis * basisDecomposer.globalPhase(); + if (bestNbasis == 2) { + // The two-basis (2x CX/CZ) template in `decomp2Supercontrolled` produces + // a sequence whose global phase is off by `pi` relative to the target; + // compensate here so the emitted sequence reproduces the target unitary + // exactly, not just up to sign. + globalPhase += std::numbers::pi; + } + // large global phases can be generated by the decomposition, thus limit + // it to [0, +2*pi) + globalPhase = remEuclid(globalPhase, 2.0 * std::numbers::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 { + // may not work for z != 0 and c != 0 (not always in Weyl chamber) + return SmallVector{ + basisDecomposer.k2r().adjoint() * target.k2r(), + basisDecomposer.k2l().adjoint() * target.k2l(), + target.k1r() * basisDecomposer.k1r().adjoint(), + target.k1l() * basisDecomposer.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 { + // Returns the Hilbert-Schmidt traces between the target canonical gate and + // the best candidate reachable with `0, 1, 2, 3` uses of the basis gate, + // respectively. Fed into `traceToFidelity` by `getDefaultNbasis` to pick + // the best basis-gate count. The closed-form expressions specialize + // `TwoQubitWeylDecomposition::getTrace(a, b, c, ap, bp, cp)` for: + // i == 0: no basis gate (ap == bp == cp == 0) + // i == 1: one basis use (ap == pi/4, bp == basis.b, cp == 0) + // i == 2: two basis uses (ap == 0, bp == 0, cp == -target.c) + // i == 3: three basis uses (target reachable exactly -> trace == 4) + // so the array has length 4 and is indexed by the number of basis uses. + 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((std::numbers::pi / 4.0) - target.a()) * + std::cos(basisDecomposer.b() - target.b()) * + std::cos(target.c()), + std::sin((std::numbers::pi / 4.0) - target.a()) * + std::sin(basisDecomposer.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); + const auto weyl = TwoQubitWeylDecomposition::create(target, std::nullopt); + return decomposer.twoQubitDecompose(weyl, numBasisUses); +} + +} // 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..454012d489 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp @@ -104,18 +104,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 = 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..d25df3c2da --- /dev/null +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -0,0 +1,256 @@ +/* + * 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 + +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(); }); +} + +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())); +} + +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; +}; + +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)); + } +} + +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()); + +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)); + } +} + +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)); +} + +INSTANTIATE_TEST_SUITE_P(ProductTwoQubitMatrices, BasisDecomposerTest, + testing::Combine(cxBasisCases(), + productMatrixCases())); +INSTANTIATE_TEST_SUITE_P(TwoQubitMatrices, BasisDecomposerTest, + testing::Combine(cxBasisCases(), + entangledMatrixCases())); From 2f5ef130541474e3524677a0df0b33f8af04a005 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 18:31:20 +0200 Subject: [PATCH 02/16] =?UTF-8?q?=E2=9C=A8=20Add=20specialized=20matrix=20?= =?UTF-8?q?test=20cases=20for=20Weyl=20decomposition=20and=20improve=20tol?= =?UTF-8?q?erance=20checks=20in=20existing=20tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Decomposition/test_weyl_decomposition.cpp | 49 +++++++++++++++++-- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp index d25df3c2da..8a5f34d17f 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -142,12 +142,50 @@ static auto cxBasisCases() { []() { 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)); + } +} + class WeylDecompositionTest : public testing::TestWithParam {}; class BasisDecomposerTest : public testing::TestWithParam< @@ -170,7 +208,8 @@ TEST_P(WeylDecompositionTest, ReconstructsWithinRequestedFidelity) { for (const double fidelity : {1.0, 1.0 - 1e-12}) { const auto decomposition = TwoQubitWeylDecomposition::create(originalMatrix, fidelity); - EXPECT_TRUE(restoreWeyl(decomposition).isApprox(originalMatrix)); + EXPECT_TRUE( + restoreWeyl(decomposition).isApprox(originalMatrix, WEYL_TOLERANCE)); } } @@ -205,6 +244,8 @@ 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}) { @@ -213,7 +254,8 @@ TEST_P(BasisDecomposerTest, ReconstructsWithinRequestedFidelity) { const auto decomposed = decomposer.twoQubitDecompose(*targetDecomposition, std::nullopt); ASSERT_TRUE(decomposed.has_value()); - EXPECT_TRUE(restoreBasis(*decomposed, basisMatrix).isApprox(target)); + EXPECT_TRUE(restoreBasis(*decomposed, basisMatrix) + .isApprox(target, WEYL_TOLERANCE)); } } @@ -245,7 +287,8 @@ TEST(BasisDecomposerNumBasisTest, ForcesZeroBasisUsesForIdentityTarget) { 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)); + EXPECT_TRUE( + restoreBasis(*decomposed, basis).isApprox(target, WEYL_TOLERANCE)); } INSTANTIATE_TEST_SUITE_P(ProductTwoQubitMatrices, BasisDecomposerTest, From 61ed24571cc367ab0fe74519be9e7d965dd3511c Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 18:36:44 +0200 Subject: [PATCH 03/16] =?UTF-8?q?=F0=9F=93=9D=20Simplify=20comments=20and?= =?UTF-8?q?=20docstrings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Weyl.h | 36 +--- mlir/include/mlir/Dialect/QCO/Utils/Matrix.h | 12 +- .../QCO/Transforms/Decomposition/Weyl.cpp | 165 ++---------------- 3 files changed, 31 insertions(+), 182 deletions(-) diff --git a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h index ac7ef42cee..9b59eafaae 100644 --- a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h @@ -57,27 +57,13 @@ class TwoQubitWeylDecomposition { [[nodiscard]] double c() const { return c_; } [[nodiscard]] double globalPhase() const { return globalPhase_; } - /** - * @brief Left single-qubit factor after the canonical gate. - * - * ``` - * q1 - k2r - C - k1r - - * A - * q0 - k2l - N - k1l - - * ``` - */ + /** @brief Single-qubit factor on qubit 0 after the canonical gate. */ [[nodiscard]] const Matrix2x2& k1l() const { return k1l_; } - /** - * @brief Left single-qubit factor before the canonical gate. - */ + /** @brief Single-qubit factor on qubit 0 before the canonical gate. */ [[nodiscard]] const Matrix2x2& k2l() const { return k2l_; } - /** - * @brief Right single-qubit factor after the canonical gate. - */ + /** @brief Single-qubit factor on qubit 1 after the canonical gate. */ [[nodiscard]] const Matrix2x2& k1r() const { return k1r_; } - /** - * @brief Right single-qubit factor before the canonical gate. - */ + /** @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, @@ -97,13 +83,11 @@ class TwoQubitWeylDecomposition { }; /** - * @brief Single-qubit factors and entangler count for basis-gate synthesis. + * @brief Native two-qubit synthesis result. * - * 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. + * 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; @@ -175,9 +159,7 @@ class TwoQubitBasisDecomposer { SmbPrecomputed smb{}; }; -/** - * @brief Decomposes @p target with respect to a two-qubit basis gate matrix. - */ +/** @brief Weyl-decomposes @p target using @p basisMatrix as entangler. */ [[nodiscard]] std::optional decomposeTwoQubitWithBasis( const Matrix4x4& target, const Matrix4x4& basisMatrix, diff --git a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h index 60324d3873..d9d46a714e 100644 --- a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h +++ b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h @@ -795,13 +795,13 @@ struct SymmetricEigen4 { Matrix4x4 eigenvectors{}; }; -/** @brief Constant-folded `RX(theta)` unitary. */ +/** @brief `RX(theta)` unitary matrix. */ [[nodiscard]] Matrix2x2 rxMatrix(double theta); -/** @brief Constant-folded `RY(theta)` unitary. */ +/** @brief `RY(theta)` unitary matrix. */ [[nodiscard]] Matrix2x2 ryMatrix(double theta); -/** @brief Constant-folded `RZ(theta)` unitary. */ +/** @brief `RZ(theta)` unitary matrix. */ [[nodiscard]] Matrix2x2 rzMatrix(double theta); /** @brief `i` times Pauli-X. */ @@ -813,13 +813,13 @@ struct SymmetricEigen4 { /** @brief `i` times Pauli-Z. */ [[nodiscard]] const Matrix2x2& iPauliZ(); -/** @brief Constant-folded two-qubit `RXX(theta)` unitary. */ +/** @brief `RXX(theta)` unitary matrix. */ [[nodiscard]] Matrix4x4 rxxMatrix(double theta); -/** @brief Constant-folded two-qubit `RYY(theta)` unitary. */ +/** @brief `RYY(theta)` unitary matrix. */ [[nodiscard]] Matrix4x4 ryyMatrix(double theta); -/** @brief Constant-folded two-qubit `RZZ(theta)` unitary. */ +/** @brief `RZZ(theta)` unitary matrix. */ [[nodiscard]] Matrix4x4 rzzMatrix(double theta); } // namespace mlir::qco diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index 8b47c02c87..397f114113 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -269,37 +269,23 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, std::optional fidelity) { auto u = unitaryMatrix; auto detU = u.determinant(); - // Project into SU(4) by dividing out the fourth root of det(U): for a 4x4 - // unitary, |det(U)| == 1 so `det^{-1/4}` both enforces det == 1 and removes - // the global phase. The extracted phase is tracked separately in - // `globalPhase` (quarter of arg(det) to match the fourth-root choice) so the - // caller can reconstruct the original matrix exactly if needed. + // Project into SU(4): det^{-1/4} removes global phase; arg(det)/4 is tracked. auto detPow = std::pow(detU, -0.25); - u *= detPow; // remove global phase from unitary matrix + u *= detPow; auto globalPhase = std::arg(detU) / 4.; - // Numerical drift can still leave tiny determinant errors after root - // normalization. Re-normalize once more instead of aborting. auto detNormalized = u.determinant(); if (std::abs(detNormalized - Complex{1.0, 0.0}) > WEYL_TOLERANCE && std::abs(detNormalized) > WEYL_TOLERANCE) { u *= std::pow(detNormalized, -0.25); } - // transform unitary matrix to magic basis; this enables two properties: - // 1. if uP ∈ SO(4), V = A ⊗ B (SO(4) → SU(2) ⊗ SU(2)) - // 2. magic basis diagonalizes canonical gate, allowing calculation of - // canonical gate parameters later on auto uP = magicBasisTransform(u, /*outOfMagicBasis=*/true); const Matrix4x4 m2 = uP.transpose() * uP; - // diagonalization yields eigenvectors (p) and eigenvalues (d); - // p is used to calculate K1/K2 (and thus the single-qubit gates - // surrounding the canonical gate); d is used to determine the Weyl - // coordinates and thus the parameters of the canonical gate auto [p, d] = diagonalizeComplexSymmetric(m2, DIAGONALIZATION_PRECISION); - // extract Weyl coordinates from eigenvalues, map to [0, 2*pi) + // Map eigenvalue phases to Weyl coordinates in [0, 2*pi). constexpr double pi = std::numbers::pi; std::array dReal{}; for (std::size_t i = 0; i < d.size(); ++i) { @@ -311,8 +297,7 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, cs[i] = remEuclid((dReal[i] + dReal[3]) / 2.0, 2.0 * pi); } - // Reorder coordinates according to min(a, pi/2 - a) with - // a = x mod pi/2 for each Weyl coordinate x + // 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); @@ -331,13 +316,10 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, dReal[2] = reordered[2]; } - // update eigenvectors (columns of p) according to new order of - // weyl coordinates const Matrix4x4 pOrig = p; for (std::size_t i = 0; i < order.size(); ++i) { p.setColumn(i, pOrig.column(order[i])); } - // apply correction for determinant if necessary if (p.determinant().real() < 0.0) { auto lastColumn = p.column(3); for (auto& entry : lastColumn) { @@ -347,33 +329,23 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, } assert(std::abs(p.determinant() - 1.0) < WEYL_TOLERANCE); - // re-create complex eigenvalue matrix; this matrix contains the - // parameters of the canonical gate which is later used in the - // verification. Since the matrix is diagonal, the matrix exponential is - // equivalent to the element-wise exponential function. 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); - // combined matrix k1 of 1q gates after canonical gate Matrix4x4 k1 = uP * p * temp; - // k1 must be orthogonal; the tolerance matches the iterative diagonalization - // residual rather than the (much tighter) default matrix tolerance. + // 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); - // combined matrix k2 of 1q gates before canonical gate Matrix4x4 k2 = p.adjoint(); - // k2 must be orthogonal; see the tolerance note on the k1 check above. assert((k2.transpose() * k2).isIdentity(WEYL_TOLERANCE)); assert(k2.determinant().real() > 0.0); k2 = magicBasisTransform(k2, /*outOfMagicBasis=*/false); - // ensure k1 and k2 are correct (when combined with the canonical gate - // parameters in-between, they are equivalent to u) std::array tempConjDiag{}; for (std::size_t k = 0; k < tempConjDiag.size(); ++k) { tempConjDiag[k] = std::conj(tempDiag[k]); @@ -384,16 +356,13 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, k2) .isApprox(u, WEYL_TOLERANCE)); - // calculate k1 = K1l ⊗ K1r auto [K1l, K1r, phaseL] = decomposeTwoQubitProductGate(k1); - // decompose k2 = K2l ⊗ K2r auto [K2l, K2r, phaseR] = decomposeTwoQubitProductGate(k2); assert(Matrix4x4::kron(K1l, K1r).isApprox(k1, WEYL_TOLERANCE)); assert(Matrix4x4::kron(K2l, K2r).isApprox(k2, WEYL_TOLERANCE)); - // accumulate global phase globalPhase += phaseL + phaseR; - // Flip into Weyl chamber + // Map into the Weyl chamber. if (cs[0] > (pi / 2.0)) { cs[0] -= 3.0 * (pi / 2.0); K1l = K1l * iPauliY(); @@ -446,7 +415,7 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, globalPhase -= (pi / 2.0); } - // bind weyl coordinates as parameters of canonical gate + // Bind Weyl coordinates to the canonical gate. auto [a, b, c] = std::tie(cs[1], cs[0], cs[2]); TwoQubitWeylDecomposition decomposition; @@ -459,13 +428,10 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, decomposition.k1r_ = K1r; decomposition.k2r_ = K2r; - // make sure decomposition is equal to input assert((Matrix4x4::kron(K1l, K1r) * decomposition.getCanonicalMatrix() * Matrix4x4::kron(K2l, K2r) * std::exp(Complex{0.0, 1.0} * globalPhase)) .isApprox(unitaryMatrix, WEYL_TOLERANCE)); - // determine actual specialization of canonical gate so that the 1q - // matrices can potentially be simplified auto flippedFromOriginal = decomposition.applySpecialization(fidelity); auto getTraceValue = [&]() { @@ -476,13 +442,8 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, return getTrace(a, b, c, decomposition.a_, decomposition.b_, decomposition.c_); }; - // use trace to calculate fidelity of applied specialization and - // adjust global phase - auto trace = getTraceValue(); + const auto trace = getTraceValue(); const double calculatedFidelity = traceToFidelity(trace); - // final check if specialization is close enough to the original matrix to - // satisfy the requested fidelity; since no forced specialization is - // allowed, this should never fail if (fidelity && calculatedFidelity + 1.0e-13 < *fidelity) { llvm::reportFatalInternalError(llvm::formatv( "TwoQubitWeylDecomposition: Calculated fidelity of " @@ -491,7 +452,6 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, } decomposition.globalPhase_ += std::arg(trace); - // final check if decomposition is still valid after specialization assert((Matrix4x4::kron(decomposition.k1l_, decomposition.k1r_) * decomposition.getCanonicalMatrix() * Matrix4x4::kron(decomposition.k2l_, decomposition.k2r_) * @@ -503,12 +463,8 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, Matrix4x4 TwoQubitWeylDecomposition::getCanonicalMatrix(double a, double b, double c) { - // Canonical gate `U_d(a, b, c) = exp(i * (a*XX + b*YY + c*ZZ))`. XX/YY/ZZ - // commute pairwise, so any product order is equivalent; the order below is - // chosen to match common Qiskit/QuantumFlow references. The negated rotation - // angles (`-2 * a`, ...) compensate for the `RXX/RYY/RZZ` convention - // `exp(-i * theta/2 * XX)`, so that the factored angles sum back to the - // intended `+a`, `+b`, `+c`. + // U_d(a,b,c) = exp(i(a XX + b YY + c ZZ)); `-2*a/b/c` matches RXX/RYY/RZZ + // gate convention; product order matches Qiskit (RZZ · RYY · RXX). const auto xx = rxxMatrix(-2.0 * a); const auto yy = ryyMatrix(-2.0 * b); const auto zz = rzzMatrix(-2.0 * c); @@ -520,37 +476,19 @@ bool TwoQubitWeylDecomposition::applySpecialization( bool flippedFromOriginal = false; const auto newSpecialization = bestSpecialization(*this, requestedFidelity); if (newSpecialization == Specialization::General) { - // U has no special symmetry. - // - // This gate binds all 6 possible parameters, so there is no need to - // make the single-qubit pre-/post-gates canonical. return flippedFromOriginal; } if (newSpecialization == Specialization::IdEquiv) { - // :math:`U \sim U_d(0,0,0)` - // Thus, :math:`\sim Id` - // - // This gate binds 0 parameters, we make it canonical by setting: - // - // :math:`K2_l = Id` , :math:`K2_r = Id`. a_ = 0.; b_ = 0.; c_ = 0.; - // unmodified global phase k1l_ = k1l_ * k2l_; k2l_ = Matrix2x2::identity(); k1r_ = k1r_ * k2r_; k2r_ = Matrix2x2::identity(); } else if (newSpecialization == Specialization::SWAPEquiv) { - // :math:`U \sim U_d(\pi/4, \pi/4, \pi/4) \sim U(\pi/4, \pi/4, -\pi/4)` - // Thus, :math:`U \sim \text{SWAP}` - // - // This gate binds 0 parameters, we make it canonical by setting: - // - // :math:`K2_l = Id` , :math:`K2_r = Id`. if (c_ > 0.) { - // unmodified global phase k1l_ = k1l_ * k2r_; k1r_ = k1r_ * k2l_; k2l_ = Matrix2x2::identity(); @@ -568,58 +506,33 @@ bool TwoQubitWeylDecomposition::applySpecialization( b_ = (std::numbers::pi / 4.0); c_ = (std::numbers::pi / 4.0); } else if (newSpecialization == Specialization::PartialSWAPEquiv) { - // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, \alpha\pi/4)` - // Thus, :math:`U \sim \text{SWAP}^\alpha` - // - // This gate binds 3 parameters, we make it canonical by setting: - // - // :math:`K2_l = Id`. auto closest = closestPartialSwap(a_, b_, c_); auto k2lDagger = k2l_.adjoint(); a_ = closest; b_ = closest; c_ = closest; - // unmodified global phase k1l_ = k1l_ * k2l_; k1r_ = k1r_ * k2l_; k2r_ = k2lDagger * k2r_; k2l_ = Matrix2x2::identity(); } else if (newSpecialization == Specialization::PartialSWAPFlipEquiv) { - // :math:`U \sim U_d(\alpha\pi/4, \alpha\pi/4, -\alpha\pi/4)` - // Thus, :math:`U \sim \text{SWAP}^\alpha` - // - // (a non-equivalent root of SWAP from the TwoQubitWeylPartialSWAPEquiv - // similar to how :math:`x = (\pm \sqrt(x))^2`) - // - // This gate binds 3 parameters, we make it canonical by setting: - // - // :math:`K2_l = Id` auto closest = closestPartialSwap(a_, b_, -c_); auto k2lDagger = k2l_.adjoint(); a_ = closest; b_ = closest; c_ = -closest; - // unmodified global phase k1l_ = k1l_ * k2l_; k1r_ = k1r_ * iPauliZ() * k2l_ * iPauliZ(); k2r_ = iPauliZ() * k2lDagger * iPauliZ() * k2r_; k2l_ = Matrix2x2::identity(); } else if (newSpecialization == Specialization::ControlledEquiv) { - // :math:`U \sim U_d(\alpha, 0, 0)` - // Thus, :math:`U \sim \text{Ctrl-U}` - // - // This gate binds 4 parameters, we make it canonical by setting: - // - // :math:`K2_l = Ry(\theta_l) Rx(\lambda_l)` - // :math:`K2_r = Ry(\theta_r) Rx(\lambda_r)` const EulerBasis eulerBasis = EulerBasis::XYX; const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = anglesFromUnitary(k2l_, eulerBasis); const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = anglesFromUnitary(k2r_, EulerBasis::XYX); - // unmodified parameter a b_ = 0.; c_ = 0.; globalPhase_ = globalPhase_ + k2lphase + k2rphase; @@ -628,55 +541,35 @@ bool TwoQubitWeylDecomposition::applySpecialization( k1r_ = k1r_ * rxMatrix(k2rphi); k2r_ = ryMatrix(k2rtheta) * rxMatrix(k2rlambda); } else if (newSpecialization == Specialization::MirrorControlledEquiv) { - // :math:`U \sim U_d(\pi/4, \pi/4, \alpha)` - // Thus, :math:`U \sim \text{SWAP} \cdot \text{Ctrl-U}` - // - // This gate binds 4 parameters, we make it canonical by setting: - // - // :math:`K2_l = Ry(\theta_l)\cdot Rz(\lambda_l)` - // :math:`K2_r = Ry(\theta_r)\cdot Rz(\lambda_r)` const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = anglesFromUnitary(k2l_, EulerBasis::ZYZ); const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = anglesFromUnitary(k2r_, EulerBasis::ZYZ); a_ = (std::numbers::pi / 4.0); b_ = (std::numbers::pi / 4.0); - // unmodified parameter c globalPhase_ = globalPhase_ + k2lphase + k2rphase; k1l_ = k1l_ * rzMatrix(k2rphi); k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); k1r_ = k1r_ * rzMatrix(k2lphi); k2r_ = ryMatrix(k2rtheta) * rzMatrix(k2rlambda); } else if (newSpecialization == Specialization::FSimaabEquiv) { - // :math:`U \sim U_d(\alpha, \alpha, \beta), \alpha \geq |\beta|` - // - // This gate binds 5 parameters, we make it canonical by setting: - // - // :math:`K2_l = Ry(\theta_l) \cdot Rz(\lambda_l)`. auto [k2ltheta, k2lphi, k2llambda, k2lphase] = anglesFromUnitary(k2l_, EulerBasis::ZYZ); auto ab = (a_ + b_) / 2.; a_ = ab; b_ = ab; - // unmodified parameter c globalPhase_ = globalPhase_ + k2lphase; k1l_ = k1l_ * rzMatrix(k2lphi); k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); k1r_ = k1r_ * rzMatrix(k2lphi); k2r_ = rzMatrix(-k2lphi) * k2r_; } else if (newSpecialization == Specialization::FSimabbEquiv) { - // :math:`U \sim U_d(\alpha, \beta, \beta), \alpha \geq \beta \geq 0` - // - // This gate binds 5 parameters, we make it canonical by setting: - // - // :math:`K2_l = Ry(\theta_l) \cdot Rx(\lambda_l)` auto eulerBasis = EulerBasis::XYX; auto [k2ltheta, k2lphi, k2llambda, k2lphase] = anglesFromUnitary(k2l_, eulerBasis); auto bc = (b_ + c_) / 2.; - // unmodified parameter a b_ = bc; c_ = bc; globalPhase_ = globalPhase_ + k2lphase; @@ -685,17 +578,11 @@ bool TwoQubitWeylDecomposition::applySpecialization( k1r_ = k1r_ * rxMatrix(k2lphi); k2r_ = rxMatrix(-k2lphi) * k2r_; } else if (newSpecialization == Specialization::FSimabmbEquiv) { - // :math:`U \sim U_d(\alpha, \beta, -\beta), \alpha \geq \beta \geq 0` - // - // This gate binds 5 parameters, we make it canonical by setting: - // - // :math:`K2_l = Ry(\theta_l) \cdot Rx(\lambda_l)` auto eulerBasis = EulerBasis::XYX; auto [k2ltheta, k2lphi, k2llambda, k2lphase] = anglesFromUnitary(k2l_, eulerBasis); auto bc = (b_ - c_) / 2.; - // unmodified parameter a b_ = bc; c_ = -bc; globalPhase_ = globalPhase_ + k2lphase; @@ -725,8 +612,6 @@ TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, relativeEq(basisDecomposer.a(), std::numbers::pi / 4.0, 1e-13, 1e-09) && relativeEq(basisDecomposer.c(), 0.0, 1e-13, 1e-09); - // Create some useful matrices U1, U2, U3 are equivalent to the basis, - // expand as Ui = Ki1.Ubasis.Ki2 auto b = basisDecomposer.b(); Complex temp{0.5, -0.5}; const Matrix2x2 k11l = Matrix2x2::fromElements( @@ -758,7 +643,6 @@ TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, auto k1rDagger = basisDecomposer.k1r().adjoint(); auto k2lDagger = basisDecomposer.k2l().adjoint(); auto k2rDagger = basisDecomposer.k2r().adjoint(); - // Pre-build the fixed parts of the matrices used in 3-part decomposition auto u0l = k31l * k1lDagger; auto u0r = k31r * k1rDagger; auto u1l = k2lDagger * k32lK21l * k1lDagger; @@ -770,7 +654,6 @@ TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, auto u2rb = k11r * k1rDagger; auto u3l = k2lDagger * k12LArr; auto u3r = k2rDagger * k12RArr; - // Pre-build the fixed parts of the matrices used in the 2-part decomposition auto q0l = k12LArr.adjoint() * k1lDagger; auto q0r = k12RArr.adjoint() * iPauliZ() * k1rDagger; auto q1la = k2lDagger * k11l.adjoint(); @@ -812,8 +695,7 @@ TwoQubitBasisDecomposer::twoQubitDecompose( std::optional numBasisGateUses) const { auto traces = this->traces(targetDecomposition); auto getDefaultNbasis = [&]() -> std::uint8_t { - // Pick the number of basis gate uses `i ∈ {0, 1, 2, 3}` that maximizes - // expected_fidelity(i) = traceToFidelity(traces[i]) * basisFidelity^i. + // Maximize traceToFidelity(traces[i]) * basisFidelity^i over i in {0..3}. auto bestValue = std::numeric_limits::lowest(); auto bestIndex = -1; for (int i = 0; std::cmp_less(i, traces.size()); ++i) { @@ -832,11 +714,9 @@ TwoQubitBasisDecomposer::twoQubitDecompose( } return static_cast(bestIndex); }; - // number of basis gates that need to be used in the decomposition auto bestNbasis = numBasisGateUses.value_or(getDefaultNbasis()); if (bestNbasis > 1 && !isSuperControlled) { - // cannot reliably decompose with more than one basis gate and a - // non-super-controlled basis gate + // More than one basis gate requires a super-controlled basis (e.g. CX). return std::nullopt; } auto chooseDecomposition = [&]() { @@ -862,14 +742,9 @@ TwoQubitBasisDecomposer::twoQubitDecompose( double globalPhase = targetDecomposition.globalPhase(); globalPhase -= bestNbasis * basisDecomposer.globalPhase(); if (bestNbasis == 2) { - // The two-basis (2x CX/CZ) template in `decomp2Supercontrolled` produces - // a sequence whose global phase is off by `pi` relative to the target; - // compensate here so the emitted sequence reproduces the target unitary - // exactly, not just up to sign. + // Two-basis template is exact up to a pi global-phase offset. globalPhase += std::numbers::pi; } - // large global phases can be generated by the decomposition, thus limit - // it to [0, +2*pi) globalPhase = remEuclid(globalPhase, 2.0 * std::numbers::pi); return TwoQubitNativeDecomposition{ @@ -889,7 +764,7 @@ TwoQubitBasisDecomposer::decomp0(const TwoQubitWeylDecomposition& target) { SmallVector TwoQubitBasisDecomposer::decomp1( const TwoQubitWeylDecomposition& target) const { - // may not work for z != 0 and c != 0 (not always in Weyl chamber) + // One basis gate; reliable only when the target is in the Weyl chamber. return SmallVector{ basisDecomposer.k2r().adjoint() * target.k2r(), basisDecomposer.k2l().adjoint() * target.k2l(), @@ -936,16 +811,8 @@ SmallVector TwoQubitBasisDecomposer::decomp3Supercontrolled( std::array, 4> TwoQubitBasisDecomposer::traces(const TwoQubitWeylDecomposition& target) const { - // Returns the Hilbert-Schmidt traces between the target canonical gate and - // the best candidate reachable with `0, 1, 2, 3` uses of the basis gate, - // respectively. Fed into `traceToFidelity` by `getDefaultNbasis` to pick - // the best basis-gate count. The closed-form expressions specialize - // `TwoQubitWeylDecomposition::getTrace(a, b, c, ap, bp, cp)` for: - // i == 0: no basis gate (ap == bp == cp == 0) - // i == 1: one basis use (ap == pi/4, bp == basis.b, cp == 0) - // i == 2: two basis uses (ap == 0, bp == 0, cp == -target.c) - // i == 3: three basis uses (target reachable exactly -> trace == 4) - // so the array has length 4 and is indexed by the number of basis uses. + // Hilbert-Schmidt traces for 0..3 basis uses; index i == numBasisUses. + // i=0: identity; i=1: one basis gate; i=2: two uses; i=3: exact (trace=4). return { 4. * std::complex{std::cos(target.a()) * std::cos(target.b()) * std::cos(target.c()), From 5d27f045550474c900eae234e69009ca16f1299e Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 22 Jun 2026 22:02:36 +0200 Subject: [PATCH 04/16] =?UTF-8?q?=F0=9F=9A=A8=20Fix=20linter=20warnings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/test_weyl_decomposition.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp index 8a5f34d17f..8a60d7a32b 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -186,6 +186,8 @@ TEST(DecompositionHelpersTest, GateMatrixFactoriesMatchCanonicalForm) { } } +namespace { + class WeylDecompositionTest : public testing::TestWithParam {}; class BasisDecomposerTest : public testing::TestWithParam< @@ -203,6 +205,8 @@ class BasisDecomposerTest : public testing::TestWithParam< std::unique_ptr targetDecomposition; }; +} // namespace + TEST_P(WeylDecompositionTest, ReconstructsWithinRequestedFidelity) { const Matrix4x4 originalMatrix = GetParam()(); for (const double fidelity : {1.0, 1.0 - 1e-12}) { From d57dab29dda1af607aca750028667b4bb430fef6 Mon Sep 17 00:00:00 2001 From: Simon Date: Mon, 22 Jun 2026 22:02:54 +0200 Subject: [PATCH 05/16] =?UTF-8?q?=F0=9F=93=9D=20Add=20comment=20about=20fi?= =?UTF-8?q?xed=20perturbation=20coefficients=20for=20deterministic=20behav?= =?UTF-8?q?ior=20in=20Weyl=20decomposition=20diagonalization?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index 397f114113..5aef8af884 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -109,6 +109,10 @@ diagonalizeComplexSymmetric(const Matrix4x4& m, double precision) { 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; From eaa10973998f924c6e3abc7d4bc7cfcd21b7265e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 22 Jun 2026 20:04:44 +0000 Subject: [PATCH 06/16] =?UTF-8?q?=F0=9F=8E=A8=20pre-commit=20fixes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index 5aef8af884..1a042db1b2 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -112,7 +112,8 @@ diagonalizeComplexSymmetric(const Matrix4x4& m, double precision) { // 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. + // keeps behavior deterministic while later attempts sample the + // distribution. if (i == 0) { randA = 1.2602066112249388; randB = 0.22317849046722027; From 362692c29c9e58cd6ca263346233a22ad7739073 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 12:52:46 +0200 Subject: [PATCH 07/16] =?UTF-8?q?=E2=9C=A8=20Add=20project=20options=20for?= =?UTF-8?q?=20MLIRQCOTransforms=20in=20CMake=20configuration?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Transforms/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) 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) From 425a03d21a85b39f1b4525bf2a83d31d8fb4a389 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 13:34:11 +0200 Subject: [PATCH 08/16] =?UTF-8?q?=E2=98=82=EF=B8=8F=20Increase=20coverage?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Decomposition/test_weyl_decomposition.cpp | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp index 8a60d7a32b..075724d040 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -295,6 +296,80 @@ TEST(BasisDecomposerNumBasisTest, ForcesZeroBasisUsesForIdentityTarget) { 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, 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(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())); From 426f3002479b745f24248c94d410fed677b8e0ff Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 13:48:44 +0200 Subject: [PATCH 09/16] =?UTF-8?q?=E2=98=82=EF=B8=8F=20Increase=20coverage?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../test_euler_decomposition.cpp | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) 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 454012d489..da86f990d9 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 @@ -441,6 +442,80 @@ 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(EulerAnglesDeathTest, InvalidBasisIsUnreachable) { + EXPECT_DEATH(static_cast(anglesFromUnitary( + Matrix2x2::identity(), static_cast(99))), + "invalid Euler basis"); +} + +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 //===----------------------------------------------------------------------===// From 392f6086a13fd53b555dc6e5f26c80733bdc8be4 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 14:36:34 +0200 Subject: [PATCH 10/16] =?UTF-8?q?=F0=9F=8E=A8=20Align=20implementation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Weyl.cpp | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index 1a042db1b2..d4459e2788 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -52,7 +52,11 @@ enum class Specialization : std::uint8_t { } // namespace -static constexpr auto DIAGONALIZATION_PRECISION = WEYL_TOLERANCE / 10; +/** Default fidelity for Weyl specialization. */ +static constexpr double DEFAULT_FIDELITY = 1.0 - 1e-12; + +/** M2 diagonalization tolerance. */ +static constexpr double DIAGONALIZATION_PRECISION = 1e-13; /** Non-negative remainder of @p a modulo @p b (always in `[0, |b|)`). */ static double remEuclid(const double a, const double b) { @@ -279,12 +283,6 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, u *= detPow; auto globalPhase = std::arg(detU) / 4.; - auto detNormalized = u.determinant(); - if (std::abs(detNormalized - Complex{1.0, 0.0}) > WEYL_TOLERANCE && - std::abs(detNormalized) > WEYL_TOLERANCE) { - u *= std::pow(detNormalized, -0.25); - } - auto uP = magicBasisTransform(u, /*outOfMagicBasis=*/true); const Matrix4x4 m2 = uP.transpose() * uP; @@ -612,7 +610,7 @@ TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, Complex{-0.5, 0.5}, Complex{0.5, -0.5}); const auto basisDecomposer = - TwoQubitWeylDecomposition::create(basisMatrix, basisFidelity); + TwoQubitWeylDecomposition::create(basisMatrix, DEFAULT_FIDELITY); const auto isSuperControlled = relativeEq(basisDecomposer.a(), std::numbers::pi / 4.0, 1e-13, 1e-09) && relativeEq(basisDecomposer.c(), 0.0, 1e-13, 1e-09); @@ -842,7 +840,7 @@ decomposeTwoQubitWithBasis(const Matrix4x4& target, const std::optional numBasisUses) { const auto decomposer = TwoQubitBasisDecomposer::create(basisMatrix, basisFidelity); - const auto weyl = TwoQubitWeylDecomposition::create(target, std::nullopt); + const auto weyl = TwoQubitWeylDecomposition::create(target, DEFAULT_FIDELITY); return decomposer.twoQubitDecompose(weyl, numBasisUses); } From 65408d6d0222437c60b705dbf3ba09b33c7bcbc7 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 14:42:42 +0200 Subject: [PATCH 11/16] =?UTF-8?q?=F0=9F=9A=A8=20Fix=20linter=20warnings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Transforms/Decomposition/test_euler_decomposition.cpp | 6 ------ 1 file changed, 6 deletions(-) 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 da86f990d9..aedc27d80f 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_euler_decomposition.cpp @@ -452,12 +452,6 @@ TEST(EulerAnglesCoverageTest, ParamsZYZUsesOffDiagonal01When10IsNearZero) { EXPECT_TRUE(std::isfinite(angles.lambda)); } -TEST(EulerAnglesDeathTest, InvalidBasisIsUnreachable) { - EXPECT_DEATH(static_cast(anglesFromUnitary( - Matrix2x2::identity(), static_cast(99))), - "invalid Euler basis"); -} - TEST(EulerAnglesCoverageTest, PhaseOnlyDecompositionSkipsRotationGates) { TestFixture fx; fx.setUp(); From 846a4edcfee77522cc963f9adc7604afecda5ff5 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 15:11:15 +0200 Subject: [PATCH 12/16] =?UTF-8?q?=E2=9A=A1=EF=B8=8F=20Enhance=20TwoQubitBa?= =?UTF-8?q?sisDecomposer=20with=20caching=20functionality=20for=20target?= =?UTF-8?q?=20decompositions.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Weyl.h | 30 +++++++++++++-- .../QCO/Transforms/Decomposition/Weyl.cpp | 38 ++++++++++++------- .../Decomposition/test_weyl_decomposition.cpp | 23 +++++++++++ 3 files changed, 74 insertions(+), 17 deletions(-) diff --git a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h index 9b59eafaae..3b51ce902f 100644 --- a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h @@ -32,7 +32,7 @@ inline constexpr double WEYL_TOLERANCE = 100 * MATRIX_TOLERANCE; * `U_canon(a,b,c) = RXX(-2a) · RYY(-2b) · RZZ(-2c)`. * * @note Adapted from TwoQubitWeylDecomposition in the IBM Qiskit framework. - * (C) Copyright IBM 2023 + * (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 @@ -99,7 +99,7 @@ struct TwoQubitNativeDecomposition { * @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 + * (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 @@ -112,6 +112,14 @@ struct TwoQubitNativeDecomposition { */ 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); @@ -119,6 +127,16 @@ class TwoQubitBasisDecomposer { 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; @@ -159,7 +177,13 @@ class TwoQubitBasisDecomposer { SmbPrecomputed smb{}; }; -/** @brief Weyl-decomposes @p target using @p basisMatrix as entangler. */ +/** + * @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, diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index d4459e2788..f2c09d181d 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -75,22 +75,24 @@ static double traceToFidelity(const Complex& trace) { static constexpr double INV_SQRT2 = 1.0 / std::numbers::sqrt2; +static const Matrix4x4 kMagicBasisNonNormalized = Matrix4x4::fromElements( // + 1, 1i, 0, 0, // + 0, 0, 1i, 1, // + 0, 0, 1i, -1, // + 1, -1i, 0, 0); +static const Matrix4x4 kMagicBasisNonNormalizedDagger = + 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) { - const Matrix4x4 bNonNormalized = Matrix4x4::fromElements( // - 1, 1i, 0, 0, // - 0, 0, 1i, 1, // - 0, 0, 1i, -1, // - 1, -1i, 0, 0); - const Matrix4x4 bNonNormalizedDagger = 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); if (outOfMagicBasis) { - return bNonNormalizedDagger * unitary * bNonNormalized; + return kMagicBasisNonNormalizedDagger * unitary * kMagicBasisNonNormalized; } - return bNonNormalized * unitary * bNonNormalizedDagger; + return kMagicBasisNonNormalized * unitary * kMagicBasisNonNormalizedDagger; } static double closestPartialSwap(double a, double b, double c) { @@ -692,6 +694,15 @@ TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, return decomposer; } +std::optional +TwoQubitBasisDecomposer::decomposeTarget( + const Matrix4x4& targetUnitary, + const std::optional numBasisGateUses) const { + const auto targetWeyl = + TwoQubitWeylDecomposition::create(targetUnitary, DEFAULT_FIDELITY); + return twoQubitDecompose(targetWeyl, numBasisGateUses); +} + std::optional TwoQubitBasisDecomposer::twoQubitDecompose( const TwoQubitWeylDecomposition& targetDecomposition, @@ -840,8 +851,7 @@ decomposeTwoQubitWithBasis(const Matrix4x4& target, const std::optional numBasisUses) { const auto decomposer = TwoQubitBasisDecomposer::create(basisMatrix, basisFidelity); - const auto weyl = TwoQubitWeylDecomposition::create(target, DEFAULT_FIDELITY); - return decomposer.twoQubitDecompose(weyl, numBasisUses); + return decomposer.decomposeTarget(target, numBasisUses); } } // namespace mlir::qco::decomposition diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp index 075724d040..bcdc627c60 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -308,6 +308,29 @@ TEST(BasisDecomposerTest, DecomposeTwoQubitWithBasisReconstructsTarget) { 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); From 7e4e69acad077b080b5b3408c9e16210c5235de5 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 15:48:23 +0200 Subject: [PATCH 13/16] =?UTF-8?q?=E2=9C=A8=20Move=20basis=20decomposition?= =?UTF-8?q?=20into=20TwoQubitBasisDecomposer=20and=20enhance=20Weyl=20deco?= =?UTF-8?q?mposition=20with=20new=20specialization=20phase=20finalization?= =?UTF-8?q?=20and=20improved=20fidelity=20calculations.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Weyl.h | 8 +- .../Decomposition/BasisDecomposer.cpp | 307 +++++++++ .../QCO/Transforms/Decomposition/Weyl.cpp | 633 ++++++------------ 3 files changed, 512 insertions(+), 436 deletions(-) create mode 100644 mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp diff --git a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h index 3b51ce902f..5a0483628f 100644 --- a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h @@ -72,6 +72,12 @@ class TwoQubitWeylDecomposition { 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_{}; @@ -172,7 +178,7 @@ class TwoQubitBasisDecomposer { traces(const TwoQubitWeylDecomposition& target) const; double basisFidelity{}; - TwoQubitWeylDecomposition basisDecomposer; + TwoQubitWeylDecomposition basisWeyl; bool isSuperControlled{}; SmbPrecomputed smb{}; }; 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..30c5f5709b --- /dev/null +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp @@ -0,0 +1,307 @@ +/* + * 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) { + 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: + llvm::reportFatalInternalError(llvm::formatv( + "Invalid number of basis gates to use in basis decomposition ({0})!", + bestNbasis)); + llvm_unreachable(""); + } + + 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/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index f2c09d181d..0ba70c6b9a 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -24,7 +24,6 @@ #include #include #include -#include #include #include #include @@ -50,15 +49,10 @@ enum class Specialization : std::uint8_t { FSimabmbEquiv, }; -} // namespace - -/** Default fidelity for Weyl specialization. */ -static constexpr double DEFAULT_FIDELITY = 1.0 - 1e-12; +static constexpr double kDiagonalizationPrecision = 1e-13; +static constexpr double kPi = std::numbers::pi; +static constexpr double kPiOver4 = kPi / 4.0; -/** M2 diagonalization tolerance. */ -static constexpr double DIAGONALIZATION_PRECISION = 1e-13; - -/** Non-negative remainder of @p a modulo @p b (always in `[0, |b|)`). */ static double remEuclid(const double a, const double b) { if (b == 0.0) { llvm::reportFatalInternalError("remEuclid expects non-zero divisor"); @@ -73,8 +67,6 @@ static double traceToFidelity(const Complex& trace) { return (4.0 + (traceAbs * traceAbs)) / 20.0; } -static constexpr double INV_SQRT2 = 1.0 / std::numbers::sqrt2; - static const Matrix4x4 kMagicBasisNonNormalized = Matrix4x4::fromElements( // 1, 1i, 0, 0, // 0, 0, 1i, 1, // @@ -96,9 +88,9 @@ static Matrix4x4 magicBasisTransform(const Matrix4x4& unitary, } static double closestPartialSwap(double a, double b, double c) { - auto m = (a + b + c) / 3.; - auto [am, bm, cm] = std::array{a - m, b - m, c - m}; - auto [ab, bc, ca] = std::array{a - b, b - c, c - a}; + 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.); } @@ -175,7 +167,7 @@ decomposeTwoQubitProductGate(const Matrix4x4& specialUnitary) { r *= (1.0 / std::sqrt(detR)); const Matrix2x2 rTConj = r.adjoint(); - Matrix4x4 temp = + const Matrix4x4 temp = specialUnitary * Matrix4x4::kron(Matrix2x2::identity(), rTConj); Matrix2x2 l = @@ -186,16 +178,16 @@ decomposeTwoQubitProductGate(const Matrix4x4& specialUnitary) { "decomposeTwoQubitProductGate: unable to decompose: detL < 0.9"); } l *= (1.0 / std::sqrt(detL)); - auto phase = std::arg(detL) / 2.; + 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) { - auto da = a - ap; - auto db = b - bp; - auto dc = c - 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)}; } @@ -204,26 +196,24 @@ static Specialization bestSpecialization(const TwoQubitWeylDecomposition& decomposition, const std::optional& requestedFidelity) { auto isClose = [&](double ap, double bp, double cp) -> bool { - auto tr = getTrace(decomposition.a(), decomposition.b(), decomposition.c(), - ap, bp, cp); + const auto tr = getTrace(decomposition.a(), decomposition.b(), + decomposition.c(), ap, bp, cp); if (requestedFidelity) { return traceToFidelity(tr) >= *requestedFidelity; } return false; }; - auto closestAbc = closestPartialSwap(decomposition.a(), decomposition.b(), - decomposition.c()); - auto closestAbMinusC = closestPartialSwap( + 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((std::numbers::pi / 4.0), (std::numbers::pi / 4.0), - (std::numbers::pi / 4.0)) || - isClose((std::numbers::pi / 4.0), (std::numbers::pi / 4.0), - -(std::numbers::pi / 4.0))) { + if (isClose(kPiOver4, kPiOver4, kPiOver4) || + isClose(kPiOver4, kPiOver4, -kPiOver4)) { return Specialization::SWAPEquiv; } if (isClose(closestAbc, closestAbc, closestAbc)) { @@ -235,8 +225,7 @@ bestSpecialization(const TwoQubitWeylDecomposition& decomposition, if (isClose(decomposition.a(), 0., 0.)) { return Specialization::ControlledEquiv; } - if (isClose((std::numbers::pi / 4.0), (std::numbers::pi / 4.0), - decomposition.c())) { + if (isClose(kPiOver4, kPiOver4, decomposition.c())) { return Specialization::MirrorControlledEquiv; } if (isClose((decomposition.a() + decomposition.b()) / 2., @@ -255,58 +244,48 @@ bestSpecialization(const TwoQubitWeylDecomposition& decomposition, return Specialization::General; } -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; - } - auto absDiff = std::abs(lhs - rhs); - if (absDiff <= epsilon) { - return true; - } - auto absLhs = std::abs(lhs); - auto absRhs = std::abs(rhs); - if (absRhs > absLhs) { - return absDiff <= absRhs * maxRelative; - } - return absDiff <= absLhs * maxRelative; +struct ChamberState { + std::array cs{}; + Matrix2x2 k1l; + Matrix2x2 k1r; + Matrix2x2 k2l; + Matrix2x2 k2r; + double globalPhase{}; + double a{}; + double b{}; + double c{}; +}; + +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}; } -TwoQubitWeylDecomposition -TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, - std::optional fidelity) { - auto u = unitaryMatrix; - auto detU = u.determinant(); - // Project into SU(4): det^{-1/4} removes global phase; arg(det)/4 is tracked. - auto detPow = std::pow(detU, -0.25); - u *= detPow; - auto globalPhase = std::arg(detU) / 4.; - - auto uP = magicBasisTransform(u, /*outOfMagicBasis=*/true); +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, kDiagonalizationPrecision); - auto [p, d] = diagonalizeComplexSymmetric(m2, DIAGONALIZATION_PRECISION); - - // Map eigenvalue phases to Weyl coordinates in [0, 2*pi). - constexpr double pi = std::numbers::pi; 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); + cs[i] = remEuclid((dReal[i] + dReal[3]) / 2.0, 2.0 * kPi); } // 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); + const auto tmp = remEuclid(cs[i], kPi / 2.0); + cstemp[i] = std::min(tmp, (kPi / 2.0) - tmp); } std::array order{0, 1, 2}; std::ranges::stable_sort( @@ -334,6 +313,13 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, } 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]); @@ -361,93 +347,93 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, 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)); + 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; - // Map into the Weyl chamber. - 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[0] > (kPi / 2.0)) { + cs[0] -= 3.0 * (kPi / 2.0); + k1l = k1l * iPauliY(); + k1r = k1r * iPauliY(); + globalPhase += (kPi / 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); + if (cs[1] > (kPi / 2.0)) { + cs[1] -= 3.0 * (kPi / 2.0); + k1l = k1l * iPauliX(); + k1r = k1r * iPauliX(); + globalPhase += (kPi / 2.0); } auto conjs = 0; - if (cs[0] > (pi / 4.0)) { - cs[0] = (pi / 2.0) - cs[0]; - K1l = K1l * iPauliY(); - K2r = iPauliY() * K2r; + if (cs[0] > kPiOver4) { + cs[0] = (kPi / 2.0) - cs[0]; + k1l = k1l * iPauliY(); + k2r = iPauliY() * k2r; conjs += 1; - globalPhase -= (pi / 2.0); + globalPhase -= (kPi / 2.0); } - if (cs[1] > (pi / 4.0)) { - cs[1] = (pi / 2.0) - cs[1]; - K1l = K1l * iPauliX(); - K2r = iPauliX() * K2r; + if (cs[1] > kPiOver4) { + cs[1] = (kPi / 2.0) - cs[1]; + k1l = k1l * iPauliX(); + k2r = iPauliX() * k2r; conjs += 1; - globalPhase += (pi / 2.0); + globalPhase += (kPi / 2.0); if (conjs == 1) { - globalPhase -= pi; + globalPhase -= kPi; } } - if (cs[2] > (pi / 2.0)) { - cs[2] -= 3.0 * (pi / 2.0); - K1l = K1l * iPauliZ(); - K1r = K1r * iPauliZ(); - globalPhase += (pi / 2.0); + if (cs[2] > (kPi / 2.0)) { + cs[2] -= 3.0 * (kPi / 2.0); + k1l = k1l * iPauliZ(); + k1r = k1r * iPauliZ(); + globalPhase += (kPi / 2.0); if (conjs == 1) { - globalPhase -= pi; + globalPhase -= kPi; } } if (conjs == 1) { - cs[2] = (pi / 2.0) - cs[2]; - K1l = K1l * iPauliZ(); - K2r = iPauliZ() * K2r; - globalPhase += (pi / 2.0); - } - if (cs[2] > (pi / 4.0)) { - cs[2] -= (pi / 2.0); - K1l = K1l * iPauliZ(); - K1r = K1r * iPauliZ(); - globalPhase -= (pi / 2.0); - } - - // Bind Weyl coordinates to the canonical gate. - auto [a, b, c] = std::tie(cs[1], cs[0], cs[2]); - - TwoQubitWeylDecomposition decomposition; - decomposition.a_ = a; - decomposition.b_ = b; - decomposition.c_ = c; - decomposition.globalPhase_ = globalPhase; - decomposition.k1l_ = K1l; - decomposition.k2l_ = K2l; - decomposition.k1r_ = K1r; - decomposition.k2r_ = K2r; - - assert((Matrix4x4::kron(K1l, K1r) * decomposition.getCanonicalMatrix() * - Matrix4x4::kron(K2l, K2r) * std::exp(Complex{0.0, 1.0} * globalPhase)) - .isApprox(unitaryMatrix, WEYL_TOLERANCE)); + cs[2] = (kPi / 2.0) - cs[2]; + k1l = k1l * iPauliZ(); + k2r = iPauliZ() * k2r; + globalPhase += (kPi / 2.0); + } + if (cs[2] > kPiOver4) { + cs[2] -= (kPi / 2.0); + k1l = k1l * iPauliZ(); + k1r = k1r * iPauliZ(); + globalPhase -= (kPi / 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; +} - auto flippedFromOriginal = decomposition.applySpecialization(fidelity); +} // namespace - auto getTraceValue = [&]() { - if (flippedFromOriginal) { - return getTrace((pi / 2.0) - a, b, -c, decomposition.a_, decomposition.b_, - decomposition.c_); - } - return getTrace(a, b, c, decomposition.a_, decomposition.b_, - decomposition.c_); - }; - const auto trace = getTraceValue(); +//===----------------------------------------------------------------------===// +// TwoQubitWeylDecomposition +//===----------------------------------------------------------------------===// + +void TwoQubitWeylDecomposition::finalizeSpecializationPhase( + bool flippedFromOriginal, double preSpecializationA, + double preSpecializationB, double preSpecializationC, + const std::optional& fidelity) { + const auto trace = + flippedFromOriginal + ? getTrace((kPi / 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( @@ -455,7 +441,34 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, "specialization is worse than requested fidelity ({0:F4} vs {1:F4})!", calculatedFidelity, *fidelity)); } - decomposition.globalPhase_ += std::arg(trace); + globalPhase_ += std::arg(trace); +} + +TwoQubitWeylDecomposition +TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, + std::optional 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); assert((Matrix4x4::kron(decomposition.k1l_, decomposition.k1r_) * decomposition.getCanonicalMatrix() * @@ -468,12 +481,7 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, Matrix4x4 TwoQubitWeylDecomposition::getCanonicalMatrix(double a, double b, double c) { - // U_d(a,b,c) = exp(i(a XX + b YY + c ZZ)); `-2*a/b/c` matches RXX/RYY/RZZ - // gate convention; product order matches Qiskit (RZZ · RYY · RXX). - const auto xx = rxxMatrix(-2.0 * a); - const auto yy = ryyMatrix(-2.0 * b); - const auto zz = rzzMatrix(-2.0 * c); - return zz * yy * xx; + return rzzMatrix(-2.0 * c) * ryyMatrix(-2.0 * b) * rxxMatrix(-2.0 * a); } bool TwoQubitWeylDecomposition::applySpecialization( @@ -484,7 +492,8 @@ bool TwoQubitWeylDecomposition::applySpecialization( return flippedFromOriginal; } - if (newSpecialization == Specialization::IdEquiv) { + switch (newSpecialization) { + case Specialization::IdEquiv: a_ = 0.; b_ = 0.; c_ = 0.; @@ -492,7 +501,8 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2l_ = Matrix2x2::identity(); k1r_ = k1r_ * k2r_; k2r_ = Matrix2x2::identity(); - } else if (newSpecialization == Specialization::SWAPEquiv) { + break; + case Specialization::SWAPEquiv: if (c_ > 0.) { k1l_ = k1l_ * k2r_; k1r_ = k1r_ * k2l_; @@ -500,20 +510,19 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2r_ = Matrix2x2::identity(); } else { flippedFromOriginal = true; - - globalPhase_ += (std::numbers::pi / 2.0); + globalPhase_ += (kPi / 2.0); k1l_ = k1l_ * iPauliZ() * k2r_; k1r_ = k1r_ * iPauliZ() * k2l_; k2l_ = Matrix2x2::identity(); k2r_ = Matrix2x2::identity(); } - a_ = (std::numbers::pi / 4.0); - b_ = (std::numbers::pi / 4.0); - c_ = (std::numbers::pi / 4.0); - } else if (newSpecialization == Specialization::PartialSWAPEquiv) { - auto closest = closestPartialSwap(a_, b_, c_); - auto k2lDagger = k2l_.adjoint(); - + a_ = kPiOver4; + b_ = kPiOver4; + c_ = kPiOver4; + break; + case Specialization::PartialSWAPEquiv: { + const auto closest = closestPartialSwap(a_, b_, c_); + const auto k2lDagger = k2l_.adjoint(); a_ = closest; b_ = closest; c_ = closest; @@ -521,10 +530,11 @@ bool TwoQubitWeylDecomposition::applySpecialization( k1r_ = k1r_ * k2l_; k2r_ = k2lDagger * k2r_; k2l_ = Matrix2x2::identity(); - } else if (newSpecialization == Specialization::PartialSWAPFlipEquiv) { - auto closest = closestPartialSwap(a_, b_, -c_); - auto k2lDagger = k2l_.adjoint(); - + break; + } + case Specialization::PartialSWAPFlipEquiv: { + const auto closest = closestPartialSwap(a_, b_, -c_); + const auto k2lDagger = k2l_.adjoint(); a_ = closest; b_ = closest; c_ = -closest; @@ -532,10 +542,11 @@ bool TwoQubitWeylDecomposition::applySpecialization( k1r_ = k1r_ * iPauliZ() * k2l_ * iPauliZ(); k2r_ = iPauliZ() * k2lDagger * iPauliZ() * k2r_; k2l_ = Matrix2x2::identity(); - } else if (newSpecialization == Specialization::ControlledEquiv) { - const EulerBasis eulerBasis = EulerBasis::XYX; + break; + } + case Specialization::ControlledEquiv: { const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = - anglesFromUnitary(k2l_, eulerBasis); + anglesFromUnitary(k2l_, EulerBasis::XYX); const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = anglesFromUnitary(k2r_, EulerBasis::XYX); b_ = 0.; @@ -545,23 +556,26 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); k1r_ = k1r_ * rxMatrix(k2rphi); k2r_ = ryMatrix(k2rtheta) * rxMatrix(k2rlambda); - } else if (newSpecialization == Specialization::MirrorControlledEquiv) { + 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_ = (std::numbers::pi / 4.0); - b_ = (std::numbers::pi / 4.0); + a_ = kPiOver4; + b_ = kPiOver4; globalPhase_ = globalPhase_ + k2lphase + k2rphase; k1l_ = k1l_ * rzMatrix(k2rphi); k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); k1r_ = k1r_ * rzMatrix(k2lphi); k2r_ = ryMatrix(k2rtheta) * rzMatrix(k2rlambda); - } else if (newSpecialization == Specialization::FSimaabEquiv) { - auto [k2ltheta, k2lphi, k2llambda, k2lphase] = + break; + } + case Specialization::FSimaabEquiv: { + const auto [k2ltheta, k2lphi, k2llambda, k2lphase] = anglesFromUnitary(k2l_, EulerBasis::ZYZ); - auto ab = (a_ + b_) / 2.; - + const auto ab = (a_ + b_) / 2.; a_ = ab; b_ = ab; globalPhase_ = globalPhase_ + k2lphase; @@ -569,12 +583,12 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); k1r_ = k1r_ * rzMatrix(k2lphi); k2r_ = rzMatrix(-k2lphi) * k2r_; - } else if (newSpecialization == Specialization::FSimabbEquiv) { - auto eulerBasis = EulerBasis::XYX; - auto [k2ltheta, k2lphi, k2llambda, k2lphase] = - anglesFromUnitary(k2l_, eulerBasis); - auto bc = (b_ + c_) / 2.; - + 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; @@ -582,12 +596,12 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); k1r_ = k1r_ * rxMatrix(k2lphi); k2r_ = rxMatrix(-k2lphi) * k2r_; - } else if (newSpecialization == Specialization::FSimabmbEquiv) { - auto eulerBasis = EulerBasis::XYX; - auto [k2ltheta, k2lphi, k2llambda, k2lphase] = - anglesFromUnitary(k2l_, eulerBasis); - auto bc = (b_ - c_) / 2.; - + 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; @@ -595,263 +609,12 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2l_ = ryMatrix(k2ltheta) * rxMatrix(k2llambda); k1r_ = k1r_ * iPauliZ() * rxMatrix(k2lphi) * iPauliZ(); k2r_ = iPauliZ() * rxMatrix(-k2lphi) * iPauliZ() * k2r_; - } else { - llvm::reportFatalInternalError( - "Unknown specialization for Weyl decomposition!"); + break; } - return flippedFromOriginal; -} - -TwoQubitBasisDecomposer -TwoQubitBasisDecomposer::create(const Matrix4x4& basisMatrix, - double basisFidelity) { - const Matrix2x2 k12RArr = Matrix2x2::fromElements( - 1i * INV_SQRT2, INV_SQRT2, -INV_SQRT2, -1i * INV_SQRT2); - const Matrix2x2 k12LArr = - Matrix2x2::fromElements(Complex{0.5, 0.5}, Complex{0.5, 0.5}, - Complex{-0.5, 0.5}, Complex{0.5, -0.5}); - - const auto basisDecomposer = - TwoQubitWeylDecomposition::create(basisMatrix, DEFAULT_FIDELITY); - const auto isSuperControlled = - relativeEq(basisDecomposer.a(), std::numbers::pi / 4.0, 1e-13, 1e-09) && - relativeEq(basisDecomposer.c(), 0.0, 1e-13, 1e-09); - - auto b = basisDecomposer.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 k22l = - Matrix2x2::fromElements(INV_SQRT2, -INV_SQRT2, INV_SQRT2, INV_SQRT2); - const Matrix2x2 k22r = Matrix2x2::fromElements(0, 1, -1, 0); - 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))); - auto k1lDagger = basisDecomposer.k1l().adjoint(); - auto k1rDagger = basisDecomposer.k1r().adjoint(); - auto k2lDagger = basisDecomposer.k2l().adjoint(); - auto k2rDagger = basisDecomposer.k2r().adjoint(); - auto u0l = k31l * k1lDagger; - auto u0r = k31r * k1rDagger; - auto u1l = k2lDagger * k32lK21l * k1lDagger; - auto u1ra = k2rDagger * k32r; - auto u1rb = k21r * k1rDagger; - auto u2la = k2lDagger * k22l; - auto u2lb = k11l * k1lDagger; - auto u2ra = k2rDagger * k22r; - auto u2rb = k11r * k1rDagger; - auto u3l = k2lDagger * k12LArr; - auto u3r = k2rDagger * k12RArr; - auto q0l = k12LArr.adjoint() * k1lDagger; - auto q0r = k12RArr.adjoint() * iPauliZ() * k1rDagger; - auto q1la = k2lDagger * k11l.adjoint(); - auto q1lb = k11l * k1lDagger; - auto q1ra = k2rDagger * iPauliZ() * k11r.adjoint(); - auto q1rb = k11r * k1rDagger; - auto q2l = k2lDagger * k12LArr; - auto q2r = k2rDagger * k12RArr; - - TwoQubitBasisDecomposer decomposer; - decomposer.basisFidelity = basisFidelity; - decomposer.basisDecomposer = basisDecomposer; - decomposer.isSuperControlled = isSuperControlled; - decomposer.smb.u0l = u0l; - decomposer.smb.u0r = u0r; - decomposer.smb.u1l = u1l; - decomposer.smb.u1ra = u1ra; - decomposer.smb.u1rb = u1rb; - decomposer.smb.u2la = u2la; - decomposer.smb.u2lb = u2lb; - decomposer.smb.u2ra = u2ra; - decomposer.smb.u2rb = u2rb; - decomposer.smb.u3l = u3l; - decomposer.smb.u3r = u3r; - decomposer.smb.q0l = q0l; - decomposer.smb.q0r = q0r; - decomposer.smb.q1la = q1la; - decomposer.smb.q1lb = q1lb; - decomposer.smb.q1ra = q1ra; - decomposer.smb.q1rb = q1rb; - decomposer.smb.q2l = q2l; - decomposer.smb.q2r = q2r; - return decomposer; -} - -std::optional -TwoQubitBasisDecomposer::decomposeTarget( - const Matrix4x4& targetUnitary, - const std::optional numBasisGateUses) const { - const auto targetWeyl = - TwoQubitWeylDecomposition::create(targetUnitary, DEFAULT_FIDELITY); - return twoQubitDecompose(targetWeyl, numBasisGateUses); -} - -std::optional -TwoQubitBasisDecomposer::twoQubitDecompose( - const TwoQubitWeylDecomposition& targetDecomposition, - std::optional numBasisGateUses) const { - auto traces = this->traces(targetDecomposition); - auto getDefaultNbasis = [&]() -> std::uint8_t { - // Maximize traceToFidelity(traces[i]) * basisFidelity^i over i in {0..3}. - auto bestValue = std::numeric_limits::lowest(); - auto bestIndex = -1; - for (int i = 0; std::cmp_less(i, traces.size()); ++i) { - auto value = traceToFidelity(traces[i]) * std::pow(basisFidelity, i); - 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); - }; - auto bestNbasis = numBasisGateUses.value_or(getDefaultNbasis()); - if (bestNbasis > 1 && !isSuperControlled) { - // More than one basis gate requires a super-controlled basis (e.g. CX). - return std::nullopt; - } - auto chooseDecomposition = [&]() { - if (bestNbasis == 0) { - return decomp0(targetDecomposition); - } - if (bestNbasis == 1) { - return decomp1(targetDecomposition); - } - if (bestNbasis == 2) { - return decomp2Supercontrolled(targetDecomposition); - } - if (bestNbasis == 3) { - return decomp3Supercontrolled(targetDecomposition); - } - llvm::reportFatalInternalError( - "Invalid number of basis gates to use in basis decomposition (" + - llvm::Twine(bestNbasis) + ")!"); - llvm_unreachable(""); - }; - SmallVector factors = chooseDecomposition(); - - double globalPhase = targetDecomposition.globalPhase(); - globalPhase -= bestNbasis * basisDecomposer.globalPhase(); - if (bestNbasis == 2) { - // Two-basis template is exact up to a pi global-phase offset. - globalPhase += std::numbers::pi; + case Specialization::General: + llvm_unreachable("unreachable specialization"); } - globalPhase = remEuclid(globalPhase, 2.0 * std::numbers::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 { - // One basis gate; reliable only when the target is in the Weyl chamber. - return SmallVector{ - basisDecomposer.k2r().adjoint() * target.k2r(), - basisDecomposer.k2l().adjoint() * target.k2l(), - target.k1r() * basisDecomposer.k1r().adjoint(), - target.k1l() * basisDecomposer.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 { - // Hilbert-Schmidt traces for 0..3 basis uses; index i == numBasisUses. - // i=0: identity; i=1: one basis gate; i=2: two uses; i=3: exact (trace=4). - 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((std::numbers::pi / 4.0) - target.a()) * - std::cos(basisDecomposer.b() - target.b()) * - std::cos(target.c()), - std::sin((std::numbers::pi / 4.0) - target.a()) * - std::sin(basisDecomposer.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); + return flippedFromOriginal; } } // namespace mlir::qco::decomposition From 341cfc88cc3b612451c3677dccd5ca9985653ac7 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 15:58:16 +0200 Subject: [PATCH 14/16] =?UTF-8?q?=F0=9F=9A=A8=20Fix=20linter=20warnings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Weyl.cpp | 122 +++++++++--------- 1 file changed, 62 insertions(+), 60 deletions(-) diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index 0ba70c6b9a..32a0ba391f 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -49,9 +49,23 @@ enum class Specialization : std::uint8_t { FSimabmbEquiv, }; -static constexpr double kDiagonalizationPrecision = 1e-13; -static constexpr double kPi = std::numbers::pi; -static constexpr double kPiOver4 = kPi / 4.0; +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) { @@ -67,12 +81,12 @@ static double traceToFidelity(const Complex& trace) { return (4.0 + (traceAbs * traceAbs)) / 20.0; } -static const Matrix4x4 kMagicBasisNonNormalized = Matrix4x4::fromElements( // - 1, 1i, 0, 0, // - 0, 0, 1i, 1, // - 0, 0, 1i, -1, // +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 kMagicBasisNonNormalizedDagger = +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}, // @@ -82,9 +96,11 @@ static const Matrix4x4 kMagicBasisNonNormalizedDagger = static Matrix4x4 magicBasisTransform(const Matrix4x4& unitary, bool outOfMagicBasis) { if (outOfMagicBasis) { - return kMagicBasisNonNormalizedDagger * unitary * kMagicBasisNonNormalized; + return MAGIC_BASIS_NON_NORMALIZED_DAGGER * unitary * + MAGIC_BASIS_NON_NORMALIZED; } - return kMagicBasisNonNormalized * unitary * kMagicBasisNonNormalizedDagger; + return MAGIC_BASIS_NON_NORMALIZED * unitary * + MAGIC_BASIS_NON_NORMALIZED_DAGGER; } static double closestPartialSwap(double a, double b, double c) { @@ -212,8 +228,8 @@ bestSpecialization(const TwoQubitWeylDecomposition& decomposition, if (isClose(0., 0., 0.)) { return Specialization::IdEquiv; } - if (isClose(kPiOver4, kPiOver4, kPiOver4) || - isClose(kPiOver4, kPiOver4, -kPiOver4)) { + 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)) { @@ -225,7 +241,7 @@ bestSpecialization(const TwoQubitWeylDecomposition& decomposition, if (isClose(decomposition.a(), 0., 0.)) { return Specialization::ControlledEquiv; } - if (isClose(kPiOver4, kPiOver4, decomposition.c())) { + if (isClose(PI_OVER_4, PI_OVER_4, decomposition.c())) { return Specialization::MirrorControlledEquiv; } if (isClose((decomposition.a() + decomposition.b()) / 2., @@ -244,18 +260,6 @@ bestSpecialization(const TwoQubitWeylDecomposition& decomposition, return Specialization::General; } -struct ChamberState { - std::array cs{}; - Matrix2x2 k1l; - Matrix2x2 k1r; - Matrix2x2 k2l; - Matrix2x2 k2r; - double globalPhase{}; - double a{}; - double b{}; - double c{}; -}; - static std::pair projectToSU4(const Matrix4x4& unitary) { auto u = unitary; const auto detU = u.determinant(); @@ -268,7 +272,7 @@ static std::tuple, computeOrderedWeylCoordinates(const Matrix4x4& u) { const auto uP = magicBasisTransform(u, /*outOfMagicBasis=*/true); const Matrix4x4 m2 = uP.transpose() * uP; - auto [p, d] = diagonalizeComplexSymmetric(m2, kDiagonalizationPrecision); + auto [p, d] = diagonalizeComplexSymmetric(m2, DIAGONALIZATION_PRECISION); std::array dReal{}; for (std::size_t i = 0; i < d.size(); ++i) { @@ -278,14 +282,14 @@ computeOrderedWeylCoordinates(const Matrix4x4& u) { std::array cs{}; for (std::size_t i = 0; i < cs.size(); ++i) { - cs[i] = remEuclid((dReal[i] + dReal[3]) / 2.0, 2.0 * kPi); + 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], kPi / 2.0); - cstemp[i] = std::min(tmp, (kPi / 2.0) - tmp); + 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( @@ -353,56 +357,56 @@ static ChamberState buildChamberState(const Matrix4x4& u, const Matrix4x4& uP, assert(Matrix4x4::kron(k2l, k2r).isApprox(k2, WEYL_TOLERANCE)); globalPhase += phaseL + phaseR; - if (cs[0] > (kPi / 2.0)) { - cs[0] -= 3.0 * (kPi / 2.0); + if (cs[0] > (PI / 2.0)) { + cs[0] -= 3.0 * (PI / 2.0); k1l = k1l * iPauliY(); k1r = k1r * iPauliY(); - globalPhase += (kPi / 2.0); + globalPhase += (PI / 2.0); } - if (cs[1] > (kPi / 2.0)) { - cs[1] -= 3.0 * (kPi / 2.0); + if (cs[1] > (PI / 2.0)) { + cs[1] -= 3.0 * (PI / 2.0); k1l = k1l * iPauliX(); k1r = k1r * iPauliX(); - globalPhase += (kPi / 2.0); + globalPhase += (PI / 2.0); } auto conjs = 0; - if (cs[0] > kPiOver4) { - cs[0] = (kPi / 2.0) - cs[0]; + if (cs[0] > PI_OVER_4) { + cs[0] = (PI / 2.0) - cs[0]; k1l = k1l * iPauliY(); k2r = iPauliY() * k2r; conjs += 1; - globalPhase -= (kPi / 2.0); + globalPhase -= (PI / 2.0); } - if (cs[1] > kPiOver4) { - cs[1] = (kPi / 2.0) - cs[1]; + if (cs[1] > PI_OVER_4) { + cs[1] = (PI / 2.0) - cs[1]; k1l = k1l * iPauliX(); k2r = iPauliX() * k2r; conjs += 1; - globalPhase += (kPi / 2.0); + globalPhase += (PI / 2.0); if (conjs == 1) { - globalPhase -= kPi; + globalPhase -= PI; } } - if (cs[2] > (kPi / 2.0)) { - cs[2] -= 3.0 * (kPi / 2.0); + if (cs[2] > (PI / 2.0)) { + cs[2] -= 3.0 * (PI / 2.0); k1l = k1l * iPauliZ(); k1r = k1r * iPauliZ(); - globalPhase += (kPi / 2.0); + globalPhase += (PI / 2.0); if (conjs == 1) { - globalPhase -= kPi; + globalPhase -= PI; } } if (conjs == 1) { - cs[2] = (kPi / 2.0) - cs[2]; + cs[2] = (PI / 2.0) - cs[2]; k1l = k1l * iPauliZ(); k2r = iPauliZ() * k2r; - globalPhase += (kPi / 2.0); + globalPhase += (PI / 2.0); } - if (cs[2] > kPiOver4) { - cs[2] -= (kPi / 2.0); + if (cs[2] > PI_OVER_4) { + cs[2] -= (PI / 2.0); k1l = k1l * iPauliZ(); k1r = k1r * iPauliZ(); - globalPhase -= (kPi / 2.0); + globalPhase -= (PI / 2.0); } ChamberState chamber; @@ -418,8 +422,6 @@ static ChamberState buildChamberState(const Matrix4x4& u, const Matrix4x4& uP, return chamber; } -} // namespace - //===----------------------------------------------------------------------===// // TwoQubitWeylDecomposition //===----------------------------------------------------------------------===// @@ -430,7 +432,7 @@ void TwoQubitWeylDecomposition::finalizeSpecializationPhase( const std::optional& fidelity) { const auto trace = flippedFromOriginal - ? getTrace((kPi / 2.0) - preSpecializationA, preSpecializationB, + ? getTrace((PI / 2.0) - preSpecializationA, preSpecializationB, -preSpecializationC, a_, b_, c_) : getTrace(preSpecializationA, preSpecializationB, preSpecializationC, a_, b_, c_); @@ -510,15 +512,15 @@ bool TwoQubitWeylDecomposition::applySpecialization( k2r_ = Matrix2x2::identity(); } else { flippedFromOriginal = true; - globalPhase_ += (kPi / 2.0); + globalPhase_ += (PI / 2.0); k1l_ = k1l_ * iPauliZ() * k2r_; k1r_ = k1r_ * iPauliZ() * k2l_; k2l_ = Matrix2x2::identity(); k2r_ = Matrix2x2::identity(); } - a_ = kPiOver4; - b_ = kPiOver4; - c_ = kPiOver4; + a_ = PI_OVER_4; + b_ = PI_OVER_4; + c_ = PI_OVER_4; break; case Specialization::PartialSWAPEquiv: { const auto closest = closestPartialSwap(a_, b_, c_); @@ -563,8 +565,8 @@ bool TwoQubitWeylDecomposition::applySpecialization( anglesFromUnitary(k2l_, EulerBasis::ZYZ); const auto [k2rtheta, k2rphi, k2rlambda, k2rphase] = anglesFromUnitary(k2r_, EulerBasis::ZYZ); - a_ = kPiOver4; - b_ = kPiOver4; + a_ = PI_OVER_4; + b_ = PI_OVER_4; globalPhase_ = globalPhase_ + k2lphase + k2rphase; k1l_ = k1l_ * rzMatrix(k2rphi); k2l_ = ryMatrix(k2ltheta) * rzMatrix(k2llambda); From 8114583941829d0548273ea126b5e82d92f61024 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 16:34:43 +0200 Subject: [PATCH 15/16] =?UTF-8?q?=F0=9F=90=87=20Address=20rabbit's=20comme?= =?UTF-8?q?nts?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../QCO/Transforms/Decomposition/Weyl.h | 15 ++++++++++++ .../Decomposition/BasisDecomposer.cpp | 13 +++++++---- .../QCO/Transforms/Decomposition/Weyl.cpp | 23 +++++++++++++++---- .../Decomposition/test_weyl_decomposition.cpp | 8 +++++++ 4 files changed, 50 insertions(+), 9 deletions(-) diff --git a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h index 5a0483628f..d5db42398d 100644 --- a/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h +++ b/mlir/include/mlir/Dialect/QCO/Transforms/Decomposition/Weyl.h @@ -129,6 +129,19 @@ class TwoQubitBasisDecomposer { [[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; @@ -177,6 +190,8 @@ class TwoQubitBasisDecomposer { [[nodiscard]] std::array, 4> traces(const TwoQubitWeylDecomposition& target) const; + TwoQubitBasisDecomposer() = default; + double basisFidelity{}; TwoQubitWeylDecomposition basisWeyl; bool isSuperControlled{}; diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp index 30c5f5709b..814702544d 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/BasisDecomposer.cpp @@ -82,6 +82,14 @@ static double traceToFidelity(const Complex& trace) { 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 = @@ -202,10 +210,7 @@ TwoQubitBasisDecomposer::twoQubitDecompose( factors = decomp3Supercontrolled(targetDecomposition); break; default: - llvm::reportFatalInternalError(llvm::formatv( - "Invalid number of basis gates to use in basis decomposition ({0})!", - bestNbasis)); - llvm_unreachable(""); + return std::nullopt; } double globalPhase = targetDecomposition.globalPhase(); diff --git a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp index 32a0ba391f..a857fdc9cc 100644 --- a/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp +++ b/mlir/lib/Dialect/QCO/Transforms/Decomposition/Weyl.cpp @@ -449,6 +449,14 @@ void TwoQubitWeylDecomposition::finalizeSpecializationPhase( 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); @@ -472,11 +480,16 @@ TwoQubitWeylDecomposition::create(const Matrix4x4& unitaryMatrix, decomposition.finalizeSpecializationPhase(flippedFromOriginal, chamber.a, chamber.b, chamber.c, fidelity); - 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 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; } diff --git a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp index bcdc627c60..bb11c6fb2b 100644 --- a/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp +++ b/mlir/unittests/Dialect/QCO/Transforms/Decomposition/test_weyl_decomposition.cpp @@ -339,6 +339,14 @@ TEST(BasisDecomposerTest, RejectsMultipleBasisUsesForNonSuperControlledBasis) { 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); From 004ee4041879296e88a84a5385a3d2360d1bbadd Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Tue, 23 Jun 2026 16:49:59 +0200 Subject: [PATCH 16/16] =?UTF-8?q?=F0=9F=93=9C=20Update=20CHANGELOG.md=20to?= =?UTF-8?q?=20include=20PR=20#1803?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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