Add all standard OpenQASM gates

graphix/ops.py

@@ -4,6 +4,7 @@
 
 from functools import reduce
 from itertools import product
+from math import pi
 from typing import TYPE_CHECKING, ClassVar, overload
 
 import numpy as np
@@ -23,6 +24,27 @@
     from graphix.parameter import ExpressionOrComplex
 
 
+@overload
+def controlled(gate: npt.NDArray[np.complex128]) -> npt.NDArray[np.complex128]: ...
+
+
+@overload
+def controlled(gate: npt.NDArray[np.object_]) -> npt.NDArray[np.object_]: ...
+
+
+def controlled(
+    gate: npt.NDArray[np.complex128] | npt.NDArray[np.object_],
+) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+    """Return the controlled version of a gate."""
+    n = gate.shape[0]
+    return np.block(
+        [
+            [np.eye(n), np.zeros((n, n))],
+            [np.zeros((n, n)), gate],
+        ]
+    )
+
+
 class Ops:
     """Basic single- and two-qubits operators."""
 
@@ -32,7 +54,21 @@ class Ops:
     Z: ClassVar[npt.NDArray[np.complex128]] = utils.lock(np.asarray([[1, 0], [0, -1]]))
     S: ClassVar[npt.NDArray[np.complex128]] = utils.lock(np.asarray([[1, 0], [0, 1j]]))
     SDG: ClassVar[npt.NDArray[np.complex128]] = utils.lock(np.asarray([[1, 0], [0, -1j]]))
+    T: ClassVar[npt.NDArray[np.complex128]] = utils.lock(np.asarray([[1, 0], [0, exp(1j * pi / 4)]]))
+    TDG: ClassVar[npt.NDArray[np.complex128]] = utils.lock(np.asarray([[1, 0], [0, exp(-1j * pi / 4)]]))
+    SX: ClassVar[npt.NDArray[np.complex128]] = utils.lock(1 / 2 * np.asarray([[1 + 1j, 1 - 1j], [1 - 1j, 1 + 1j]]))
+    SXDG: ClassVar[npt.NDArray[np.complex128]] = utils.lock(1 / 2 * np.asarray([[1 - 1j, 1 + 1j], [1 + 1j, 1 - 1j]]))
     H: ClassVar[npt.NDArray[np.complex128]] = utils.lock(np.asarray([[1, 1], [1, -1]]) / np.sqrt(2))
+    CY: ClassVar[npt.NDArray[np.complex128]] = utils.lock(
+        np.asarray(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, 0, -1j],
+                [0, 0, 1j, 0],
+            ],
+        )
+    )
     CZ: ClassVar[npt.NDArray[np.complex128]] = utils.lock(
         np.asarray(
             [
@@ -96,6 +132,31 @@ def _cast_array(
             return np.asarray(array, dtype=np.object_)
         return np.asarray(array, dtype=np.complex128)
 
+    @overload
+    @staticmethod
+    def p(theta: Angle) -> npt.NDArray[np.complex128]: ...
+
+    @overload
+    @staticmethod
+    def p(theta: Expression) -> npt.NDArray[np.object_]: ...
+
+    @staticmethod
+    def p(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        r"""Phase gate.
+
+        We have :math:`P(\theta) = \mathrm e^{\theta/2} R_Z(\theta)`.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 2*2 np.asarray
+        """
+        return Ops._cast_array([[1, 0], [0, exp(1j * angle_to_rad(theta))]], theta)
+
     @overload
     @staticmethod
     def rx(theta: Angle) -> npt.NDArray[np.complex128]: ...
@@ -167,6 +228,187 @@ def rz(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np
         """
         return Ops._cast_array([[exp(-1j * angle_to_rad(theta) / 2), 0], [0, exp(1j * angle_to_rad(theta) / 2)]], theta)
 
+    @staticmethod
+    def u(
+        theta: ParameterizedAngle, phi: ParameterizedAngle, lambda_: ParameterizedAngle
+    ) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Universal single-qubit gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+        phi : Angle | Expression
+            rotation angle in units of π
+        lambda_ : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 2*2 np.asarray
+        """
+        cos, sin = cos_sin(angle_to_rad(theta) / 2)
+        phi_rad = angle_to_rad(phi)
+        lambda_rad = angle_to_rad(lambda_)
+        return Ops._cast_array(
+            [[cos, -exp(1j * lambda_rad) * sin], [exp(1j * phi_rad) * sin, exp(1j * (phi_rad + lambda_rad)) * cos]],
+            theta,
+        )
+
+    @staticmethod
+    def cu(
+        theta: ParameterizedAngle, phi: ParameterizedAngle, lambda_: ParameterizedAngle, gamma: ParameterizedAngle
+    ) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Universal controlled single-qubit gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+        phi : Angle | Expression
+            rotation angle in units of π
+        lambda_ : Angle | Expression
+            rotation angle in units of π
+        gamma : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 4*4 np.asarray
+        """
+        cos, sin = cos_sin(angle_to_rad(theta) / 2)
+        phi_rad = angle_to_rad(phi)
+        lambda_rad = angle_to_rad(lambda_)
+        gamma_rad = angle_to_rad(gamma)
+        return Ops._cast_array(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, exp(1j * gamma_rad) * cos, -exp(1j * (gamma_rad + lambda_rad)) * sin],
+                [0, 0, exp(1j * (gamma_rad + phi_rad)) * sin, exp(1j * (gamma_rad + phi_rad + lambda_rad)) * cos],
+            ],
+            theta,
+        )
+
+    CH: ClassVar[npt.NDArray[np.complex128]] = controlled(H)
+
+    CSWAP: ClassVar[npt.NDArray[np.complex128]] = controlled(SWAP)
+
+    @overload
+    @staticmethod
+    def cj(theta: Angle) -> npt.NDArray[np.complex128]: ...
+
+    @overload
+    @staticmethod
+    def cj(theta: Expression) -> npt.NDArray[np.object_]: ...
+
+    @staticmethod
+    def cj(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Controlled-J gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 4*4 np.asarray
+        """
+        return controlled(Ops.j(theta))
+
+    @overload
+    @staticmethod
+    def cp(theta: Angle) -> npt.NDArray[np.complex128]: ...
+
+    @overload
+    @staticmethod
+    def cp(theta: Expression) -> npt.NDArray[np.object_]: ...
+
+    @staticmethod
+    def cp(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Controlled-phase gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 4*4 np.asarray
+        """
+        return controlled(Ops.p(theta))
+
+    @overload
+    @staticmethod
+    def crx(theta: Angle) -> npt.NDArray[np.complex128]: ...
+
+    @overload
+    @staticmethod
+    def crx(theta: Expression) -> npt.NDArray[np.object_]: ...
+
+    @staticmethod
+    def crx(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Controlled-RX gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 4*4 np.asarray
+        """
+        return controlled(Ops.rx(theta))
+
+    @overload
+    @staticmethod
+    def cry(theta: Angle) -> npt.NDArray[np.complex128]: ...
+
+    @overload
+    @staticmethod
+    def cry(theta: Expression) -> npt.NDArray[np.object_]: ...
+
+    @staticmethod
+    def cry(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Controlled-RY gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 4*4 np.asarray
+        """
+        return controlled(Ops.ry(theta))
+
+    @overload
+    @staticmethod
+    def crz(theta: Angle) -> npt.NDArray[np.complex128]: ...
+
+    @overload
+    @staticmethod
+    def crz(theta: Expression) -> npt.NDArray[np.object_]: ...
+
+    @staticmethod
+    def crz(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
+        """Controlled-RZ gate.
+
+        Parameters
+        ----------
+        theta : Angle | Expression
+            rotation angle in units of π
+
+        Returns
+        -------
+        operator : 4*4 np.asarray
+        """
+        return controlled(Ops.rz(theta))
+
     @overload
     @staticmethod
     def j(theta: Angle) -> npt.NDArray[np.complex128]: ...

graphix/qasm3_exporter.py

@@ -62,7 +62,7 @@ def circuit_to_qasm3_lines(circuit: Circuit, *, transpile: bool = True) -> Itera
         An iterator over the OpenQASM 3.0 lines that represent the circuit.
     """
     if transpile:
-        circuit = circuit.transpile_j_to_rzh().transpile_measurements_to_z_axis()
+        circuit = circuit.transpile_cj().transpile_rzz().transpile_j_to_rzh().transpile_measurements_to_z_axis()
     yield "OPENQASM 3;"
     yield 'include "stdgates.inc";'
     yield f"qubit[{circuit.width}] q;"
@@ -118,29 +118,57 @@ def instruction_to_qasm3(instruction: InstructionType) -> str:
                     "OpenQASM3 only supports measurements on Z axis. Use `Circuit.transpile_measurements_to_z_axis` to rewrite measurements on X and Y axes, or setting `transpile=True`."
                 )
             return f"b[{instruction.target}] = measure q[{instruction.target}]"
-        case InstructionKind.RX | InstructionKind.RY | InstructionKind.RZ:
+        case InstructionKind.RX | InstructionKind.RY | InstructionKind.RZ | InstructionKind.P:
             angle = angle_to_qasm3(instruction.angle)
             return qasm3_gate_call(
                 instruction.kind.name.lower(), args=[angle], operands=[qasm3_qubit(instruction.target)]
             )
+        case InstructionKind.CRX | InstructionKind.CRY | InstructionKind.CRZ | InstructionKind.CP:
+            angle = angle_to_qasm3(instruction.angle)
+            return qasm3_gate_call(
+                instruction.kind.name.lower(),
+                args=[angle],
+                operands=[qasm3_qubit(instruction.control), qasm3_qubit(instruction.target)],
+            )
         case InstructionKind.J:
             raise ValueError(
                 "J gates must be decomposed before QASM3 export using `Circuit.transpile_j_to_rzh`, or setting `transpile=True`."
             )
-        case InstructionKind.H | InstructionKind.S | InstructionKind.X | InstructionKind.Y | InstructionKind.Z:
+        case InstructionKind.CJ:
+            raise ValueError(
+                "CJ gates must be decomposed before QASM3 export using `Circuit.transpile_cj`, or setting `transpile=True`."
+            )
+        case (
+            InstructionKind.H
+            | InstructionKind.S
+            | InstructionKind.SDG
+            | InstructionKind.T
+            | InstructionKind.TDG
+            | InstructionKind.SX
+            | InstructionKind.SXDG
+            | InstructionKind.X
+            | InstructionKind.Y
+            | InstructionKind.Z
+        ):
             return qasm3_gate_call(instruction.kind.name.lower(), [qasm3_qubit(instruction.target)])
         case InstructionKind.I:
             return qasm3_gate_call("id", [qasm3_qubit(instruction.target)])
         case InstructionKind.CNOT:
             return qasm3_gate_call("cx", [qasm3_qubit(instruction.control), qasm3_qubit(instruction.target)])
+        case InstructionKind.CY:
+            return qasm3_gate_call("cy", [qasm3_qubit(instruction.control), qasm3_qubit(instruction.target)])
         case InstructionKind.SWAP:
             return qasm3_gate_call("swap", [qasm3_qubit(instruction.targets[i]) for i in (0, 1)])
+        case InstructionKind.CSWAP:
+            return qasm3_gate_call(
+                "cswap",
+                [qasm3_qubit(qubit) for qubit in [instruction.control, *[instruction.targets[i] for i in (0, 1)]]],
+            )
         case InstructionKind.CZ:
             return qasm3_gate_call("cz", [qasm3_qubit(instruction.targets[i]) for i in (0, 1)])
         case InstructionKind.RZZ:
-            angle = angle_to_qasm3(instruction.angle)
-            return qasm3_gate_call(
-                "crz", args=[angle], operands=[qasm3_qubit(instruction.control), qasm3_qubit(instruction.target)]
+            raise ValueError(
+                "RZZ gates must be decomposed before QASM3 export using `Circuit.transpile_rzz`, or setting `transpile=True`."
             )
         case InstructionKind.CCX:
             return qasm3_gate_call(
@@ -151,6 +179,23 @@ def instruction_to_qasm3(instruction: InstructionType) -> str:
                     qasm3_qubit(instruction.target),
                 ],
             )
+        case InstructionKind.U:
+            theta = angle_to_qasm3(instruction.theta)
+            phi = angle_to_qasm3(instruction.phi)
+            lambda_ = angle_to_qasm3(instruction.lambda_)
+            return qasm3_gate_call("U", args=[theta, phi, lambda_], operands=[qasm3_qubit(instruction.target)])
+        case InstructionKind.CU:
+            theta = angle_to_qasm3(instruction.theta)
+            phi = angle_to_qasm3(instruction.phi)
+            lambda_ = angle_to_qasm3(instruction.lambda_)
+            gamma = angle_to_qasm3(instruction.gamma)
+            return qasm3_gate_call(
+                "cu",
+                args=[theta, phi, lambda_, gamma],
+                operands=[qasm3_qubit(instruction.control), qasm3_qubit(instruction.target)],
+            )
+        case InstructionKind.GPHASE:
+            return qasm3_gate_call("gphase", operands=[], args=[angle_to_qasm3(instruction.angle)])
         case _:
             assert_never(instruction.kind)