[CQT-297] Check that composing bloch sphere rotations is done correctly (#425)

Co-authored-by: Roberto Turrado Camblor <rturrado@gmail.com>

Author: elenbaasc

opensquirrel/passes/decomposer/cnot_decomposer.py

@@ -38,7 +38,7 @@ def decompose(self, g: Gate) -> list[Gate]:
         # See https://threeplusone.com/pubs/on_gates.pdf
 
         # Try special case first, see https://arxiv.org/pdf/quant-ph/9503016.pdf lemma 5.5
-        controlled_rotation_times_x = general_merger.compose_bloch_sphere_rotations(X(target_qubit), g.target_gate)
+        controlled_rotation_times_x = general_merger.compose_bloch_sphere_rotations(g.target_gate, X(target_qubit))
         theta0_with_x, theta1_with_x, theta2_with_x = ZYZDecomposer().get_decomposition_angles(
             controlled_rotation_times_x.axis,
             controlled_rotation_times_x.angle,

opensquirrel/passes/merger/general_merger.py

@@ -15,8 +15,16 @@
 
 def compose_bloch_sphere_rotations(bsr_a: BlochSphereRotation, bsr_b: BlochSphereRotation) -> BlochSphereRotation:
     """Computes the Bloch sphere rotation resulting from the composition of two Bloch sphere rotations.
-    The first rotation is applied and then the second.
-    If the final Bloch sphere rotation is anonymous, we will try to match it to a default gate.
+    The first rotation (A) is applied and then the second (B):
+
+    As separate gates:
+        A q
+        B q
+
+    A linear operations:
+        (B * A) q
+
+    If the final Bloch sphere rotation is anonymous, we try to match it to a default gate.
 
     Uses Rodrigues' rotation formula (see https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula).
     """
@@ -40,7 +48,7 @@ def compose_bloch_sphere_rotations(bsr_a: BlochSphereRotation, bsr_b: BlochSpher
             * (
                 sin(bsr_a.angle / 2) * cos(bsr_b.angle / 2) * bsr_a.axis.value
                 + cos(bsr_a.angle / 2) * sin(bsr_b.angle / 2) * bsr_b.axis.value
-                + sin(bsr_a.angle / 2) * sin(bsr_b.angle / 2) * np.cross(bsr_a.axis, bsr_b.axis)
+                + sin(bsr_a.angle / 2) * sin(bsr_b.angle / 2) * np.cross(bsr_b.axis, bsr_a.axis)
             )
         ),
         order_of_magnitude,

opensquirrel/passes/merger/single_qubit_gates_merger.py

@@ -25,7 +25,7 @@ def merge(self, ir: IR, qubit_register_size: int) -> None:
             instruction: Instruction = cast(Instruction, statement)
             if isinstance(instruction, BlochSphereRotation):
                 already_accumulated = accumulators_per_qubit[instruction.qubit]
-                composed = compose_bloch_sphere_rotations(instruction, already_accumulated)
+                composed = compose_bloch_sphere_rotations(already_accumulated, instruction)
                 accumulators_per_qubit[instruction.qubit] = composed
                 del ir.statements[statement_index]
                 continue

test/merger/test_general_merger.py

@@ -2,7 +2,7 @@
 
 import pytest
 
-from opensquirrel import Circuit, CircuitBuilder, H, I, Rx, Ry, X
+from opensquirrel import Circuit, CircuitBuilder, H, I, Rx, Ry, X, Y, Z
 from opensquirrel.ir import BlochSphereRotation, Float
 from opensquirrel.passes.merger.general_merger import compose_bloch_sphere_rotations, rearrange_barriers
 
@@ -19,29 +19,45 @@ def test_compose_bloch_sphere_rotations_different_axis() -> None:
     a = BlochSphereRotation(qubit=123, axis=(1, 0, 0), angle=math.pi / 2)
     b = BlochSphereRotation(qubit=123, axis=(0, 0, 1), angle=-math.pi / 2)
     c = BlochSphereRotation(qubit=123, axis=(0, 1, 0), angle=math.pi / 2)
-    composed = compose_bloch_sphere_rotations(a, compose_bloch_sphere_rotations(b, c))
+    composed = compose_bloch_sphere_rotations(compose_bloch_sphere_rotations(c, b), a)
     assert composed == BlochSphereRotation(qubit=123, axis=(1, 1, 0), angle=math.pi)
 
 
 @pytest.mark.parametrize(
     ("bsr_a", "bsr_b", "expected_result"),
     [
+        (Y(0), X(0), BlochSphereRotation(0, axis=(0, 0, 1), angle=math.pi, phase=math.pi)),
+        (X(0), Y(0), BlochSphereRotation(0, axis=(0, 0, -1), angle=math.pi, phase=math.pi)),
+        (Z(0), Y(0), BlochSphereRotation(0, axis=(1, 0, 0), angle=math.pi, phase=math.pi)),
+        (Y(0), Z(0), BlochSphereRotation(0, axis=(-1, 0, 0), angle=math.pi, phase=math.pi)),
+        (Z(0), X(0), BlochSphereRotation(0, axis=(0, -1, 0), angle=math.pi, phase=math.pi)),
+        (X(0), Z(0), BlochSphereRotation(0, axis=(0, 1, 0), angle=math.pi, phase=math.pi)),
+        (
+            Rx(0, theta=math.pi),
+            Ry(0, theta=-math.pi / 2),
+            BlochSphereRotation(0, axis=(0.70711, 0.0, 0.70711), angle=math.pi, phase=0.0),
+        ),
         (I(0), Rx(0, theta=math.pi), Rx(0, theta=math.pi)),
         (Rx(0, theta=math.pi), I(0), Rx(0, theta=math.pi)),
         (X(0), X(0), I(0)),
         (H(0), H(0), I(0)),
         (Rx(0, theta=math.pi), Rx(0, theta=math.pi), I(0)),
         (Rx(0, theta=math.pi / 2), Rx(0, theta=math.pi / 2), Rx(0, theta=math.pi)),
-        (Rx(0, theta=math.pi), Ry(0, theta=math.pi / 2), BlochSphereRotation(0, axis=(1, 0, 1), angle=math.pi)),
     ],
     ids=[
+        "[bsr_a = Y, bsr_b = X] X * Y = iZ",
+        "[bsr_a = X, bsr_b = Y] Y * X = -iZ",
+        "[bsr_a = Z, bsr_b = Y] Y * Z = iX",
+        "[bsr_a = Y, bsr_b = Z] Z * Y = -iX",
+        "[bsr_a = Z, bsr_b = X] X * Z = -iY",
+        "[bsr_a = X, bsr_b = Z] Z * X = iY",
+        "[bsr_a = Rx(pi), bsr_b = Ry(-pi/2)] Ry(-pi/2) * Rx(pi) ~ H",
         "[bsr_a == I]",
         "[bsr_b == I]",
         "[bsr_a.generator == bsr_b.generator] X * X == I",
         "[bsr_a.generator == bsr_b.generator] H * H == I",
         "[bsr_a.generator == bsr_b.generator] Rx(pi) * Rx(pi) == I",
         "[bsr_a.generator == bsr_b.generator] Rx(pi/2) * Rx(pi/2) = Rx(pi) ~ X",
-        "[bsr_a.generator != bsr_b.generator] Rx(pi) * Ry(pi/2) ~ H",
     ],
 )
 def test_compose_bloch_sphere_rotations(