[bugfix] Handle missing edge case in decomposition of Rotation (#1208)

pytket/conanfile.py

@@ -32,7 +32,7 @@ def package(self):
         cmake.install()
 
     def requirements(self):
-        self.requires("tket/1.2.83@tket/stable")
+        self.requires("tket/1.2.84@tket/stable")
         self.requires("tklog/0.3.3@tket/stable")
         self.requires("tkrng/0.3.3@tket/stable")
         self.requires("tkassert/0.3.4@tket/stable")

pytket/docs/changelog.rst

@@ -1,6 +1,13 @@
 Changelog
 =========
 
+Unreleased
+----------
+
+Fixes:
+
+* Handle a missing edge case in decomposition of single-qubit rotations.
+
 1.23.0 (January 2024)
 ---------------------
 

tket/conanfile.py

@@ -23,7 +23,7 @@
 
 class TketConan(ConanFile):
     name = "tket"
-    version = "1.2.83"
+    version = "1.2.84"
     package_type = "library"
     license = "Apache 2"
     homepage = "https://github.com/CQCL/tket"

tket/src/Gate/Rotation.cpp

@@ -14,7 +14,10 @@
 
 #include "tket/Gate/Rotation.hpp"
 
-#include "symengine/constants.h"
+#include <symengine/constants.h>
+
+#include <tkassert/Assert.hpp>
+
 #include "tket/OpType/OpDesc.hpp"
 #include "tket/OpType/OpType.hpp"
 #include "tket/Utils/Expression.hpp"
@@ -103,6 +106,22 @@ static std::tuple<Expr, Expr, Expr> xyx_angles_from_coeffs(
   if (i_zero && j_zero) return {-0.5, 2 * atan2_bypi(k, s), 0.5};
   if (i_zero && k_zero) return {0, 2 * atan2_bypi(j, s), 0};
   if (j_zero && k_zero) return {2 * atan2_bypi(i, s), 0, 0};
+  if (s_zero) {
+    // We want to return (a, b, c) s.t.
+    // Rx(c) Ry(b) Rx(a) = i*i_ + j*j_ + k*k_ (where i_, j_, k_ are the
+    // quaternionic roots of -1. Multiplying out the LHS gives
+    // (Cc + i_ Sc)(Cb + j_ Sb)(Ca + i_ Sa)
+    // where for brevity Sa = sin(pi*a/2), Ca = cos(pi*a/2) etc.
+    // Expanding and using the sin/cos addition identities gives
+    // Cb cos(pi*l) = 0, Cb sin(pi*l) = i, Sb cos(pi*m) = j, Sb sin(pi*m) = k
+    // where l = (c+a)/2 and m = (c-a)/2.
+    // Since we have already handled the case s=i=0 we know that Cb != 0, so
+    // cos(pi*l) = 0. A solution to these equations is given by
+    // l = 1/2, m = atan2(k,j)/pi, b = (2/pi)*acos(i)
+    // and finally substituting a=l-m, b=l+m gives the following solution:
+    Expr m = atan2_bypi(k, j);
+    return {0.5 - m, 2 * acos_bypi(i), 0.5 + m};
+  }
 
   // This is a (partial) workaround for
   // https://github.com/symengine/symengine/issues/1806
@@ -124,6 +143,7 @@ static std::tuple<Expr, Expr, Expr> xyx_angles_from_coeffs(
   // 2 * atan2(B, A).
   // Finally, note that u must be well-defined because we have already dealt
   // with all cases where s = 0.
+  TKET_ASSERT(!s_zero);
   if (approx_0(SymEngine::expand(i * j + s * k))) {
     Expr u = expr_div(i, s);
     if (SymEngine::free_symbols(u).empty()) {

tket/test/CMakeLists.txt

@@ -105,6 +105,7 @@ add_executable(test-tket
     src/OpType/test_OpTypeFunctions.cpp
     src/Passes/test_SynthesiseTK.cpp
     src/Passes/test_SynthesiseTket.cpp
+    src/Gate/test_Rotation.cpp
     src/Gate/test_GateUnitaryMatrix.cpp
     src/Simulation/test_CircuitSimulator.cpp
     src/Simulation/test_PauliExpBoxUnitaryCalculator.cpp

tket/test/src/Gate/test_Rotation.cpp

@@ -0,0 +1,58 @@
+// Copyright 2019-2023 Cambridge Quantum Computing
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include <catch2/catch_test_macros.hpp>
+#include <vector>
+
+#include "../testutil.hpp"
+#include "tket/Circuit/Circuit.hpp"
+#include "tket/Gate/Rotation.hpp"
+#include "tket/OpType/OpType.hpp"
+#include "tket/Utils/Expression.hpp"
+
+namespace tket {
+
+SCENARIO("xyx decomposition") {
+  const std::vector<std::tuple<Expr, Expr, Expr>> angles = {
+      {0.2, 0.3, 0.4},
+      {0.4, 0.8, 1.4},
+      {0.4, 0.8, 0.6},
+  };
+  for (const std::tuple<Expr, Expr, Expr> &abc : angles) {
+    Expr a = std::get<0>(abc);
+    Expr b = std::get<1>(abc);
+    Expr c = std::get<2>(abc);
+    Rotation rxa(OpType::Rx, a);
+    Rotation ryb(OpType::Ry, b);
+    Rotation rxc(OpType::Rx, c);
+    Rotation r(rxa);
+    r.apply(ryb);
+    r.apply(rxc);
+    std::tuple<Expr, Expr, Expr> abc1 = r.to_pqp(OpType::Rx, OpType::Ry);
+    Expr a1 = std::get<0>(abc1);
+    Expr b1 = std::get<1>(abc1);
+    Expr c1 = std::get<2>(abc1);
+    Circuit circ(1);
+    circ.add_op<unsigned>(OpType::Rx, a, {0});
+    circ.add_op<unsigned>(OpType::Ry, b, {0});
+    circ.add_op<unsigned>(OpType::Rx, c, {0});
+    Circuit circ1(1);
+    circ1.add_op<unsigned>(OpType::Rx, a1, {0});
+    circ1.add_op<unsigned>(OpType::Ry, b1, {0});
+    circ1.add_op<unsigned>(OpType::Rx, c1, {0});
+    REQUIRE(test_unitary_comparison(circ, circ1));
+  }
+}
+
+}  // namespace tket