Add special case for single-CNOT decomposition of controlled-U

opensquirrel/cnot_decomposer.py

@@ -1,7 +1,8 @@
 import math
 
+from opensquirrel import merger
 from opensquirrel.common import ATOL
-from opensquirrel.default_gates import cnot, ry, rz
+from opensquirrel.default_gates import cnot, ry, rz, x
 from opensquirrel.identity_filter import filter_out_identities
 from opensquirrel.replacer import Decomposer
 from opensquirrel.squirrel_ir import BlochSphereRotation, ControlledGate, Float, Gate
@@ -29,30 +30,36 @@ def decompose(g: Gate) -> [Gate]:
             # ControlledGate's with 2+ control qubits are ignored.
             return [g]
 
+        target_qubit = g.target_gate.qubit
+
         # Perform ZYZ decomposition on the target gate.
         # This gives us an ABC decomposition (U = AXBXC, ABC = I) of the target gate.
         # See https://threeplusone.com/pubs/on_gates.pdf
-        theta0, theta1, theta2 = get_zyz_decomposition_angles(g.target_gate.angle, g.target_gate.axis)
-        target_qubit = g.target_gate.qubit
 
-        # First try to see if we can get away with a single CNOT.
-        # FIXME: see https://github.com/QuTech-Delft/OpenSquirrel/issues/99 this could be extended, I believe.
-        if abs(abs(theta0 + theta2) % (2 * math.pi)) < ATOL and abs(abs(theta1 - math.pi) % (2 * math.pi)) < ATOL:
-            # g == rz(theta0) Y rz(theta2) == rz(theta0 - pi / 2) X rz(theta2 + pi / 2)
-            # theta0 + theta2 == 0
+        # Try special case first, see https://arxiv.org/pdf/quant-ph/9503016.pdf lemma 5.5
+        controlled_rotation_times_x = merger.compose_bloch_sphere_rotations(x(target_qubit), g.target_gate)
+        theta0_with_x, theta1_with_x, theta2_with_x = get_zyz_decomposition_angles(
+            controlled_rotation_times_x.angle, controlled_rotation_times_x.axis
+        )
+        if abs((theta0_with_x - theta2_with_x) % (2 * math.pi)) < ATOL:
+            # The decomposition can use a single CNOT according to the lemma.
+
+            A = [ry(q=target_qubit, theta=Float(-theta1_with_x / 2)), rz(q=target_qubit, theta=Float(-theta2_with_x))]
 
-            alpha0 = theta0 - math.pi / 2
-            alpha2 = theta2 + math.pi / 2
+            B = [
+                rz(q=target_qubit, theta=Float(theta2_with_x)),
+                ry(q=target_qubit, theta=Float(theta1_with_x / 2)),
+            ]
 
             return filter_out_identities(
-                [
-                    rz(q=target_qubit, theta=Float(alpha2)),
-                    cnot(control=g.control_qubit, target=target_qubit),
-                    rz(q=target_qubit, theta=Float(alpha0)),
-                    rz(q=g.control_qubit, theta=Float(g.target_gate.phase - math.pi / 2)),
-                ]
+                B
+                + [cnot(control=g.control_qubit, target=target_qubit)]
+                + A
+                + [rz(q=g.control_qubit, theta=Float(g.target_gate.phase - math.pi / 2))]
             )
 
+        theta0, theta1, theta2 = get_zyz_decomposition_angles(g.target_gate.angle, g.target_gate.axis)
+
         A = [ry(q=target_qubit, theta=Float(theta1 / 2)), rz(q=target_qubit, theta=Float(theta2))]
 
         B = [

opensquirrel/merger.py

@@ -16,9 +16,11 @@ def compose_bloch_sphere_rotations(a: BlochSphereRotation, b: BlochSphereRotatio
     """
     assert a.qubit == b.qubit, "Cannot merge two BlochSphereRotation's on different qubits"
 
-    combined_angle = 2 * acos(
-        cos(a.angle / 2) * cos(b.angle / 2) - sin(a.angle / 2) * sin(b.angle / 2) * np.dot(a.axis, b.axis)
-    )
+    acos_argument = cos(a.angle / 2) * cos(b.angle / 2) - sin(a.angle / 2) * sin(b.angle / 2) * np.dot(a.axis, b.axis)
+    # This fixes float approximations like 1.0000000000002 which acos doesn't like.
+    acos_argument = max(min(acos_argument, 1.0), -1.0)
+
+    combined_angle = 2 * acos(acos_argument)
 
     if abs(sin(combined_angle / 2)) < ATOL:
         return BlochSphereRotation.identity(a.qubit)

test/test_cnot_decomposer.py

@@ -29,13 +29,12 @@ def test_cnot(self):
     def test_cz(self):
         self.assertEqual(
             CNOTDecomposer.decompose(cz(Qubit(0), Qubit(1))),
-            # FIXME: this should only be H-CNOT-H no? Check https://github.com/QuTech-Delft/OpenSquirrel/issues/99
             [
-                rz(Qubit(1), Float(math.pi / 2)),
+                rz(Qubit(1), Float(math.pi)),
+                ry(Qubit(1), Float(math.pi / 2)),
                 cnot(Qubit(0), Qubit(1)),
-                rz(Qubit(1), Float(-math.pi / 2)),
-                cnot(Qubit(0), Qubit(1)),
-                rz(Qubit(0), Float(math.pi / 2)),
+                ry(Qubit(1), Float(-math.pi / 2)),
+                rz(Qubit(1), Float(math.pi)),
             ],
         )
 

test/test_integration.py

@@ -166,9 +166,6 @@ def test_libqasm_error(self):
                 use_libqasm=True,
             )
 
-    @unittest.skipIf(
-        importlib.util.find_spec("quantify_scheduler") is None, reason="quantify_scheduler is not installed"
-    )
     def test_export_quantify_scheduler(self):
         myCircuit = Circuit.from_string(
             """
@@ -177,12 +174,16 @@ def test_export_quantify_scheduler(self):
                 qubit[3] qreg
 
                 h qreg[1]
+                cz qreg[0], qreg[1]
+                cnot qreg[0], qreg[1]
                 crk qreg[0], qreg[1], 4
                 h qreg[0]
             """
         )
 
         myCircuit.decompose(decomposer=CNOTDecomposer)
+
+        # Quantify-scheduler prefers CZ.
         myCircuit.replace(
             cnot,
             lambda control, target: [
@@ -191,34 +192,51 @@ def test_export_quantify_scheduler(self):
                 h(target),
             ],
         )
-        myCircuit.merge_single_qubit_gates()
-        myCircuit.decompose(decomposer=ZYZDecomposer)  # FIXME: for best gate count we need a Z-XY decomposer.
-
-        exported_schedule = myCircuit.export(format=ExportFormat.QUANTIFY_SCHEDULER)
-
-        self.assertEqual(exported_schedule.name, "Exported OpenSquirrel circuit")
 
-        operation_ids = [v["operation_id"] for k, v in exported_schedule.schedulables.items()]
-        operations = [exported_schedule.operations[operation_id].name for operation_id in operation_ids]
+        # Diminish gate count by single-qubit gate fusion.
+        myCircuit.merge_single_qubit_gates()
 
-        self.assertEqual(
-            operations,
-            [
-                "Rz(1.5707963, 'qreg[1]')",
-                "Rxy(0.19634954, 1.5707963, 'qreg[1]')",
-                "Rz(-1.5707963, 'qreg[1]')",
-                "CZ (qreg[0], qreg[1])",
-                "Rz(1.5707963, 'qreg[1]')",
-                "Rxy(-0.19634954, 1.5707963, 'qreg[1]')",
-                "Rz(-1.5707963, 'qreg[1]')",
-                "CZ (qreg[0], qreg[1])",
-                "Rz(0.19634954, 'qreg[0]')",
-                "Rxy(-1.5707963, 1.5707963, 'qreg[0]')",
-                "Rz(3.1415927, 'qreg[0]')",
-                "Rz(3.1415927, 'qreg[1]')",
-                "Rxy(1.5707963, 1.5707963, 'qreg[1]')",
-            ],
-        )
+        # FIXME: for best gate count we need a Z-XY decomposer.
+        # See https://github.com/QuTech-Delft/OpenSquirrel/issues/98
+        myCircuit.decompose(decomposer=ZYZDecomposer)
+
+        if importlib.util.find_spec("quantify_scheduler") is None:
+            with self.assertRaisesRegex(
+                Exception, "quantify-scheduler is not installed, or cannot be installed on " "your system"
+            ):
+                myCircuit.export(format=ExportFormat.QUANTIFY_SCHEDULER)
+        else:
+            exported_schedule = myCircuit.export(format=ExportFormat.QUANTIFY_SCHEDULER)
+
+            self.assertEqual(exported_schedule.name, "Exported OpenSquirrel circuit")
+
+            operation_ids = [v["operation_id"] for k, v in exported_schedule.schedulables.items()]
+            operations = [exported_schedule.operations[operation_id].name for operation_id in operation_ids]
+
+            self.assertEqual(
+                operations,
+                [
+                    "Rz(3.1415927, 'qreg[1]')",
+                    "Rxy(1.5707963, 1.5707963, 'qreg[1]')",
+                    "CZ (qreg[0], qreg[1])",
+                    "Rz(3.1415927, 'qreg[1]')",
+                    "Rxy(1.5707963, 1.5707963, 'qreg[1]')",
+                    "CZ (qreg[0], qreg[1])",
+                    "Rz(1.5707963, 'qreg[1]')",
+                    "Rxy(0.19634954, 1.5707963, 'qreg[1]')",
+                    "Rz(-1.5707963, 'qreg[1]')",
+                    "CZ (qreg[0], qreg[1])",
+                    "Rz(1.5707963, 'qreg[1]')",
+                    "Rxy(-0.19634954, 1.5707963, 'qreg[1]')",
+                    "Rz(-1.5707963, 'qreg[1]')",
+                    "CZ (qreg[0], qreg[1])",
+                    "Rz(0.19634954, 'qreg[0]')",
+                    "Rxy(-1.5707963, 1.5707963, 'qreg[0]')",
+                    "Rz(3.1415927, 'qreg[0]')",
+                    "Rz(3.1415927, 'qreg[1]')",
+                    "Rxy(1.5707963, 1.5707963, 'qreg[1]')",
+                ],
+            )
 
 
 if __name__ == "__main__":