From bd3a574de6c52176d958c4079c635a563f693283 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 25 Aug 2026 18:05:46 +0200
Subject: [PATCH 01/18] Add all standard OpenQASM gates

This commit adds the following gates specified in the OpenQASM 3.0
specification: SDG, T, TDG, SX, SXDG, CY, P, U, CP, CRX, CRY, CRZ,
CU, CSWAP, and GPHASE.

Support for GPHASE is currently limited: the gate is ignored by the
circuit simulator and the transpiler. It is useful, however, to have
it internally in order to add control to the decompositions of the
other gates.

Additionally, this commit adds the CJ gate, which is useful internally
for adding control to decompositions that involve J gates.

These additions were motivated by the issue
TeamGraphix/graphix-qasm-parser#14. The OpenQASM parser plugin can now
parse all the gates and map them to native Graphix gates without
translation, except for `u1`, `u2`, and `u3`, which are maintained for
compatibility with OpenQASM 2.0 and are translated.
---
 graphix/instruction.py                        | 514 +++++++++++++++-
 graphix/ops.py                                | 242 ++++++++
 graphix/qasm3_exporter.py                     |  57 +-
 graphix/transpiler.py                         | 578 +++++++++++++++++-
 tests/test_instruction.py                     | 196 ++++--
 .../test_qasm3_exporter_to_graphix_parser.py  |  23 +-
 tests/test_qasm3_exporter_to_qiskit.py        |  48 +-
 tests/test_transpiler.py                      | 136 ++++-
 8 files changed, 1644 insertions(+), 150 deletions(-)

diff --git a/graphix/instruction.py b/graphix/instruction.py
index c70d58968..eb5927024 100644
--- a/graphix/instruction.py
+++ b/graphix/instruction.py
@@ -55,6 +55,22 @@ class InstructionKind(Enum):
     RX = enum.auto()
     RY = enum.auto()
     RZ = enum.auto()
+    SDG = enum.auto()
+    T = enum.auto()
+    TDG = enum.auto()
+    SX = enum.auto()
+    SXDG = enum.auto()
+    CY = enum.auto()
+    P = enum.auto()
+    U = enum.auto()
+    CJ = enum.auto()
+    CP = enum.auto()
+    CRX = enum.auto()
+    CRY = enum.auto()
+    CRZ = enum.auto()
+    CU = enum.auto()
+    CSWAP = enum.auto()
+    GPHASE = enum.auto()
 
 
 class _KindChecker:
@@ -159,24 +175,39 @@ def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> RZZ:
 
 
 @dataclass(repr=False)
-class CNOT(_KindChecker, BaseInstruction):
-    """CNOT circuit instruction."""
+class ControlledSingleTargetInstruction(BaseInstruction):
+    """Base class for controlled single-target circuit instructions."""
 
     target: int
     control: int
-    kind: ClassVar[Literal[InstructionKind.CNOT]] = field(default=InstructionKind.CNOT, init=False)
 
     @override
-    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> CNOT:
+    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
         target = visitor.visit_qubit(self.target)
         control = visitor.visit_qubit(self.control)
         if copy:
-            return CNOT(target, control)
+            return type(self)(target, control)
         self.target = target
         self.control = control
         return self
 
 
+@dataclass(repr=False)
+class CY(_KindChecker, ControlledSingleTargetInstruction):
+    """CY circuit instruction."""
+
+    kind: ClassVar[Literal[InstructionKind.CY]] = field(default=InstructionKind.CY, init=False)
+
+
+@dataclass(repr=False)
+class CNOT(_KindChecker, ControlledSingleTargetInstruction):
+    """CNOT circuit instruction."""
+
+    kind: ClassVar[Literal[InstructionKind.CNOT]] = field(default=InstructionKind.CNOT, init=False)
+
+
+# CZ is not defined as a ControlledSingleTargetInstruction because of
+# the symmetry between the control and the target.
 @dataclass(repr=False)
 class CZ(_KindChecker, BaseInstruction):
     """CZ circuit instruction."""
@@ -211,6 +242,51 @@ def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> SWAP:
         return self
 
 
+@dataclass(repr=False)
+class CSWAP(_KindChecker, BaseInstruction):
+    r"""CSWAP circuit instruction.
+
+    The CSWAP gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & \cos \frac \theta 2 & -\mathrm i \sin \frac \theta 2\\
+        0 & 0 & -\mathrm i \sin \frac \theta 2 & \cos \frac \theta 2
+      \end{matrix}\right]
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\
+        0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\
+        0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
+        0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\
+        0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\
+        0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\
+        0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
+      \end{matrix}\right]
+    """
+
+    control: int
+    targets: tuple[int, int]
+    kind: ClassVar[Literal[InstructionKind.CSWAP]] = field(default=InstructionKind.CSWAP, init=False)
+
+    @override
+    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> CSWAP:
+        control = visitor.visit_qubit(self.control)
+        u, v = self.targets
+        targets = (visitor.visit_qubit(u), visitor.visit_qubit(v))
+        if copy:
+            return CSWAP(control, targets)
+        self.control = control
+        self.targets = targets
+        return self
+
+
 @dataclass(repr=False)
 class SingleTargetInstruction(BaseInstruction):
     """Base class for single-target circuit instructions."""
@@ -240,6 +316,71 @@ class S(_KindChecker, SingleTargetInstruction):
     kind: ClassVar[Literal[InstructionKind.S]] = field(default=InstructionKind.S, init=False)
 
 
+@dataclass(repr=False)
+class SDG(_KindChecker, SingleTargetInstruction):
+    r"""SDG circuit instruction.
+
+    The :math:`S^\dagger` gate applies the matrix
+    :math:`\left[\begin{matrix}1 & 0\\0 & - \mathrm i\end{matrix}\right]`.
+
+    We have :math:`S^\dagger = \mathrm e^{\mathrm i \frac \pi 4} R_Z(-\frac \pi 2)`.
+    """
+
+    kind: ClassVar[Literal[InstructionKind.SDG]] = field(default=InstructionKind.SDG, init=False)
+
+
+@dataclass(repr=False)
+class T(_KindChecker, SingleTargetInstruction):
+    r"""T circuit instruction.
+
+    The :math:`T` gate applies the matrix
+    :math:`\left[\begin{matrix}1 & 0\\0 & \mathrm e^{\mathrm i \frac \pi 4}\end{matrix}\right]`.
+
+    We have :math:`T = \mathrm e^{\mathrm i \frac \pi 8} R_Z(\frac \pi 4)`.
+    """
+
+    kind: ClassVar[Literal[InstructionKind.T]] = field(default=InstructionKind.T, init=False)
+
+
+@dataclass(repr=False)
+class TDG(_KindChecker, SingleTargetInstruction):
+    r"""TDG circuit instruction.
+
+    The :math:`T^\dagger` gate applies the matrix
+    :math:`\left[\begin{matrix}1 & 0\\0 & \mathrm e^{- \mathrm i \frac \pi 4}\end{matrix}\right]`.
+
+    We have :math:`T^\dagger = \mathrm e^{\mathrm i \frac \pi 8} R_Z(- \frac \pi 4)`.
+    """
+
+    kind: ClassVar[Literal[InstructionKind.TDG]] = field(default=InstructionKind.TDG, init=False)
+
+
+@dataclass(repr=False)
+class SX(_KindChecker, SingleTargetInstruction):
+    r"""SX circuit instruction.
+
+    The :math:`SX` (:math:`\sqrt X`) gate applies the matrix
+    :math:`\frac 1 2 \left[\begin{matrix}1 + \mathrm i & 1 - \mathrm i\\1 - \mathrm i & 1 + \mathrm i\end{matrix}\right]`.
+
+    We have :math:`SX = \mathrm e^{\mathrm i \frac \pi 4} R_X(\frac \pi 2)`.
+    """
+
+    kind: ClassVar[Literal[InstructionKind.SX]] = field(default=InstructionKind.SX, init=False)
+
+
+@dataclass(repr=False)
+class SXDG(_KindChecker, SingleTargetInstruction):
+    r"""SXDG circuit instruction.
+
+    The :math:`SX^\dagger` (:math:`{\sqrt X}^\dagger`) gate applies the matrix
+    :math:`\frac 1 2 \left[\begin{matrix}1 - \mathrm i & 1 + \mathrm i\\1 + \mathrm i & 1 - \mathrm i\end{matrix}\right]`.
+
+    We have :math:`SX^\dagger = \mathrm e^{\mathrm i \frac \pi 4} R_X(-\frac \pi 2)`.
+    """
+
+    kind: ClassVar[Literal[InstructionKind.SXDG]] = field(default=InstructionKind.SXDG, init=False)
+
+
 @dataclass(repr=False)
 class X(_KindChecker, SingleTargetInstruction):
     """X circuit instruction."""
@@ -305,6 +446,25 @@ def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
         return self
 
 
+@dataclass(repr=False)
+class P(_KindChecker, RotationInstruction):
+    r"""P rotation circuit instruction.
+
+    The :math:`P(\theta)` gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0\\
+        0 & \mathrm e^{\mathrm i \theta}
+      \end{matrix}\right]
+
+    We have :math:`P(\theta) = \mathrm e^{\theta/2} R_Z(\theta)`.
+    """
+
+    kind: ClassVar[Literal[InstructionKind.P]] = field(default=InstructionKind.P, init=False)
+
+
 @dataclass(repr=False)
 class RX(_KindChecker, RotationInstruction):
     """X rotation circuit instruction."""
@@ -333,36 +493,265 @@ class J(_KindChecker, RotationInstruction):
     kind: ClassVar[Literal[InstructionKind.J]] = field(default=InstructionKind.J, init=False)
 
 
-class InstructionWithoutRZZ:
-    """Grouping of all instructions except RZZ for namespace exposure.
+@dataclass(repr=False)
+class U(_KindChecker, BaseInstruction):
+    r"""U circuit instruction.
 
-    Notes
-    -----
-    This class is not meant to be instantiated, but rather serves as a namespace for all instructions except RZZ.
-    The type alias for "any command" is :data:`InstructionKind`.
+    The :math:`U(\theta, \phi, \lambda)` gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        \cos \frac \theta 2 & - \mathrm e^{\mathrm i\lambda} \sin \frac \theta 2 \\
+        \mathrm e^{\mathrm i \phi} \sin \frac \theta 2 & \mathrm e^{\mathrm i (\phi + \lambda)} \cos \frac \theta 2
+      \end{matrix}\right]
+
+    It can be decomposed as
+
+    .. math::
+
+      U(\theta, \phi, \lambda) = \mathrm e^{\mathrm i\frac{\phi + \lambda}{2}} R_Z(\phi) R_Y(\theta) R_Z(\lambda)
+      = \mathrm e^{\mathrm i\frac{\theta}{2}}
+        H J\left(\phi + \frac{\pi}{2}\right)
+        J(\theta)
+        J\left(\lambda - \frac{\pi}{2}\right)
     """
 
-    CCX: TypeAlias = CCX
-    CNOT: TypeAlias = CNOT
-    CZ: TypeAlias = CZ
-    SWAP: TypeAlias = SWAP
-    H: TypeAlias = H
-    S: TypeAlias = S
-    X: TypeAlias = X
-    Y: TypeAlias = Y
-    Z: TypeAlias = Z
-    I: TypeAlias = I
-    M: TypeAlias = M
-    RX: TypeAlias = RX
-    RY: TypeAlias = RY
-    RZ: TypeAlias = RZ
-    J: TypeAlias = J
+    target: int
+    theta: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    phi: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    lambda_: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    kind: ClassVar[Literal[InstructionKind.U]] = field(default=InstructionKind.U, init=False)
 
-    def __init__(self) -> None:
-        raise TypeError("InstructionWithoutRZZ is a namespace, not a class.")
+    @override
+    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
+        target = visitor.visit_qubit(self.target)
+        theta = visitor.visit_angle(self.theta)
+        phi = visitor.visit_angle(self.phi)
+        lambda_ = visitor.visit_angle(self.lambda_)
+        if copy:
+            return type(self)(target, theta, phi, lambda_)
+        self.target = target
+        self.theta = theta
+        self.phi = phi
+        self.lambda_ = lambda_
+        return self
 
 
-class Instruction(InstructionWithoutRZZ):
+@dataclass(repr=False)
+class CU(_KindChecker, BaseInstruction):
+    r"""Controlled-U circuit instruction.
+
+    The :math:`CU(\theta, \phi, \lambda, \gamma)` gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0 \\
+        0 & 1 & 0 & 0 \\
+        0 & 0 & \mathrm e^{\mathrm i\gamma}
+          \cos\left(\frac{\theta}{2}\right) &
+          -\mathrm e^{\mathrm i(\gamma + \lambda)}
+          \sin\left(\frac{\theta}{2}\right) \\
+        0 & 0 & \mathrm e^{\mathrm i(\gamma + \phi)}
+          \sin\left(\frac{\theta}{2}\right) &
+          \mathrm e^{\mathrm i(\gamma + \phi + \lambda)}
+          \cos\left(\frac{\theta}{2}\right)
+      \end{matrix}\right]
+
+    It can be decomposed as
+
+    .. math::
+
+      CU(\theta, \phi, \lambda, \gamma) =
+        \left(P\left(\frac{\gamma - \theta} 2\right) \otimes I)
+        CJ(0) CJ\left(\phi + \frac \pi 2\right)
+        CJ(\theta) CJ\left(\lambda - \frac \pi 2\right)
+    """
+
+    control: int
+    target: int
+    theta: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    phi: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    lambda_: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    gamma: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    kind: ClassVar[Literal[InstructionKind.CU]] = field(default=InstructionKind.CU, init=False)
+
+    @override
+    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
+        control = visitor.visit_qubit(self.control)
+        target = visitor.visit_qubit(self.target)
+        theta = visitor.visit_angle(self.theta)
+        phi = visitor.visit_angle(self.phi)
+        lambda_ = visitor.visit_angle(self.lambda_)
+        gamma = visitor.visit_angle(self.gamma)
+        if copy:
+            return type(self)(control, target, theta, phi, lambda_, gamma)
+        self.control = control
+        self.target = target
+        self.theta = theta
+        self.phi = phi
+        self.lambda_ = lambda_
+        self.gamma = gamma
+        return self
+
+
+@dataclass(repr=False)
+class ControlledRotationInstruction(BaseInstruction):
+    """Base class for rotation instructions."""
+
+    target: int
+    control: int
+    angle: ParameterizedAngle = field(metadata={"repr": repr_angle})
+
+    @override
+    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
+        target = visitor.visit_qubit(self.target)
+        control = visitor.visit_qubit(self.control)
+        angle = visitor.visit_angle(self.angle)
+        if copy:
+            return type(self)(target, control, angle)
+        self.target = target
+        self.control = control
+        self.angle = angle
+        return self
+
+
+@dataclass(repr=False)
+class CP(_KindChecker, ControlledRotationInstruction):
+    r"""Controlled-P rotation circuit instruction.
+
+    The :math:`CP(\theta)` gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & 1 & 0\\
+        0 & 0 & 0 & \mathrm e^{\mathrm i \theta}
+      \end{matrix}\right]
+    """
+
+    kind: ClassVar[Literal[InstructionKind.CP]] = field(default=InstructionKind.CP, init=False)
+
+
+@dataclass(repr=False)
+class CRX(_KindChecker, ControlledRotationInstruction):
+    r"""Controlled-X rotation circuit instruction.
+
+    The :math:`CRX(\theta)` gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & \cos \frac \theta 2 & -\mathrm i \sin \frac \theta 2\\
+        0 & 0 & -\mathrm i \sin \frac \theta 2 & \cos \frac \theta 2
+      \end{matrix}\right]
+    """
+
+    kind: ClassVar[Literal[InstructionKind.CRX]] = field(default=InstructionKind.CRX, init=False)
+
+
+@dataclass(repr=False)
+class CRY(_KindChecker, ControlledRotationInstruction):
+    r"""Controlled-Y rotation circuit instruction.
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & \cos \frac \theta 2 & - \sin \frac \theta 2\\
+        0 & 0 & \sin \frac \theta 2 & \cos \frac \theta 2
+      \end{matrix}\right]
+    """
+
+    kind: ClassVar[Literal[InstructionKind.CRY]] = field(default=InstructionKind.CRY, init=False)
+
+
+@dataclass(repr=False)
+class CRZ(_KindChecker, ControlledRotationInstruction):
+    r"""Controlled-Z rotation circuit instruction.
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & \mathrm e^{-\mathrm i \frac \theta 2} & 0\\
+        0 & 0 & 0 & \mathrm e^{\mathrm i \frac \theta 2}
+      \end{matrix}\right]
+    """
+
+    kind: ClassVar[Literal[InstructionKind.CRZ]] = field(default=InstructionKind.CRZ, init=False)
+
+
+@dataclass(repr=False)
+class CJ(_KindChecker, ControlledRotationInstruction):
+    r"""Controlled-J circuit instruction.
+
+    The :math:`CJ(\alpha)` gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & \frac 1 {\sqrt 2} & \frac 1 {\sqrt 2} \mathrm e^{\mathrm i \alpha}\\
+        0 & 0 & \frac 1 {\sqrt 2} & - \frac 1 {\sqrt 2} \mathrm e^{\mathrm i \alpha}
+      \end{matrix}\right]
+
+    Following Lemmas 4.3 and 5.1 of Barenco et al. (1995), we define:
+
+    .. math::
+
+      \begin{aligned}
+        A &= R_Y\left(\frac \pi 4\right),\\
+        B &= R_Y\left(- \frac \pi 4\right) R_Z(- \delta),\\
+        C &= R_Z(\delta),\\
+        \delta &= \frac {\alpha + \pi} 2
+      \end{aligned}
+
+    These operators satisfy :math:`ABC = I` and
+    :math:`AXBXC = \mathrm e^{-\mathrm i \delta} J(\alpha)` with
+    :math:``.
+
+    Consequently, :math:`CJ(\alpha)` can be decomposed as:
+
+    .. math::
+
+      CJ(\alpha) = (P(\delta) \otimes I) \, (I \otimes A) \, CX \, (I \otimes B) \, CX \, (I \otimes C)
+
+    References
+    ----------
+    Barenco, A., Bennett, C. H., Cleve, R., DiVincenzo, D. P., Margolus, N., Shor, P., Sleator, T., Smolin, J. A., & Weinfurter, H. (1995).
+    Elementary gates for quantum computation. Physical Review A, 52(5), 3457-3467.
+    https://doi.org/10.1103/physreva.52.3457
+    """
+
+    kind: ClassVar[Literal[InstructionKind.CJ]] = field(default=InstructionKind.CJ, init=False)
+
+
+@dataclass(repr=False)
+class GPHASE(_KindChecker, BaseInstruction):
+    """GPHASE circuit instruction."""
+
+    angle: ParameterizedAngle = field(metadata={"repr": repr_angle})
+    kind: ClassVar[Literal[InstructionKind.GPHASE]] = field(default=InstructionKind.GPHASE, init=False)
+
+    @override
+    def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> GPHASE:
+        angle = visitor.visit_angle(self.angle)
+        if copy:
+            return GPHASE(angle)
+        self.angle = angle
+        return self
+
+
+class Instruction:
     """Grouping of all instructions for namespace exposure.
 
     Notes
@@ -371,12 +760,75 @@ class Instruction(InstructionWithoutRZZ):
     The type alias for "any command" is :data:`InstructionKind`.
     """
 
+    I: TypeAlias = I
+    X: TypeAlias = X
+    Y: TypeAlias = Y
+    Z: TypeAlias = Z
+    H: TypeAlias = H
+    S: TypeAlias = S
+    SDG: TypeAlias = SDG
+    T: TypeAlias = T
+    TDG: TypeAlias = TDG
+    SX: TypeAlias = SX
+    SXDG: TypeAlias = SXDG
+    J: TypeAlias = J
+    P: TypeAlias = P
+    RX: TypeAlias = RX
+    RY: TypeAlias = RY
+    RZ: TypeAlias = RZ
+    U: TypeAlias = U
+    CJ: TypeAlias = CJ
+    CP: TypeAlias = CP
+    CRX: TypeAlias = CRX
+    CRY: TypeAlias = CRY
+    CRZ: TypeAlias = CRZ
+    CU: TypeAlias = CU
+    CNOT: TypeAlias = CNOT
+    CY: TypeAlias = CY
+    CZ: TypeAlias = CZ
+    CCX: TypeAlias = CCX
     RZZ: TypeAlias = RZZ
+    SWAP: TypeAlias = SWAP
+    CSWAP: TypeAlias = CSWAP
+    M: TypeAlias = M
+    GPHASE: TypeAlias = GPHASE
 
     def __init__(self) -> None:
         raise TypeError("Instruction is a namespace, not a class.")
 
 
 if TYPE_CHECKING:
-    InstructionTypeWithoutRZZ = CCX | CNOT | SWAP | CZ | H | S | X | Y | Z | I | M | RX | RY | RZ | J
-    InstructionType = InstructionTypeWithoutRZZ | RZZ
+    InstructionType = (
+        I
+        | X
+        | Y
+        | Z
+        | H
+        | S
+        | SDG
+        | T
+        | TDG
+        | SX
+        | SXDG
+        | J
+        | P
+        | RX
+        | RY
+        | RZ
+        | U
+        | CJ
+        | CP
+        | CRX
+        | CRY
+        | CRZ
+        | CU
+        | CNOT
+        | CY
+        | CZ
+        | CCX
+        | RZZ
+        | SWAP
+        | CSWAP
+        | M
+        | GPHASE
+    )
diff --git a/graphix/ops.py b/graphix/ops.py
index 381a212df..5cfcecb7f 100644
--- a/graphix/ops.py
+++ b/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]: ...
diff --git a/graphix/qasm3_exporter.py b/graphix/qasm3_exporter.py
index 19d152f64..c6a8b4778 100644
--- a/graphix/qasm3_exporter.py
+++ b/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", [angle_to_qasm3(instruction.angle)])
         case _:
             assert_never(instruction.kind)
 
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 9613b9e9d..d3b3a972b 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -17,7 +17,7 @@
 # override introduced in Python 3.12
 from typing_extensions import assert_never, override
 
-from graphix import command, instruction, parameter
+from graphix import Instruction, command, instruction, parameter
 from graphix.branch_selector import BranchSelector, RandomBranchSelector
 from graphix.flow.core import CausalFlow, _corrections_to_partial_order_layers
 from graphix.fundamentals import ANGLE_PI, Axis
@@ -170,6 +170,38 @@ def add(self, instr: InstructionType) -> None:
                 self.rz(instr.target, instr.angle)
             case InstructionKind.J:
                 self.j(instr.target, instr.angle)
+            case InstructionKind.SDG:
+                self.sdg(instr.target)
+            case InstructionKind.T:
+                self.t(instr.target)
+            case InstructionKind.TDG:
+                self.tdg(instr.target)
+            case InstructionKind.SX:
+                self.sx(instr.target)
+            case InstructionKind.SXDG:
+                self.sxdg(instr.target)
+            case InstructionKind.CY:
+                self.cy(instr.control, instr.target)
+            case InstructionKind.P:
+                self.p(instr.target, instr.angle)
+            case InstructionKind.U:
+                self.u(instr.target, instr.theta, instr.phi, instr.lambda_)
+            case InstructionKind.CJ:
+                self.cj(instr.control, instr.target, instr.angle)
+            case InstructionKind.CP:
+                self.cp(instr.control, instr.target, instr.angle)
+            case InstructionKind.CRX:
+                self.crx(instr.control, instr.target, instr.angle)
+            case InstructionKind.CRY:
+                self.cry(instr.control, instr.target, instr.angle)
+            case InstructionKind.CRZ:
+                self.crz(instr.control, instr.target, instr.angle)
+            case InstructionKind.CU:
+                self.cu(instr.control, instr.target, instr.theta, instr.phi, instr.lambda_, instr.gamma)
+            case InstructionKind.CSWAP:
+                self.cswap(instr.control, instr.targets[0], instr.targets[1])
+            case InstructionKind.GPHASE:
+                self.gphase(instr.angle)
             case _:
                 assert_never(instr.kind)
 
@@ -427,6 +459,310 @@ def m(self, qubit: int, axis: Axis) -> None:
         self.instruction.append(instruction.M(target=qubit, axis=axis))
         self.active_qubits.remove(qubit)
 
+    def sdg(self, qubit: int) -> None:
+        """Apply an SDG gate.
+
+        See :class:`~graphix.instruction.SDG` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.SDG(target=qubit))
+
+    def t(self, qubit: int) -> None:
+        """Apply a T gate.
+
+        See :class:`~graphix.instruction.T` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.T(target=qubit))
+
+    def tdg(self, qubit: int) -> None:
+        """Apply a TDG gate.
+
+        See :class:`~graphix.instruction.TDG` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.TDG(target=qubit))
+
+    def sx(self, qubit: int) -> None:
+        """Apply an SX gate.
+
+        See :class:`~graphix.instruction.SX` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.SX(target=qubit))
+
+    def sxdg(self, qubit: int) -> None:
+        """Apply an SXDG gate.
+
+        See :class:`~graphix.instruction.SXDG` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.SXDG(target=qubit))
+
+    def cy(self, control: int, target: int) -> None:
+        """Apply a Controlled-Y gate.
+
+        See :class:`~graphix.instruction.CY` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(Instruction.CY(control=control, target=target))
+
+    def p(self, qubit: int, angle: ParameterizedAngle) -> None:
+        """Apply a Phase rotation gate.
+
+        See :class:`~graphix.instruction.P` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.P(target=qubit, angle=angle))
+
+    def u(self, qubit: int, theta: ParameterizedAngle, phi: ParameterizedAngle, lambda_: ParameterizedAngle) -> None:
+        """Apply an U gate.
+
+        See :class:`~graphix.instruction.U` for more information.
+
+        Parameters
+        ----------
+        qubit : int
+            target qubit
+        theta : ParameterizedAngle
+            rotation angle in units of π
+        phi : ParameterizedAngle
+            rotation angle in units of π
+        lambda_ : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert qubit in self.active_qubits
+        self.instruction.append(Instruction.U(target=qubit, theta=theta, phi=phi, lambda_=lambda_))
+
+    def cj(self, control: int, target: int, angle: ParameterizedAngle) -> None:
+        """Apply a controlled-J rotation gate.
+
+        See :class:`~graphix.instruction.CJ` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(Instruction.CJ(control=control, target=target, angle=angle))
+
+    def cp(self, control: int, target: int, angle: ParameterizedAngle) -> None:
+        """Apply a controlled-P rotation gate.
+
+        See :class:`~graphix.instruction.CP` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(Instruction.CP(control=control, target=target, angle=angle))
+
+    def crx(self, control: int, target: int, angle: ParameterizedAngle) -> None:
+        """Apply an controlled-X rotation gate.
+
+        See :class:`~graphix.instruction.CRX` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(Instruction.CRX(control=control, target=target, angle=angle))
+
+    def cry(self, control: int, target: int, angle: ParameterizedAngle) -> None:
+        """Apply a controlled-Y rotation gate.
+
+        See :class:`~graphix.instruction.CRY` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        angle : ParameterizedAngle
+            angle in units of π
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(Instruction.CRY(control=control, target=target, angle=angle))
+
+    def crz(self, control: int, target: int, angle: ParameterizedAngle) -> None:
+        """Apply a controlled-Z rotation gate.
+
+        See :class:`~graphix.instruction.CRZ` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(Instruction.CRZ(control=control, target=target, angle=angle))
+
+    def cr(self, control: int, target: int, axis: Axis, angle: ParameterizedAngle) -> None:
+        """Apply a controlled-rotation gate on the given axis.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        axis : Axis
+            rotation axis
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        match axis:
+            case Axis.X:
+                self.crx(control, target, angle)
+            case Axis.Y:
+                self.cry(control, target, angle)
+            case Axis.Z:
+                self.crz(control, target, angle)
+            case _:
+                assert_never(axis)
+
+    def cu(
+        self,
+        control: int,
+        target: int,
+        theta: ParameterizedAngle,
+        phi: ParameterizedAngle,
+        lambda_: ParameterizedAngle,
+        gamma: ParameterizedAngle,
+    ) -> None:
+        """Apply a controlled-U gate.
+
+        See :class:`~graphix.instruction.CU` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        target : int
+            target qubit
+        theta : ParameterizedAngle
+            rotation angle in units of π
+        phi : ParameterizedAngle
+            rotation angle in units of π
+        lambda_ : ParameterizedAngle
+            rotation angle in units of π
+        gamma : ParameterizedAngle
+            rotation angle in units of π
+        """
+        assert control in self.active_qubits
+        assert target in self.active_qubits
+        assert control != target
+        self.instruction.append(
+            Instruction.CU(control=control, target=target, theta=theta, phi=phi, lambda_=lambda_, gamma=gamma)
+        )
+
+    def cswap(self, control: int, qubit1: int, qubit2: int) -> None:
+        """Apply a CSWAP gate.
+
+        See :class:`~graphix.instruction.CSWAP` for more information.
+
+        Parameters
+        ----------
+        control : int
+            control qubit
+        qubit1 : int
+            first qubit to be swapped
+        qubit2 : int
+            second qubit to be swapped
+        """
+        assert control in self.active_qubits
+        assert qubit1 in self.active_qubits
+        assert qubit2 in self.active_qubits
+        assert control != qubit1
+        assert control != qubit2
+        assert qubit1 != qubit2
+        self.instruction.append(Instruction.CSWAP(control=control, targets=(qubit1, qubit2)))
+
+    def gphase(self, angle: ParameterizedAngle) -> None:
+        r"""Apply a global phase.
+
+        See :class:`~graphix.instruction.GPHASE` for more information.
+
+        Parameters
+        ----------
+        angle : ParameterizedAngle
+            rotation angle in units of π
+        """
+        self.instruction.append(Instruction.GPHASE(angle))
+
     def transpile_to_causalflow(self) -> TranspiledFlow:
         """Transpile a circuit via J-∧z decomposition to a causal flow.
 
@@ -476,6 +812,9 @@ def transpile_to_causalflow(self) -> TranspiledFlow:
                     else:
                         graph.add_edge(i0, i1)
                     continue
+                case InstructionKind.GPHASE:
+                    # Global phase is currently ignored
+                    pass
                 case _:
                     assert_never(instr.kind)
         outputs = [i for i in indices if i is not None]
@@ -590,6 +929,8 @@ def simulate(
 
         classical_measures: list[Outcome] = []
 
+        gphase: ParameterizedAngle = 0
+
         for i in range(len(self.instruction)):
             instr = self.instruction[i]
 
@@ -605,6 +946,8 @@ def evolve(op: Matrix, qargs: Iterable[int]) -> None:
                 case instruction.InstructionKind.SWAP:
                     u, v = instr.targets
                     _backend.state.swap((_backend.node_index.index(u), _backend.node_index.index(v)))
+                case instruction.InstructionKind.CY:
+                    evolve(Ops.CY, [instr.control, instr.target])
                 case instruction.InstructionKind.CZ:
                     u, v = instr.targets
                     _backend.state.entangle((_backend.node_index.index(u), _backend.node_index.index(v)))
@@ -612,6 +955,16 @@ def evolve(op: Matrix, qargs: Iterable[int]) -> None:
                     pass
                 case instruction.InstructionKind.S:
                     evolve_single(Ops.S, instr.target)
+                case instruction.InstructionKind.SDG:
+                    evolve_single(Ops.SDG, instr.target)
+                case instruction.InstructionKind.T:
+                    evolve_single(Ops.T, instr.target)
+                case instruction.InstructionKind.TDG:
+                    evolve_single(Ops.TDG, instr.target)
+                case instruction.InstructionKind.SX:
+                    evolve_single(Ops.SX, instr.target)
+                case instruction.InstructionKind.SXDG:
+                    evolve_single(Ops.SXDG, instr.target)
                 case instruction.InstructionKind.H:
                     evolve_single(Ops.H, instr.target)
                 case instruction.InstructionKind.X:
@@ -620,6 +973,8 @@ def evolve(op: Matrix, qargs: Iterable[int]) -> None:
                     evolve_single(Ops.Y, instr.target)
                 case instruction.InstructionKind.Z:
                     evolve_single(Ops.Z, instr.target)
+                case instruction.InstructionKind.P:
+                    evolve_single(Ops.p(instr.angle), instr.target)
                 case instruction.InstructionKind.RX:
                     evolve_single(Ops.rx(instr.angle), instr.target)
                 case instruction.InstructionKind.RY:
@@ -628,17 +983,36 @@ def evolve(op: Matrix, qargs: Iterable[int]) -> None:
                     evolve_single(Ops.rz(instr.angle), instr.target)
                 case instruction.InstructionKind.J:
                     evolve_single(Ops.j(instr.angle), instr.target)
+                case instruction.InstructionKind.CJ:
+                    evolve(Ops.cj(instr.angle), [instr.control, instr.target])
+                case instruction.InstructionKind.U:
+                    evolve_single(Ops.u(instr.theta, instr.phi, instr.lambda_), instr.target)
+                case instruction.InstructionKind.CU:
+                    evolve(Ops.cu(instr.theta, instr.phi, instr.lambda_, instr.gamma), [instr.control, instr.target])
+                case instruction.InstructionKind.CP:
+                    evolve(Ops.cp(instr.angle), [instr.control, instr.target])
+                case instruction.InstructionKind.CRX:
+                    evolve(Ops.crx(instr.angle), [instr.control, instr.target])
+                case instruction.InstructionKind.CRY:
+                    evolve(Ops.cry(instr.angle), [instr.control, instr.target])
+                case instruction.InstructionKind.CRZ:
+                    evolve(Ops.crz(instr.angle), [instr.control, instr.target])
                 case instruction.InstructionKind.RZZ:
                     evolve(Ops.rzz(instr.angle), [instr.control, instr.target])
                 case instruction.InstructionKind.CCX:
                     evolve(Ops.CCX, [instr.controls[0], instr.controls[1], instr.target])
+                case instruction.InstructionKind.CSWAP:
+                    evolve(Ops.CSWAP, [instr.control, instr.targets[0], instr.targets[1]])
                 case instruction.InstructionKind.M:
                     result = _backend.measure(
                         instr.target, PauliMeasurement(instr.axis), rng=rng, stacklevel=stacklevel + 1
                     )
                     classical_measures.append(result)
+                case InstructionKind.GPHASE:
+                    gphase += instr.angle
                 case _:
-                    raise ValueError(f"Unknown instruction: {instr}")
+                    assert_never(instr.kind)
+        # Global phase is currently ignored
         return SimulateResult(_backend.state, tuple(classical_measures))
 
     def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Circuit:
@@ -749,6 +1123,28 @@ def transpile_j_to_rzh(self) -> Circuit:
                     new_circuit.add(instr)
         return new_circuit
 
+    def transpile_cj(self) -> Circuit:
+        """Return an equivalent circuit where all CJ gates have been replaced with OpenQASM gates."""
+        new_circuit = Circuit(self.width)
+        for instr in self.instruction:
+            match instr.kind:
+                case InstructionKind.CJ:
+                    new_circuit.extend(decompose_cj(instr))
+                case _:
+                    new_circuit.add(instr)
+        return new_circuit
+
+    def transpile_rzz(self) -> Circuit:
+        """Return an equivalent circuit where all RZZ gates have been replaced with OpenQASM gates."""
+        new_circuit = Circuit(self.width)
+        for instr in self.instruction:
+            match instr.kind:
+                case InstructionKind.RZZ:
+                    new_circuit.extend(decompose_rzz(instr))
+                case _:
+                    new_circuit.add(instr)
+        return new_circuit
+
 
 def decompose_rzz(instr: instruction.RZZ) -> Iterator[instruction.CNOT | instruction.RZ]:
     """Yield a decomposition of RZZ(α) gate as CNOT(control, target)·Rz(target, α)·CNOT(control, target).
@@ -847,8 +1243,8 @@ def decompose_swap(instr: instruction.SWAP) -> Iterator[instruction.CNOT]:
     yield instruction.CNOT(control=instr.targets[0], target=instr.targets[1])
 
 
-def decompose_y(instr: instruction.Y) -> Iterator[instruction.X | instruction.Z]:
-    """Return a decomposition of the Y gate as X·Z.
+def decompose_y(instr: instruction.Y) -> Iterator[instruction.X | instruction.Z | Instruction.GPHASE]:
+    r"""Return a decomposition of the Y gate as :math:`\mathrm e^{\mathrm i \frac \pi 2} X Z`.
 
     Parameters
     ----------
@@ -861,9 +1257,10 @@ def decompose_y(instr: instruction.Y) -> Iterator[instruction.X | instruction.Z]
     """
     yield instruction.Z(instr.target)
     yield instruction.X(instr.target)
+    yield instruction.GPHASE(ANGLE_PI / 2)
 
 
-def decompose_rx(instr: instruction.RX) -> Iterator[instruction.J]:
+def decompose_rx(instr: instruction.RX) -> Iterator[instruction.J | Instruction.GPHASE]:
     """Yield a J decomposition of the RX gate.
 
     The Rx(α) gate is decomposed into J(α)·H (that is to say, J(α)·J(0)).
@@ -880,9 +1277,10 @@ def decompose_rx(instr: instruction.RX) -> Iterator[instruction.J]:
     """
     yield instruction.J(instr.target, 0)
     yield instruction.J(instr.target, instr.angle)
+    yield instruction.GPHASE(-instr.angle / 2)
 
 
-def decompose_ry(instr: instruction.RY) -> Iterator[instruction.J]:
+def decompose_ry(instr: instruction.RY) -> Iterator[instruction.J | Instruction.GPHASE]:
     """Yield a J decomposition of the RY gate.
 
     The Ry(α) gate is decomposed into J(0)·J(π/2)·J(α)·J(-π/2).
@@ -902,9 +1300,10 @@ def decompose_ry(instr: instruction.RY) -> Iterator[instruction.J]:
     yield instruction.J(target=instr.target, angle=instr.angle)
     yield instruction.J(target=instr.target, angle=ANGLE_PI / 2)
     yield instruction.J(target=instr.target, angle=0)
+    yield instruction.GPHASE(-instr.angle / 2)
 
 
-def decompose_rz(instr: instruction.RZ) -> Iterator[instruction.J]:
+def decompose_rz(instr: instruction.RZ) -> Iterator[instruction.J | Instruction.GPHASE]:
     """Yield a J decomposition of the RZ gate.
 
     The Rz(α) gate is decomposed into H·J(α) (that is to say, J(0)·J(α)).
@@ -921,9 +1320,129 @@ def decompose_rz(instr: instruction.RZ) -> Iterator[instruction.J]:
     """
     yield instruction.J(target=instr.target, angle=instr.angle)
     yield instruction.J(target=instr.target, angle=0)
+    yield instruction.GPHASE(-instr.angle / 2)
+
+
+def decompose_u(instr: instruction.U) -> Iterator[instruction.J | Instruction.GPHASE]:
+    """Yield a J decomposition of the U gate.
+
+    The U(θ, φ, λ) gate is decomposed into H·J(φ + 𝜋/2)·J(θ)·J(λ - 𝜋/2) (that is to say, J(0)·J(φ + 𝜋/2)·J(θ)·J(λ - 𝜋/2)).
+
+    Parameters
+    ----------
+        instr: the U instruction to decompose.
+
+    Returns
+    -------
+        the decomposition.
+
+    """
+    yield Instruction.J(instr.target, instr.lambda_ - ANGLE_PI / 2)
+    yield Instruction.J(instr.target, instr.theta)
+    yield Instruction.J(instr.target, instr.phi + ANGLE_PI / 2)
+    yield Instruction.J(instr.target, 0)
+    yield Instruction.GPHASE(-instr.theta / 2)
+
+
+def decompose_cu(instr: instruction.CU) -> Iterator[instruction.CJ | Instruction.P]:
+    """Yield a J decomposition of the U gate.
+
+    The U(θ, φ, λ) gate is decomposed into H·J(φ + 𝜋/2)·J(θ)·J(λ - 𝜋/2) (that is to say, J(0)·J(φ + 𝜋/2)·J(θ)·J(λ - 𝜋/2)).
+
+    Parameters
+    ----------
+        instr: the U instruction to decompose.
+
+    Returns
+    -------
+        the decomposition.
+
+    """
+    yield Instruction.CJ(control=instr.control, target=instr.target, angle=instr.lambda_ - ANGLE_PI / 2)
+    yield Instruction.CJ(control=instr.control, target=instr.target, angle=instr.theta)
+    yield Instruction.CJ(control=instr.control, target=instr.target, angle=instr.phi + ANGLE_PI / 2)
+    yield Instruction.CJ(control=instr.control, target=instr.target, angle=0)
+    yield Instruction.P(target=instr.control, angle=instr.gamma - instr.theta / 2)
+
+
+def insert_control(
+    control: int,
+    instrs: Iterable[
+        Instruction.GPHASE
+        | Instruction.X
+        | Instruction.Z
+        | Instruction.J
+        | Instruction.CZ
+        | Instruction.CNOT
+        | Instruction.RZ
+    ],
+) -> Iterable[InstructionType]:
+    """Yield a controlled gate sequence from a gate sequence.
+
+    Parameters
+    ----------
+    control: int
+        The control qubit.
+    instrs: Iterable[Instruction.J | Instruction.CZ]
+        The J-∧z decomposition.
+
+    Yields
+    ------
+    InstructionType
+        The controlled gate sequence.
+    """
+    gphase: ParameterizedAngle = 0
+    for instr in instrs:
+        match instr.kind:
+            case InstructionKind.X:
+                yield instruction.CNOT(control=control, target=instr.target)
+            case InstructionKind.Z:
+                yield instruction.CZ((control, instr.target))
+            case InstructionKind.J:
+                yield instruction.CJ(control=control, target=instr.target, angle=instr.angle)
+            case InstructionKind.CZ:
+                u, v = instr.targets
+                yield instruction.H(v)
+                yield instruction.CCX(target=v, controls=(control, u))
+                yield instruction.H(v)
+            case InstructionKind.CNOT:
+                yield instruction.CCX(target=instr.target, controls=(control, instr.control))
+            case InstructionKind.RZ:
+                yield instruction.CRZ(control=control, target=instr.target, angle=instr.angle)
+            case InstructionKind.GPHASE:
+                gphase += instr.angle
+            case _:
+                assert_never(instr.kind)
+    yield Instruction.P(target=control, angle=gphase)
+
+
+def decompose_cj(instr: Instruction.CJ) -> Iterator[InstructionType]:
+    """Yield a decomposed gate sequence of the CJ gate.
+
+    See :class:`~graphix.instruction.CJ` for more information.
+    """
+    delta = (instr.angle + ANGLE_PI) / 2
+    yield instruction.RZ(target=instr.target, angle=delta)
+    yield instruction.CNOT(control=instr.control, target=instr.target)
+    yield instruction.RZ(target=instr.target, angle=-delta)
+    yield instruction.RY(target=instr.target, angle=-ANGLE_PI / 4)
+    yield instruction.CNOT(control=instr.control, target=instr.target)
+    yield instruction.RY(target=instr.target, angle=ANGLE_PI / 4)
+    yield instruction.P(target=instr.control, angle=delta)
+
+
+def decompose_p(instr: Instruction.P) -> Iterator[Instruction.RZ | Instruction.GPHASE]:
+    """Yield a decomposed gate sequence of the P gate.
+
+    See :class:`~graphix.instruction.P` for more information.
+    """
+    yield Instruction.RZ(instr.target, instr.angle)
+    yield Instruction.GPHASE(instr.angle / 2)
 
 
-def instructions_to_jcz(instrs: Iterable[InstructionType]) -> Iterator[instruction.J | instruction.CZ | instruction.M]:
+def instructions_to_jcz(
+    instrs: Iterable[InstructionType],
+) -> Iterator[instruction.J | instruction.CZ | instruction.M | Instruction.GPHASE]:
     """Yield a J-∧z decomposition of the instruction.
 
     Parameters
@@ -945,6 +1464,16 @@ def instructions_to_jcz(instrs: Iterable[InstructionType]) -> Iterator[instructi
                 yield instruction.J(instr.target, 0)
             case InstructionKind.S:
                 yield from decompose_rz(instruction.RZ(instr.target, ANGLE_PI / 2))
+            case InstructionKind.SDG:
+                yield from decompose_rz(instruction.RZ(instr.target, -ANGLE_PI / 2))
+            case InstructionKind.T:
+                yield from decompose_rz(instruction.RZ(instr.target, ANGLE_PI / 4))
+            case InstructionKind.TDG:
+                yield from decompose_rz(instruction.RZ(instr.target, -ANGLE_PI / 4))
+            case InstructionKind.SX:
+                yield from decompose_rx(instruction.RX(instr.target, ANGLE_PI / 2))
+            case InstructionKind.SXDG:
+                yield from decompose_rx(instruction.RX(instr.target, -ANGLE_PI / 2))
             case InstructionKind.X:
                 yield from decompose_rx(instruction.RX(instr.target, ANGLE_PI))
             case InstructionKind.Y:
@@ -957,6 +1486,10 @@ def instructions_to_jcz(instrs: Iterable[InstructionType]) -> Iterator[instructi
                 yield from decompose_ry(instr)
             case InstructionKind.RZ:
                 yield from decompose_rz(instr)
+            case InstructionKind.P:
+                yield from instructions_to_jcz(decompose_p(instr))
+            case InstructionKind.U:
+                yield from decompose_u(instr)
             case InstructionKind.CCX:
                 yield from instructions_to_jcz(decompose_ccx(instr))
             case InstructionKind.RZZ:
@@ -965,6 +1498,35 @@ def instructions_to_jcz(instrs: Iterable[InstructionType]) -> Iterator[instructi
                 yield from instructions_to_jcz(decompose_cnot(instr))
             case InstructionKind.SWAP:
                 yield from instructions_to_jcz(decompose_swap(instr))
+            case InstructionKind.CY:
+                yield from instructions_to_jcz(insert_control(instr.control, decompose_y(Instruction.Y(instr.target))))
+            case InstructionKind.CJ:
+                yield from instructions_to_jcz(decompose_cj(instr))
+            case InstructionKind.CP:
+                yield from instructions_to_jcz(
+                    insert_control(instr.control, decompose_p(Instruction.P(instr.target, instr.angle)))
+                )
+            case InstructionKind.CRX:
+                yield from instructions_to_jcz(
+                    insert_control(instr.control, decompose_rx(Instruction.RX(instr.target, instr.angle)))
+                )
+            case InstructionKind.CRY:
+                yield from instructions_to_jcz(
+                    insert_control(instr.control, decompose_ry(Instruction.RY(instr.target, instr.angle)))
+                )
+            case InstructionKind.CRZ:
+                yield from instructions_to_jcz(
+                    insert_control(instr.control, decompose_rz(Instruction.RZ(instr.target, instr.angle)))
+                )
+            case InstructionKind.CU:
+                yield from instructions_to_jcz(decompose_cu(instr))
+            case InstructionKind.CSWAP:
+                yield from instructions_to_jcz(
+                    insert_control(instr.control, decompose_swap(Instruction.SWAP(instr.targets)))
+                )
+            case InstructionKind.GPHASE:
+                # Global phase is currently ignored
+                pass
             case _:
                 assert_never(instr.kind)
 
diff --git a/tests/test_instruction.py b/tests/test_instruction.py
index 4228ece67..c2c520a2c 100644
--- a/tests/test_instruction.py
+++ b/tests/test_instruction.py
@@ -1,38 +1,82 @@
 from __future__ import annotations
 
 from copy import copy
+from dataclasses import dataclass
 from typing import TYPE_CHECKING
 
+import numpy as np
 import pytest
 
 # override introduced in Python 3.12
 from typing_extensions import override
 
-from graphix import ANGLE_PI, Axis, Clifford
-from graphix.instruction import Instruction, InstructionVisitor
+from graphix import ANGLE_PI, Axis, Clifford, Instruction
+from graphix.fundamentals import angle_to_rad
+from graphix.instruction import InstructionVisitor
+from graphix.ops import Ops
 
 if TYPE_CHECKING:
+    from collections.abc import Callable
+
+    from numpy.random import Generator
+
     from graphix.fundamentals import ParameterizedAngle
     from graphix.instruction import InstructionType
 
-ALL_INSTRUCTIONS = [
-    Instruction.CCX(target=0, controls=(1, 2)),
-    Instruction.RZZ(target=0, control=1, angle=ANGLE_PI / 4),
-    Instruction.CNOT(target=0, control=1),
-    Instruction.SWAP(targets=(0, 1)),
-    Instruction.CZ(targets=(0, 1)),
-    Instruction.H(target=0),
-    Instruction.S(target=0),
-    Instruction.X(target=0),
-    Instruction.Y(target=0),
-    Instruction.Z(target=0),
-    Instruction.I(target=0),
-    Instruction.RX(target=0, angle=ANGLE_PI / 4),
-    Instruction.RY(target=0, angle=ANGLE_PI / 4),
-    Instruction.RZ(target=0, angle=ANGLE_PI / 4),
-    Instruction.J(target=0, angle=ANGLE_PI / 4),
-    Instruction.M(target=0, axis=Axis.X),
-]
+
+@dataclass(frozen=True)
+class InstructionTestCase:
+    name: str
+    instruction: Callable[[Generator], InstructionType]
+
+
+INSTRUCTION_TEST_CASES: tuple[InstructionTestCase, ...] = (
+    InstructionTestCase("CCX", lambda _rng: Instruction.CCX(0, (1, 2))),
+    InstructionTestCase("RZZ", lambda rng: Instruction.RZZ(0, 1, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("CZ", lambda _rng: Instruction.CZ((0, 1))),
+    InstructionTestCase("CNOT", lambda _rng: Instruction.CNOT(0, 1)),
+    InstructionTestCase("SWAP", lambda _rng: Instruction.SWAP((0, 1))),
+    InstructionTestCase("H", lambda _rng: Instruction.H(0)),
+    InstructionTestCase("S", lambda _rng: Instruction.S(0)),
+    InstructionTestCase("SDG", lambda _rng: Instruction.SDG(0)),
+    InstructionTestCase("T", lambda _rng: Instruction.T(0)),
+    InstructionTestCase("TDG", lambda _rng: Instruction.TDG(0)),
+    InstructionTestCase("SX", lambda _rng: Instruction.SX(0)),
+    InstructionTestCase("SXDG", lambda _rng: Instruction.SXDG(0)),
+    InstructionTestCase("X", lambda _rng: Instruction.X(0)),
+    InstructionTestCase("Y", lambda _rng: Instruction.Y(0)),
+    InstructionTestCase("Z", lambda _rng: Instruction.Z(0)),
+    InstructionTestCase("I", lambda _rng: Instruction.I(0)),
+    InstructionTestCase("RX", lambda rng: Instruction.RX(0, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("RY", lambda rng: Instruction.RY(0, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("RZ", lambda rng: Instruction.RZ(0, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("J", lambda rng: Instruction.J(0, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("P", lambda rng: Instruction.P(0, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase(
+        "U",
+        lambda rng: Instruction.U(
+            0, rng.random() * 2 * ANGLE_PI, rng.random() * 2 * ANGLE_PI, rng.random() * 2 * ANGLE_PI
+        ),
+    ),
+    InstructionTestCase("CY", lambda _rng: Instruction.CY(0, 1)),
+    InstructionTestCase("CJ", lambda rng: Instruction.CJ(0, 1, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("CP", lambda rng: Instruction.CP(0, 1, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("CRX", lambda rng: Instruction.CRX(0, 1, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("CRY", lambda rng: Instruction.CRY(0, 1, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase("CRZ", lambda rng: Instruction.CRZ(0, 1, rng.random() * 2 * ANGLE_PI)),
+    InstructionTestCase(
+        "CU",
+        lambda rng: Instruction.CU(
+            0,
+            1,
+            rng.random() * 2 * ANGLE_PI,
+            rng.random() * 2 * ANGLE_PI,
+            rng.random() * 2 * ANGLE_PI,
+            rng.random() * 2 * ANGLE_PI,
+        ),
+    ),
+    InstructionTestCase("CSWAP", lambda _rng: Instruction.CSWAP(0, (1, 2))),
+)
 
 
 class VisitQubit(InstructionVisitor):
@@ -53,43 +97,85 @@ def visit_axis(self, axis: Axis) -> Axis:
         return axis.clifford(Clifford.H)
 
 
-@pytest.mark.parametrize("instruction", ALL_INSTRUCTIONS)
-def test_visit_qubit(instruction: InstructionType) -> None:
-    # Copy the instruction to keep ALL_INSTRUCTIONS unmodified
-    instr_copy = copy(instruction)
+@pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+def test_visit_qubit(fx_rng: Generator, test_case: InstructionTestCase) -> None:
+    instr = test_case.instruction(fx_rng)
+    instr_copy = copy(instr)
     visitor = VisitQubit()
-    instr_visited = instr_copy.visit(visitor, copy=True)
-    assert instr_copy == instruction
-    assert instr_visited != instruction
-    instr_copy.visit(visitor, copy=False)
-    assert instr_copy != instruction
-    assert instr_visited == instr_copy
-
-
-@pytest.mark.parametrize("instruction", ALL_INSTRUCTIONS)
-def test_visit_angle(instruction: InstructionType) -> None:
-    if not hasattr(instruction, "angle"):
+    instr_visited = instr.visit(visitor, copy=True)
+    assert instr == instr_copy
+    assert instr_visited != instr_copy
+    instr.visit(visitor, copy=False)
+    assert instr != instr_copy
+    assert instr_visited == instr
+
+
+@pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+def test_visit_angle(fx_rng: Generator, test_case: InstructionTestCase) -> None:
+    instr = test_case.instruction(fx_rng)
+    if not hasattr(instr, "angle"):
         pytest.skip()
-    # Copy the instruction to keep ALL_INSTRUCTIONS unmodified
-    instr_copy = copy(instruction)
+    instr_copy = copy(instr)
     visitor = VisitAngle()
-    instr_visited = instr_copy.visit(visitor, copy=True)
-    assert instr_copy == instruction
-    assert instr_visited != instruction
-    instr_copy.visit(visitor, copy=False)
-    assert instr_copy != instruction
-    assert instr_visited == instr_copy
-
-
-@pytest.mark.parametrize("instruction", ALL_INSTRUCTIONS)
-def test_visit_axis(instruction: InstructionType) -> None:
-    if not hasattr(instruction, "axis"):
+    instr_visited = instr.visit(visitor, copy=True)
+    assert instr == instr_copy
+    assert instr_visited != instr_copy
+    instr.visit(visitor, copy=False)
+    assert instr != instr_copy
+    assert instr_visited == instr
+
+
+@pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+def test_visit_axis(fx_rng: Generator, test_case: InstructionTestCase) -> None:
+    instr = test_case.instruction(fx_rng)
+    if not hasattr(instr, "axis"):
         pytest.skip()
-    instr_copy = copy(instruction)
+    instr_copy = copy(instr)
     visitor = VisitAxis()
-    instr_visited = instr_copy.visit(visitor, copy=True)
-    assert instr_copy == instruction
-    assert instr_visited != instruction
-    instr_copy.visit(visitor, copy=False)
-    assert instr_copy != instruction
-    assert instr_visited == instr_copy
+    instr_visited = instr.visit(visitor, copy=True)
+    assert instr == instr_copy
+    assert instr_visited != instr_copy
+    instr.visit(visitor, copy=False)
+    assert instr != instr_copy
+    assert instr_visited == instr
+
+
+def test_u(fx_rng: Generator) -> None:
+    theta = fx_rng.random()
+    phi = fx_rng.random()
+    lambda_ = fx_rng.random()
+    np.testing.assert_allclose(
+        Ops.u(theta, phi, lambda_),
+        np.exp(-1j * angle_to_rad(theta) / 2)
+        * (Ops.H @ Ops.j(phi + ANGLE_PI / 2) @ Ops.j(theta) @ Ops.j(lambda_ - ANGLE_PI / 2)),
+    )
+
+
+def test_cj(fx_rng: Generator) -> None:
+    alpha = fx_rng.random()
+    delta = (alpha + ANGLE_PI) / 2
+    a = Ops.ry(ANGLE_PI / 4)
+    b = Ops.ry(-ANGLE_PI / 4) @ Ops.rz(-delta)
+    c = Ops.rz(delta)
+    np.testing.assert_allclose(a @ b @ c, Ops.I, atol=1e-15)
+    np.testing.assert_allclose(a @ Ops.X @ b @ Ops.X @ c, np.exp(-1j * angle_to_rad(delta)) * Ops.j(alpha))
+    np.testing.assert_allclose(
+        Ops.cj(alpha),
+        np.kron(Ops.p(delta), Ops.I) @ np.kron(Ops.I, a) @ Ops.CNOT @ np.kron(Ops.I, b) @ Ops.CNOT @ np.kron(Ops.I, c),
+        atol=1e-15,
+    )
+
+
+def test_cu(fx_rng: Generator) -> None:
+    theta = fx_rng.random()
+    phi = fx_rng.random()
+    lambda_ = fx_rng.random()
+    gamma = fx_rng.random()
+    np.testing.assert_allclose(
+        Ops.cu(theta, phi, lambda_, gamma),
+        np.kron(Ops.p(gamma - theta / 2), Ops.I)
+        @ Ops.cj(0)
+        @ Ops.cj(phi + ANGLE_PI / 2)
+        @ Ops.cj(theta)
+        @ Ops.cj(lambda_ - ANGLE_PI / 2),
+    )
diff --git a/tests/test_qasm3_exporter_to_graphix_parser.py b/tests/test_qasm3_exporter_to_graphix_parser.py
index 87dcd7d16..8133c18ff 100644
--- a/tests/test_qasm3_exporter_to_graphix_parser.py
+++ b/tests/test_qasm3_exporter_to_graphix_parser.py
@@ -2,6 +2,8 @@
 
 from __future__ import annotations
 
+import dataclasses
+import math
 from typing import TYPE_CHECKING
 
 import pytest
@@ -12,10 +14,10 @@
 from graphix.instruction import InstructionKind
 from graphix.qasm3_exporter import circuit_to_qasm3
 from graphix.random_objects import rand_circuit
-from tests.test_instruction import ALL_INSTRUCTIONS
+from tests.test_instruction import INSTRUCTION_TEST_CASES
 
 if TYPE_CHECKING:
-    from graphix.instruction import InstructionType
+    from tests.test_instruction import InstructionTestCase
 
 try:
     from graphix_qasm_parser import OpenQASMParser  # type: ignore[import-not-found, unused-ignore]
@@ -36,7 +38,13 @@ def check_round_trip(circuit: Circuit) -> None:
     check_circuit = circuit.transpile_j_to_rzh()
     parser = OpenQASMParser()
     parsed_circuit = parser.parse_str(qasm)
-    assert parsed_circuit.instruction == check_circuit.instruction
+    for parsed_instr, instr in zip(parsed_circuit.instruction, check_circuit.instruction, strict=True):
+        assert parsed_instr.kind == instr.kind
+        assert all(
+            math.isclose(x, y) if isinstance(x, float) and isinstance(y, float) else x == y
+            for field in dataclasses.fields(parsed_instr)
+            for x, y in [(getattr(parsed_instr, field.name), getattr(instr, field.name))]
+        )
 
 
 @pytest.mark.parametrize("jumps", range(1, 11))
@@ -48,11 +56,12 @@ def test_circuit_to_qasm3(fx_bg: PCG64, jumps: int) -> None:
     check_round_trip(rand_circuit(nqubits, depth, rng, use_j=True, use_cz=True))
 
 
-@pytest.mark.parametrize("instruction", ALL_INSTRUCTIONS)
-def test_instruction_to_qasm3(instruction: InstructionType) -> None:
-    if instruction.kind == InstructionKind.M:
+@pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+def test_instruction_to_qasm3(fx_rng: Generator, test_case: InstructionTestCase) -> None:
+    instr = test_case.instruction(fx_rng)
+    if instr.kind in {InstructionKind.CJ, InstructionKind.RZZ, InstructionKind.M}:
         pytest.skip()
-    check_round_trip(Circuit(3, instr=[instruction]))
+    check_round_trip(Circuit(3, instr=[instr]))
 
 
 def test_j_to_qasm3() -> None:
diff --git a/tests/test_qasm3_exporter_to_qiskit.py b/tests/test_qasm3_exporter_to_qiskit.py
index 9b0ff7028..c8f107dde 100644
--- a/tests/test_qasm3_exporter_to_qiskit.py
+++ b/tests/test_qasm3_exporter_to_qiskit.py
@@ -13,16 +13,19 @@
 from graphix.clifford import Clifford
 from graphix.command import C, CommandKind, E, M, N
 from graphix.fundamentals import Plane
+from graphix.instruction import InstructionKind
 from graphix.measurements import BlochMeasurement, Measurement, outcome
 from graphix.optimization import single_qubit_domains
-from graphix.qasm3_exporter import pattern_to_qasm3
+from graphix.qasm3_exporter import circuit_to_qasm3, pattern_to_qasm3
 from graphix.random_objects import rand_circuit
 from graphix.sim.statevec import StatevectorBackend
 from graphix.states import BasicStates
+from tests.test_instruction import INSTRUCTION_TEST_CASES
 
 if TYPE_CHECKING:
     from graphix.measurements import Outcome
     from graphix.states import State
+    from tests.test_instruction import InstructionTestCase
 
 try:
     import qiskit
@@ -40,7 +43,7 @@
         sys.exit(1)
 
 
-def check_qasm3(pattern: Pattern) -> None:
+def check_qasm3_pattern(pattern: Pattern) -> None:
     """Check that we obtain equivalent statevectors whether we simulate the pattern with Graphix or we use Qiskit AER simulator."""
     qasm3 = pattern_to_qasm3(pattern)
     qc = qiskit_qasm3_import.parse(qasm3)
@@ -81,13 +84,13 @@ def check_qasm3(pattern: Pattern) -> None:
 
 
 def test_to_qasm3_qubits_preparation() -> None:
-    check_qasm3(Pattern(cmds=[N(0), N(1)]))
-    check_qasm3(Pattern(input_nodes=[0], cmds=[N(1)]))
+    check_qasm3_pattern(Pattern(cmds=[N(0), N(1)]))
+    check_qasm3_pattern(Pattern(input_nodes=[0], cmds=[N(1)]))
 
 
 def test_to_qasm3_entanglement() -> None:
-    check_qasm3(Pattern(input_nodes=[0, 1], cmds=[E((0, 1))]))
-    check_qasm3(Pattern(input_nodes=[0, 1], cmds=[N(2), E((1, 2))]))
+    check_qasm3_pattern(Pattern(input_nodes=[0, 1], cmds=[E((0, 1))]))
+    check_qasm3_pattern(Pattern(input_nodes=[0, 1], cmds=[N(2), E((1, 2))]))
 
 
 @pytest.mark.parametrize("clifford", Clifford)
@@ -95,21 +98,21 @@ def test_to_qasm3_entanglement() -> None:
     "state", [BasicStates.ZERO, BasicStates.PLUS, pytest.param(BasicStates.MINUS, marks=pytest.mark.xfail)]
 )
 def test_to_qasm3_clifford(clifford: Clifford, state: State) -> None:
-    check_qasm3(Pattern(cmds=[N(0, state), C(0, clifford)]))
+    check_qasm3_pattern(Pattern(cmds=[N(0, state), C(0, clifford)]))
 
 
 @pytest.mark.parametrize("state", [BasicStates.ZERO, BasicStates.PLUS])
 @pytest.mark.parametrize("plane", list(Plane))
 @pytest.mark.parametrize("angle", [0, 0.25, 1.75])
 def test_to_qasm3_measurement(state: State, plane: Plane, angle: float) -> None:
-    check_qasm3(Pattern(cmds=[N(0, state), N(1), E((0, 1)), M(0, BlochMeasurement(angle, plane))]))
+    check_qasm3_pattern(Pattern(cmds=[N(0, state), N(1), E((0, 1)), M(0, BlochMeasurement(angle, plane))]))
 
 
 def test_to_qasm3_hadamard() -> None:
     circuit = Circuit(1)
     circuit.h(0)
     pattern = circuit.transpile().pattern
-    check_qasm3(pattern)
+    check_qasm3_pattern(pattern)
 
 
 @pytest.mark.parametrize("jumps", range(1, 11))
@@ -126,4 +129,29 @@ def test_to_qasm3_random_circuit(fx_bg: PCG64, jumps: int) -> None:
     # qiskit_qasm3_import.exceptions.ConversionError: unhandled binary operator '^'
     pattern = single_qubit_domains(pattern)
 
-    check_qasm3(pattern)
+    check_qasm3_pattern(pattern)
+
+
+def check_qasm3_circuit(circuit: Circuit) -> None:
+    """Check that we obtain equivalent statevectors whether we simulate the circuit with Graphix or we use Qiskit AER simulator."""
+    qasm3 = circuit_to_qasm3(circuit)
+    qc = qiskit_qasm3_import.parse(qasm3)
+    qc.save_statevector()  # type:ignore[attr-defined]
+    aer_backend = AerSimulator(method="statevector")
+    transpiled = qiskit.transpile(qc, aer_backend)
+    result = aer_backend.run(transpiled, shots=1, memory=True).result()
+    state_qiskit = result.get_statevector()
+    n = int(np.log2(len(state_qiskit)))
+    state_qiskit = state_qiskit.reshape((2,) * n).transpose(*reversed(range(n))).reshape(-1)
+    state_graphix = circuit.simulate(input_state=BasicStates.ZERO).state
+    assert state_graphix.isclose(state_qiskit)
+
+
+@pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+def test_instruction_to_qasm3(fx_rng: Generator, test_case: InstructionTestCase) -> None:
+    instr = test_case.instruction(fx_rng)
+    if instr.kind in {InstructionKind.CJ, InstructionKind.RZZ, InstructionKind.M}:
+        pytest.skip()
+    if instr.kind in {InstructionKind.SXDG, InstructionKind.U}:
+        pytest.skip("qiskit_qasm3_import.exceptions.ConversionError: gate 'sxdg'/'u' is not defined.")
+    check_qasm3_circuit(Circuit(3, instr=[instr]))
diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py
index 0624ce506..f41ff6dc4 100644
--- a/tests/test_transpiler.py
+++ b/tests/test_transpiler.py
@@ -7,63 +7,55 @@
 import pytest
 from numpy.random import PCG64, Generator
 
-from graphix import instruction
+from graphix import Instruction, instruction
 from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
-from graphix.fundamentals import ANGLE_PI, Axis, Sign
+from graphix.fundamentals import Axis, Sign
 from graphix.instruction import I, InstructionKind
 from graphix.random_objects import rand_circuit, rand_gate, rand_state_vector
 from graphix.sim.density_matrix import DensityMatrix
 from graphix.sim.statevec import Statevector, StatevectorBackend
 from graphix.simulator import DefaultMeasureMethod
 from graphix.states import BasicStates
-from graphix.transpiler import Circuit, OutputIndex, OutputKind, decompose_ccx, transpile_swaps
+from graphix.transpiler import (
+    Circuit,
+    OutputIndex,
+    OutputKind,
+    decompose_ccx,
+    decompose_cu,
+    decompose_p,
+    decompose_rx,
+    decompose_rz,
+    decompose_y,
+    insert_control,
+    instructions_to_jcz,
+    transpile_swaps,
+)
 from tests.test_branch_selector import CheckedBranchSelector
-from tests.test_instruction import VisitAngle
+from tests.test_instruction import INSTRUCTION_TEST_CASES, VisitAngle
 
 if TYPE_CHECKING:
-    from collections.abc import Callable
-    from typing import Literal, TypeAlias
+    from typing import Literal
 
-    from graphix.instruction import InstructionType
     from graphix.measurements import Outcome
+    from tests.test_instruction import InstructionTestCase
 
-    InstructionTestCase: TypeAlias = Callable[[Generator], InstructionType]
     _DenseStateBackendLiteral = Literal["statevector", "densitymatrix"]
 
 
-INSTRUCTION_TEST_CASES: list[InstructionTestCase] = [
-    lambda _rng: instruction.CCX(0, (1, 2)),
-    lambda rng: instruction.RZZ(0, 1, rng.random() * 2 * ANGLE_PI),
-    lambda _rng: instruction.CZ((0, 1)),
-    lambda _rng: instruction.CNOT(0, 1),
-    lambda _rng: instruction.SWAP((0, 1)),
-    lambda _rng: instruction.H(0),
-    lambda _rng: instruction.S(0),
-    lambda _rng: instruction.X(0),
-    lambda _rng: instruction.Y(0),
-    lambda _rng: instruction.Z(0),
-    lambda _rng: instruction.I(0),
-    lambda rng: instruction.RX(0, rng.random() * 2 * ANGLE_PI),
-    lambda rng: instruction.RY(0, rng.random() * 2 * ANGLE_PI),
-    lambda rng: instruction.RZ(0, rng.random() * 2 * ANGLE_PI),
-    lambda rng: instruction.J(0, rng.random() * 2 * ANGLE_PI),
-]
-
-
 class TestTranspilerUnitGates:
-    @pytest.mark.parametrize("instruction", INSTRUCTION_TEST_CASES)
-    def test_instruction_flow(self, fx_rng: Generator, instruction: InstructionTestCase) -> None:
-        circuit = Circuit(3, instr=[instruction(fx_rng)])
+    @pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+    def test_instruction_flow(self, fx_rng: Generator, test_case: InstructionTestCase) -> None:
+        circuit = Circuit(3, instr=[test_case.instruction(fx_rng)])
         pattern = circuit.transpile().pattern
         circuit.transpile_to_causalflow().flow.check_well_formed()
         flow = pattern.to_bloch().to_causalflow()
         flow.check_well_formed()
 
     @pytest.mark.parametrize("jumps", range(1, 11))
-    @pytest.mark.parametrize("instruction", INSTRUCTION_TEST_CASES)
-    def test_instructions(self, fx_bg: PCG64, jumps: int, instruction: InstructionTestCase) -> None:
+    @pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+    def test_instructions(self, fx_bg: PCG64, jumps: int, test_case: InstructionTestCase) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = Circuit(3, instr=[instruction(rng)])
+        circuit = Circuit(3, instr=[test_case.instruction(rng)])
         pattern = circuit.transpile().pattern
         input_state = rand_state_vector(3, rng=rng)
         state = circuit.simulate(input_state=input_state).state
@@ -414,3 +406,81 @@ def test_visit() -> None:
     assert circ.instruction != circ2.instruction
     assert circ.visit(visitor) is circ
     assert circ.instruction == circ2.instruction
+
+
+def test_transpile_cj(fx_rng: Generator) -> None:
+    alpha = fx_rng.random()
+    circuit = Circuit(2)
+    circuit.cj(0, 1, alpha)
+    decomposed_circuit = circuit.transpile_cj()
+    input_state = rand_state_vector(2, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)
+
+
+def test_decompose_cy(fx_rng: Generator) -> None:
+    circuit = Circuit(2)
+    circuit.cy(0, 1)
+    decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_y(Instruction.Y(1))))
+    input_state = rand_state_vector(2, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)
+
+
+def test_decompose_cp(fx_rng: Generator) -> None:
+    angle = fx_rng.random()
+    circuit = Circuit(2)
+    circuit.cp(0, 1, angle)
+    decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_p(Instruction.P(1, angle))))
+    input_state = rand_state_vector(2, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)
+
+
+def test_decompose_crx(fx_rng: Generator) -> None:
+    angle = fx_rng.random()
+    circuit = Circuit(2)
+    circuit.crx(0, 1, angle)
+    decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_rx(Instruction.RX(1, angle))))
+    input_state = rand_state_vector(2, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)
+
+
+def test_decompose_crz(fx_rng: Generator) -> None:
+    angle = fx_rng.random()
+    circuit = Circuit(2)
+    circuit.crz(0, 1, angle)
+    decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_rz(Instruction.RZ(1, angle))))
+    input_state = rand_state_vector(2, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)
+
+
+def test_decompose_cu(fx_rng: Generator) -> None:
+    theta = fx_rng.random()
+    phi = fx_rng.random()
+    lambda_ = fx_rng.random()
+    gamma = fx_rng.random()
+    circuit = Circuit(2)
+    circuit.cu(0, 1, theta, phi, lambda_, gamma)
+    decomposed_circuit = Circuit(2, instr=decompose_cu(Instruction.CU(0, 1, theta, phi, lambda_, gamma)))
+    input_state = rand_state_vector(2, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)
+
+
+@pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
+def test_instructions_to_jcz(fx_rng: Generator, test_case: InstructionTestCase) -> None:
+    circuit = Circuit(3, instr=[test_case.instruction(fx_rng)])
+    decomposed_circuit = Circuit(3, instr=instructions_to_jcz(circuit.instruction))
+    input_state = rand_state_vector(3, rng=fx_rng)
+    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
+    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
+    assert state.isclose(state2, atol=1e-15)

From 428b7cc538c4a93ddc1724dec3b16c01643fcbeb Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 25 Aug 2026 21:30:41 +0200
Subject: [PATCH 02/18] Update noxfile

---
 noxfile.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/noxfile.py b/noxfile.py
index 072b0ba7d..ee31e3c72 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -97,7 +97,7 @@ class ReverseDependency:
     [
         ReverseDependency("https://github.com/thierry-martinez/graphix-symbolic", branch="in-place_methods"),
         ReverseDependency("https://github.com/thierry-martinez/graphix-stim-backend", branch="rename-simulate"),
-        ReverseDependency("https://github.com/TeamGraphix/graphix-qasm-parser"),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-qasm-parser", branch="add_openqasm_gates"),
         ReverseDependency(
             "https://github.com/thierry-martinez/graphix-ibmq", doctest_modules=False, branch="rename-simulate"
         ),
@@ -109,7 +109,7 @@ class ReverseDependency:
             install_target=".[dev]",
             branch="rename-simulate",
         ),
-        ReverseDependency("https://github.com/thierry-martinez/graphix-mqtbench", branch="rename-simulate"),
+        ReverseDependency("https://github.com/thierry-martinez/graphix-mqtbench", branch="add_openqasm_gates"),
     ],
 )
 def tests_reverse_dependencies(session: Session, package: ReverseDependency) -> None:

From 575af833b05c043c15cc4020fd0e3e465d886b71 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 25 Aug 2026 22:07:56 +0200
Subject: [PATCH 03/18] Fix test coverage

---
 graphix/qasm3_exporter.py                     |  4 +-
 graphix/transpiler.py                         | 17 +-----
 tests/test_instruction.py                     |  7 ++-
 tests/test_qasm3_exporter.py                  | 28 ++++++++--
 .../test_qasm3_exporter_to_graphix_parser.py  | 29 ++++++++--
 tests/test_transpiler.py                      | 56 +++++++++----------
 6 files changed, 85 insertions(+), 56 deletions(-)

diff --git a/graphix/qasm3_exporter.py b/graphix/qasm3_exporter.py
index c6a8b4778..32af8af9d 100644
--- a/graphix/qasm3_exporter.py
+++ b/graphix/qasm3_exporter.py
@@ -183,7 +183,7 @@ def instruction_to_qasm3(instruction: InstructionType) -> str:
             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)])
+            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)
@@ -195,7 +195,7 @@ def instruction_to_qasm3(instruction: InstructionType) -> str:
                 operands=[qasm3_qubit(instruction.control), qasm3_qubit(instruction.target)],
             )
         case InstructionKind.GPHASE:
-            return qasm3_gate_call("gphase", [angle_to_qasm3(instruction.angle)])
+            return qasm3_gate_call("gphase", operands=[], args=[angle_to_qasm3(instruction.angle)])
         case _:
             assert_never(instruction.kind)
 
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index d3b3a972b..29c91f50b 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -1368,13 +1368,7 @@ def decompose_cu(instr: instruction.CU) -> Iterator[instruction.CJ | Instruction
 def insert_control(
     control: int,
     instrs: Iterable[
-        Instruction.GPHASE
-        | Instruction.X
-        | Instruction.Z
-        | Instruction.J
-        | Instruction.CZ
-        | Instruction.CNOT
-        | Instruction.RZ
+        Instruction.GPHASE | Instruction.X | Instruction.Z | Instruction.J | Instruction.CNOT | Instruction.RZ
     ],
 ) -> Iterable[InstructionType]:
     """Yield a controlled gate sequence from a gate sequence.
@@ -1383,8 +1377,8 @@ def insert_control(
     ----------
     control: int
         The control qubit.
-    instrs: Iterable[Instruction.J | Instruction.CZ]
-        The J-∧z decomposition.
+    instrs: Iterable[Instruction.X | Instruction.Z | Instruction.J | Instruction.CNOT | Instruction.RZ]
+        The gate sequence.
 
     Yields
     ------
@@ -1400,11 +1394,6 @@ def insert_control(
                 yield instruction.CZ((control, instr.target))
             case InstructionKind.J:
                 yield instruction.CJ(control=control, target=instr.target, angle=instr.angle)
-            case InstructionKind.CZ:
-                u, v = instr.targets
-                yield instruction.H(v)
-                yield instruction.CCX(target=v, controls=(control, u))
-                yield instruction.H(v)
             case InstructionKind.CNOT:
                 yield instruction.CCX(target=instr.target, controls=(control, instr.control))
             case InstructionKind.RZ:
diff --git a/tests/test_instruction.py b/tests/test_instruction.py
index c2c520a2c..32b18f2bc 100644
--- a/tests/test_instruction.py
+++ b/tests/test_instruction.py
@@ -76,6 +76,7 @@ class InstructionTestCase:
         ),
     ),
     InstructionTestCase("CSWAP", lambda _rng: Instruction.CSWAP(0, (1, 2))),
+    InstructionTestCase("GPHASE", lambda rng: Instruction.GPHASE(rng.random() * 2 * ANGLE_PI)),
 )
 
 
@@ -104,9 +105,11 @@ def test_visit_qubit(fx_rng: Generator, test_case: InstructionTestCase) -> None:
     visitor = VisitQubit()
     instr_visited = instr.visit(visitor, copy=True)
     assert instr == instr_copy
-    assert instr_visited != instr_copy
+    if test_case.name != "GPHASE":
+        assert instr_visited != instr_copy
     instr.visit(visitor, copy=False)
-    assert instr != instr_copy
+    if test_case.name != "GPHASE":
+        assert instr != instr_copy
     assert instr_visited == instr
 
 
diff --git a/tests/test_qasm3_exporter.py b/tests/test_qasm3_exporter.py
index 723e167c1..6c67be367 100644
--- a/tests/test_qasm3_exporter.py
+++ b/tests/test_qasm3_exporter.py
@@ -61,13 +61,33 @@ def test_to_qasm3_random_circuit(fx_bg: PCG64, jumps: int) -> None:
     _qasm3 = pattern_to_qasm3(pattern)
 
 
-def test_to_qasm3_failures() -> None:
-    circuit = Circuit(2)
+def test_to_qasm3_measure_on_x_axis() -> None:
+    circuit = Circuit(1)
     circuit.m(0, Axis.X)
     with pytest.raises(ValueError, match="OpenQASM3 only supports measurements on Z axis"):
         circuit_to_qasm3(circuit, transpile=False)
-    circuit = circuit.transpile_measurements_to_z_axis()
-    circuit.j(1, 0.25)
+    _qasm3 = circuit_to_qasm3(circuit)
+
+
+def test_to_qasm3_j() -> None:
+    circuit = Circuit(1)
+    circuit.j(0, 0.25)
     with pytest.raises(ValueError, match="J gates must be decomposed before QASM3 export"):
         circuit_to_qasm3(circuit, transpile=False)
     _qasm3 = circuit_to_qasm3(circuit)
+
+
+def test_to_qasm3_cj() -> None:
+    circuit = Circuit(2)
+    circuit.cj(0, 1, 0.25)
+    with pytest.raises(ValueError, match="CJ gates must be decomposed before QASM3 export"):
+        circuit_to_qasm3(circuit, transpile=False)
+    _qasm3 = circuit_to_qasm3(circuit)
+
+
+def test_to_qasm3_rzz() -> None:
+    circuit = Circuit(2)
+    circuit.rzz(0, 1, 0.25)
+    with pytest.raises(ValueError, match="RZZ gates must be decomposed before QASM3 export"):
+        circuit_to_qasm3(circuit, transpile=False)
+    _qasm3 = circuit_to_qasm3(circuit)
diff --git a/tests/test_qasm3_exporter_to_graphix_parser.py b/tests/test_qasm3_exporter_to_graphix_parser.py
index 8133c18ff..727d3b24a 100644
--- a/tests/test_qasm3_exporter_to_graphix_parser.py
+++ b/tests/test_qasm3_exporter_to_graphix_parser.py
@@ -65,14 +65,33 @@ def test_instruction_to_qasm3(fx_rng: Generator, test_case: InstructionTestCase)
 
 
 def test_j_to_qasm3() -> None:
-    circuit = Circuit(3, instr=[Instruction.J(target=0, angle=ANGLE_PI / 4)])
+    circuit = Circuit(1, instr=[Instruction.J(target=0, angle=ANGLE_PI / 4)])
     qasm = circuit_to_qasm3(circuit)
     parser = OpenQASMParser()
     parsed_circuit = parser.parse_str(qasm)
     assert parsed_circuit.instruction == circuit.transpile_j_to_rzh().instruction
 
 
-def test_j_to_qasm3_failure() -> None:
-    circuit = Circuit(3, instr=[Instruction.J(target=0, angle=ANGLE_PI / 4)])
-    with pytest.raises(ValueError):
-        circuit_to_qasm3(circuit, transpile=False)
+def test_cj_to_qasm3() -> None:
+    circuit = Circuit(2, instr=[Instruction.CJ(control=0, target=1, angle=ANGLE_PI / 4)])
+    qasm = circuit_to_qasm3(circuit)
+    parser = OpenQASMParser()
+    parsed_circuit = parser.parse_str(qasm)
+    assert parsed_circuit.instruction == circuit.transpile_cj().instruction
+
+
+def test_rzz_to_qasm3() -> None:
+    circuit = Circuit(2, instr=[Instruction.RZZ(control=0, target=1, angle=ANGLE_PI / 4)])
+    qasm = circuit_to_qasm3(circuit)
+    parser = OpenQASMParser()
+    parsed_circuit = parser.parse_str(qasm)
+    assert parsed_circuit.instruction == circuit.transpile_rzz().instruction
+
+
+def test_gphase_to_qasm3() -> None:
+    instr = Instruction.GPHASE(ANGLE_PI / 4)
+    circuit = Circuit(1, instr=[instr])
+    qasm = circuit_to_qasm3(circuit)
+    parser = OpenQASMParser()
+    parsed_circuit = parser.parse_str(qasm)
+    assert parsed_circuit.instruction == [instr]
diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py
index f41ff6dc4..696a19e71 100644
--- a/tests/test_transpiler.py
+++ b/tests/test_transpiler.py
@@ -9,7 +9,7 @@
 
 from graphix import Instruction, instruction
 from graphix.branch_selector import ConstBranchSelector, FixedBranchSelector
-from graphix.fundamentals import Axis, Sign
+from graphix.fundamentals import ANGLE_PI, Axis, Sign
 from graphix.instruction import I, InstructionKind
 from graphix.random_objects import rand_circuit, rand_gate, rand_state_vector
 from graphix.sim.density_matrix import DensityMatrix
@@ -408,25 +408,26 @@ def test_visit() -> None:
     assert circ.instruction == circ2.instruction
 
 
+def check_circuit_equivalence(circuit1: Circuit, circuit2: Circuit, rng: Generator) -> bool:
+    input_state = rand_state_vector(circuit1.width, rng=rng)
+    state1 = circuit1.simulate(input_state=input_state, rng=rng).state
+    state2 = circuit2.simulate(input_state=input_state, rng=rng).state
+    return state1.isclose(state2, atol=1e-15)
+
+
 def test_transpile_cj(fx_rng: Generator) -> None:
     alpha = fx_rng.random()
     circuit = Circuit(2)
     circuit.cj(0, 1, alpha)
     decomposed_circuit = circuit.transpile_cj()
-    input_state = rand_state_vector(2, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 
 def test_decompose_cy(fx_rng: Generator) -> None:
     circuit = Circuit(2)
     circuit.cy(0, 1)
     decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_y(Instruction.Y(1))))
-    input_state = rand_state_vector(2, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 
 def test_decompose_cp(fx_rng: Generator) -> None:
@@ -434,10 +435,7 @@ def test_decompose_cp(fx_rng: Generator) -> None:
     circuit = Circuit(2)
     circuit.cp(0, 1, angle)
     decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_p(Instruction.P(1, angle))))
-    input_state = rand_state_vector(2, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 
 def test_decompose_crx(fx_rng: Generator) -> None:
@@ -445,10 +443,7 @@ def test_decompose_crx(fx_rng: Generator) -> None:
     circuit = Circuit(2)
     circuit.crx(0, 1, angle)
     decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_rx(Instruction.RX(1, angle))))
-    input_state = rand_state_vector(2, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 
 def test_decompose_crz(fx_rng: Generator) -> None:
@@ -456,10 +451,7 @@ def test_decompose_crz(fx_rng: Generator) -> None:
     circuit = Circuit(2)
     circuit.crz(0, 1, angle)
     decomposed_circuit = Circuit(2, instr=insert_control(0, decompose_rz(Instruction.RZ(1, angle))))
-    input_state = rand_state_vector(2, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 
 def test_decompose_cu(fx_rng: Generator) -> None:
@@ -470,17 +462,23 @@ def test_decompose_cu(fx_rng: Generator) -> None:
     circuit = Circuit(2)
     circuit.cu(0, 1, theta, phi, lambda_, gamma)
     decomposed_circuit = Circuit(2, instr=decompose_cu(Instruction.CU(0, 1, theta, phi, lambda_, gamma)))
-    input_state = rand_state_vector(2, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 
 @pytest.mark.parametrize("test_case", INSTRUCTION_TEST_CASES)
 def test_instructions_to_jcz(fx_rng: Generator, test_case: InstructionTestCase) -> None:
     circuit = Circuit(3, instr=[test_case.instruction(fx_rng)])
     decomposed_circuit = Circuit(3, instr=instructions_to_jcz(circuit.instruction))
-    input_state = rand_state_vector(3, rng=fx_rng)
-    state = circuit.simulate(input_state=input_state, rng=fx_rng).state
-    state2 = decomposed_circuit.simulate(input_state=input_state, rng=fx_rng).state
-    assert state.isclose(state2, atol=1e-15)
+    assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
+
+
+def test_cr() -> None:
+    circuit = Circuit(2)
+    circuit.cr(control=0, target=1, axis=Axis.X, angle=ANGLE_PI / 2)
+    circuit.cr(control=1, target=0, axis=Axis.Y, angle=ANGLE_PI / 4)
+    circuit.cr(control=0, target=1, axis=Axis.Z, angle=ANGLE_PI / 8)
+    assert circuit.instruction == [
+        Instruction.CRX(control=0, target=1, angle=ANGLE_PI / 2),
+        Instruction.CRY(control=1, target=0, angle=ANGLE_PI / 4),
+        Instruction.CRZ(control=0, target=1, angle=ANGLE_PI / 8),
+    ]

From bbf20ce3a9fa5d3f34ae04df66ff261375d49b89 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 25 Aug 2026 22:16:11 +0200
Subject: [PATCH 04/18] Fix Qiskit test

---
 tests/test_qasm3_exporter_to_qiskit.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_qasm3_exporter_to_qiskit.py b/tests/test_qasm3_exporter_to_qiskit.py
index c8f107dde..5e30bee90 100644
--- a/tests/test_qasm3_exporter_to_qiskit.py
+++ b/tests/test_qasm3_exporter_to_qiskit.py
@@ -152,6 +152,6 @@ def test_instruction_to_qasm3(fx_rng: Generator, test_case: InstructionTestCase)
     instr = test_case.instruction(fx_rng)
     if instr.kind in {InstructionKind.CJ, InstructionKind.RZZ, InstructionKind.M}:
         pytest.skip()
-    if instr.kind in {InstructionKind.SXDG, InstructionKind.U}:
-        pytest.skip("qiskit_qasm3_import.exceptions.ConversionError: gate 'sxdg'/'u' is not defined.")
+    if instr.kind == InstructionKind.SXDG:
+        pytest.skip("qiskit_qasm3_import.exceptions.ConversionError: gate 'sxdg' is not defined.")
     check_qasm3_circuit(Circuit(3, instr=[instr]))

From a3e47548aeca31ab4b4e335e25aee305a7217462 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 25 Aug 2026 22:33:30 +0200
Subject: [PATCH 05/18] Fix qasm plugin version installed by nox

---
 noxfile.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/noxfile.py b/noxfile.py
index ee31e3c72..dde2b5cf4 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -51,7 +51,7 @@ def tests_all(session: Session) -> None:
     """Run the test suite with all dependencies."""
     session.install(".[dev]")
     # This dependency is added here to avoid circular dependencies
-    session.install("graphix-qasm-parser>=0.1.1")
+    session.install("graphix-qasm-parser@git+https://github.com/thierry-martinez/graphix-qasm-parser@add_openqasm_gates")
     run_pytest(session, doctest_modules=True, mpl=True)
 
 

From 147eaac1122a8357be708241511916fbc694cba3 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 25 Aug 2026 22:37:51 +0200
Subject: [PATCH 06/18] ruff

---
 noxfile.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/noxfile.py b/noxfile.py
index dde2b5cf4..94e420451 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -51,7 +51,9 @@ def tests_all(session: Session) -> None:
     """Run the test suite with all dependencies."""
     session.install(".[dev]")
     # This dependency is added here to avoid circular dependencies
-    session.install("graphix-qasm-parser@git+https://github.com/thierry-martinez/graphix-qasm-parser@add_openqasm_gates")
+    session.install(
+        "graphix-qasm-parser@git+https://github.com/thierry-martinez/graphix-qasm-parser@add_openqasm_gates"
+    )
     run_pytest(session, doctest_modules=True, mpl=True)
 
 

From 6672b2127babcdfcba1ee3b3f57957917bb79b92 Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Mon, 31 Aug 2026 23:01:00 +0200
Subject: [PATCH 07/18] Update graphix/ops.py

Co-authored-by: matulni <m.uldemolins@gmail.com>
---
 graphix/ops.py | 13 +------------
 1 file changed, 1 insertion(+), 12 deletions(-)

diff --git a/graphix/ops.py b/graphix/ops.py
index 5cfcecb7f..1ac2b3daf 100644
--- a/graphix/ops.py
+++ b/graphix/ops.py
@@ -276,19 +276,8 @@ def cu(
         -------
         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,
-        )
+        return controlled(exp(1j * gamma_rad) * Ops.u(theta, phi, lambda_))
 
     CH: ClassVar[npt.NDArray[np.complex128]] = controlled(H)
 

From d5141eaab90effd8ec8aff1f98366d9ec27d122a Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Mon, 31 Aug 2026 23:40:18 +0200
Subject: [PATCH 08/18] Update graphix/transpiler.py

Co-authored-by: matulni <m.uldemolins@gmail.com>
---
 graphix/transpiler.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 29c91f50b..e3271d61f 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -1405,7 +1405,7 @@ def insert_control(
     yield Instruction.P(target=control, angle=gphase)
 
 
-def decompose_cj(instr: Instruction.CJ) -> Iterator[InstructionType]:
+def decompose_cj(instr: Instruction.CJ) -> Iterator[Instruction.RZ | Instruction.CNOT | Instruction.RY | Instruction.P]:
     """Yield a decomposed gate sequence of the CJ gate.
 
     See :class:`~graphix.instruction.CJ` for more information.

From 93200dbfa3860c8bb6310cb9efe2093e7abe7a17 Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Mon, 31 Aug 2026 23:41:12 +0200
Subject: [PATCH 09/18] Update graphix/transpiler.py

Co-authored-by: matulni <m.uldemolins@gmail.com>
---
 graphix/transpiler.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index e3271d61f..1c60670e9 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -1377,7 +1377,7 @@ def insert_control(
     ----------
     control: int
         The control qubit.
-    instrs: Iterable[Instruction.X | Instruction.Z | Instruction.J | Instruction.CNOT | Instruction.RZ]
+    instrs: Iterable[Instruction.GPHASE | Instruction.X | Instruction.Z | Instruction.J | Instruction.CNOT | Instruction.RZ]
         The gate sequence.
 
     Yields

From 8aea0e36c9b8c2cc9ca0693ec66f44a78ca4f60c Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Mon, 31 Aug 2026 23:41:42 +0200
Subject: [PATCH 10/18] Update graphix/transpiler.py

Co-authored-by: matulni <m.uldemolins@gmail.com>
---
 graphix/transpiler.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 1c60670e9..c6341605b 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -1370,7 +1370,7 @@ def insert_control(
     instrs: Iterable[
         Instruction.GPHASE | Instruction.X | Instruction.Z | Instruction.J | Instruction.CNOT | Instruction.RZ
     ],
-) -> Iterable[InstructionType]:
+) -> Iterator[Instruction.CNOT | Instruction.CZ | Instruction.CJ | Instruction.CCX | Instruction.CRZ | Instruction.P]:
     """Yield a controlled gate sequence from a gate sequence.
 
     Parameters

From bcf376add9ca18c79df00f74be7bee8aa43e1932 Mon Sep 17 00:00:00 2001
From: thierry-martinez <thierry.martinez@inria.fr>
Date: Mon, 31 Aug 2026 23:43:14 +0200
Subject: [PATCH 11/18] Update graphix/transpiler.py

Co-authored-by: matulni <m.uldemolins@gmail.com>
---
 graphix/transpiler.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index c6341605b..d7ba7b2ea 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -1326,7 +1326,8 @@ def decompose_rz(instr: instruction.RZ) -> Iterator[instruction.J | Instruction.
 def decompose_u(instr: instruction.U) -> Iterator[instruction.J | Instruction.GPHASE]:
     """Yield a J decomposition of the U gate.
 
-    The U(θ, φ, λ) gate is decomposed into H·J(φ + 𝜋/2)·J(θ)·J(λ - 𝜋/2) (that is to say, J(0)·J(φ + 𝜋/2)·J(θ)·J(λ - 𝜋/2)).
+    The :math:`U(\theta, \phi, \lambda)` gate is decomposed as :math:`e^{-i \theta/2} \cdot H \cdot J(\phi + \pi/2) \cdot J(\theta) \cdot J(\lambda - \pi/2)`, or equivalently, :math:`J(0) \cdot J(\phi + \pi/2) \cdot J(\theta) \cdot J(\lambda - \pi/2)`.
+
 
     Parameters
     ----------

From 06bb6f267d0461da8f7ed4808aa8f072f393f9ad Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 1 Sep 2026 00:00:12 +0200
Subject: [PATCH 12/18] Implement Mateo's suggestions

---
 graphix/circ_ext/compilation.py               | 10 +-
 graphix/instruction.py                        |  9 --
 graphix/ops.py                                |  8 +-
 graphix/qasm3_exporter.py                     |  2 +-
 graphix/transpiler.py                         | 91 ++++++-------------
 .../test_qasm3_exporter_to_graphix_parser.py  |  8 +-
 tests/test_transpiler.py                      |  6 +-
 7 files changed, 47 insertions(+), 87 deletions(-)

diff --git a/graphix/circ_ext/compilation.py b/graphix/circ_ext/compilation.py
index 0aea02d43..661796968 100644
--- a/graphix/circ_ext/compilation.py
+++ b/graphix/circ_ext/compilation.py
@@ -8,7 +8,7 @@
 import numpy as np
 
 from graphix.fundamentals import ANGLE_PI, Axis
-from graphix.instruction import CNOT, SWAP, H, S, X, Y, Z
+from graphix.instruction import CNOT, SDG, SWAP, H, X, Y, Z
 from graphix.transpiler import Circuit
 
 if TYPE_CHECKING:
@@ -200,7 +200,7 @@ def cm_berg_pass(clifford_map: CliffordMap, circuit: Circuit) -> None:
     -----
     This pass only handles unitaries so far (Clifford maps with the same number of input and output nodes).
 
-    Gate set: H, S, CNOT, SWAP, X, Y, Z
+    Gate set: H, SDG, CNOT, SWAP, X, Y, Z
 
     This function converts a ``CliffordMap`` into a sequence of quantum
     gate instructions by operating on its binary tableau representation.
@@ -280,7 +280,7 @@ def do_step_1(tab: MatGF2, instructions: list[InstructionType], row_idx: int) ->
         col_idx_zx = np.flatnonzero(tab[row_idx, n : 2 * n])  # Don't take the sign column
         for j in col_idx_zx:
             # Each iteration sets the element `tab[row_idx, n+j]` to 0.
-            add_s(tab, instructions, int(j)) if tab[row_idx, j] else add_h(tab, instructions, int(j))
+            add_sdg(tab, instructions, int(j)) if tab[row_idx, j] else add_h(tab, instructions, int(j))
 
     def do_step_2(tab: MatGF2, instructions: list[InstructionType], row_idx: int) -> int:
         col_idx_xx = np.flatnonzero(tab[row_idx, :n])
@@ -298,11 +298,11 @@ def add_h(tab: MatGF2, instructions: list[InstructionType], q: int) -> None:
         tab[:, [q, q + n]] = tab[:, [q + n, q]]  # The usual tuple assignment `a, b = b, a` does not work here.
         instructions.append(H(q))
 
-    def add_s(tab: MatGF2, instructions: list[InstructionType], q: int) -> None:
+    def add_sdg(tab: MatGF2, instructions: list[InstructionType], q: int) -> None:
         tab[:, -1] ^= tab[:, q] & tab[:, q + n]
         tab[:, q + n] ^= tab[:, q]
         q = int(q)
-        instructions.extend((S(q), Z(q)))  # We append Sdagger to get C instead of C^dagger
+        instructions.append(SDG(q))  # We append S^dagger to get C instead of C^dagger
 
     def add_cnot(tab: MatGF2, instructions: list[InstructionType], qc: int, qt: int) -> None:
         tab[:, -1] ^= tab[:, qc] & tab[:, qt + n] & (tab[:, qt] ^ tab[:, qc + n] ^ 1)
diff --git a/graphix/instruction.py b/graphix/instruction.py
index eb5927024..1a0c68bf4 100644
--- a/graphix/instruction.py
+++ b/graphix/instruction.py
@@ -248,15 +248,6 @@ class CSWAP(_KindChecker, BaseInstruction):
 
     The CSWAP gate applies the matrix
 
-    .. math::
-
-      \left[\begin{matrix}
-        1 & 0 & 0 & 0\\
-        0 & 1 & 0 & 0\\
-        0 & 0 & \cos \frac \theta 2 & -\mathrm i \sin \frac \theta 2\\
-        0 & 0 & -\mathrm i \sin \frac \theta 2 & \cos \frac \theta 2
-      \end{matrix}\right]
-
     .. math::
 
       \left[\begin{matrix}
diff --git a/graphix/ops.py b/graphix/ops.py
index 1ac2b3daf..9d4bfd864 100644
--- a/graphix/ops.py
+++ b/graphix/ops.py
@@ -114,6 +114,10 @@ class Ops:
         )
     )
 
+    CH: ClassVar[npt.NDArray[np.complex128]] = controlled(H)
+
+    CSWAP: ClassVar[npt.NDArray[np.complex128]] = controlled(SWAP)
+
     @overload
     @staticmethod
     def _cast_array(array: Iterable[Iterable[complex]], theta: Angle) -> npt.NDArray[np.complex128]: ...
@@ -279,10 +283,6 @@ def cu(
         gamma_rad = angle_to_rad(gamma)
         return controlled(exp(1j * gamma_rad) * Ops.u(theta, phi, lambda_))
 
-    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]: ...
diff --git a/graphix/qasm3_exporter.py b/graphix/qasm3_exporter.py
index 32af8af9d..11cd9d358 100644
--- a/graphix/qasm3_exporter.py
+++ b/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_cj().transpile_rzz().transpile_j_to_rzh().transpile_measurements_to_z_axis()
+        circuit = circuit.transpile_to_qasm_gates()
     yield "OPENQASM 3;"
     yield 'include "stdgates.inc";'
     yield f"qubit[{circuit.width}] q;"
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index d7ba7b2ea..4785625b8 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -929,8 +929,6 @@ def simulate(
 
         classical_measures: list[Outcome] = []
 
-        gphase: ParameterizedAngle = 0
-
         for i in range(len(self.instruction)):
             instr = self.instruction[i]
 
@@ -1009,10 +1007,10 @@ def evolve(op: Matrix, qargs: Iterable[int]) -> None:
                     )
                     classical_measures.append(result)
                 case InstructionKind.GPHASE:
-                    gphase += instr.angle
+                    # Global phase is currently ignored
+                    pass
                 case _:
                     assert_never(instr.kind)
-        # Global phase is currently ignored
         return SimulateResult(_backend.state, tuple(classical_measures))
 
     def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Circuit:
@@ -1091,56 +1089,30 @@ def replace_parameters(
     ) -> Circuit:
         return self.apply_angle(lambda angle: parameter.with_parameters(angle, assignment), copy=copy)
 
-    def transpile_measurements_to_z_axis(self) -> Circuit:
-        """Return an equivalent circuit where all measurements are on Z axis."""
-        circuit = Circuit(width=self.width)
-        for instr in self.instruction:
-            if instr.kind == InstructionKind.M:
-                match instr.axis:
-                    case Axis.X:
-                        circuit.h(instr.target)
-                        circuit.m(instr.target, Axis.Z)
-                    case Axis.Y:
-                        circuit.rx(instr.target, ANGLE_PI / 2)
-                        circuit.m(instr.target, Axis.Z)
-                    case Axis.Z:
-                        circuit.add(instr)
-                    case _:
-                        assert_never(instr.axis)
-            else:
-                circuit.add(instr)
-        return circuit
-
-    def transpile_j_to_rzh(self) -> Circuit:
-        """Return an equivalent circuit where all J gates have been replaced with RZ and H gates."""
+    def transpile_to_qasm_gates(self) -> Circuit:
+        """Return an equivalent circuit using only the standard OpenQASM gate set."""
         new_circuit = Circuit(self.width)
         for instr in self.instruction:
             match instr.kind:
                 case InstructionKind.J:
                     new_circuit.add(instruction.RZ(target=instr.target, angle=instr.angle))
                     new_circuit.add(instruction.H(target=instr.target))
-                case _:
-                    new_circuit.add(instr)
-        return new_circuit
-
-    def transpile_cj(self) -> Circuit:
-        """Return an equivalent circuit where all CJ gates have been replaced with OpenQASM gates."""
-        new_circuit = Circuit(self.width)
-        for instr in self.instruction:
-            match instr.kind:
                 case InstructionKind.CJ:
                     new_circuit.extend(decompose_cj(instr))
-                case _:
-                    new_circuit.add(instr)
-        return new_circuit
-
-    def transpile_rzz(self) -> Circuit:
-        """Return an equivalent circuit where all RZZ gates have been replaced with OpenQASM gates."""
-        new_circuit = Circuit(self.width)
-        for instr in self.instruction:
-            match instr.kind:
                 case InstructionKind.RZZ:
                     new_circuit.extend(decompose_rzz(instr))
+                case InstructionKind.M:
+                    match instr.axis:
+                        case Axis.X:
+                            new_circuit.h(instr.target)
+                            new_circuit.m(instr.target, Axis.Z)
+                        case Axis.Y:
+                            new_circuit.rx(instr.target, ANGLE_PI / 2)
+                            new_circuit.m(instr.target, Axis.Z)
+                        case Axis.Z:
+                            new_circuit.add(instr)
+                        case _:
+                            assert_never(instr.axis)
                 case _:
                     new_circuit.add(instr)
         return new_circuit
@@ -1386,24 +1358,22 @@ def insert_control(
     InstructionType
         The controlled gate sequence.
     """
-    gphase: ParameterizedAngle = 0
     for instr in instrs:
         match instr.kind:
             case InstructionKind.X:
-                yield instruction.CNOT(control=control, target=instr.target)
+                yield Instruction.CNOT(control=control, target=instr.target)
             case InstructionKind.Z:
-                yield instruction.CZ((control, instr.target))
+                yield Instruction.CZ((control, instr.target))
             case InstructionKind.J:
-                yield instruction.CJ(control=control, target=instr.target, angle=instr.angle)
+                yield Instruction.CJ(control=control, target=instr.target, angle=instr.angle)
             case InstructionKind.CNOT:
-                yield instruction.CCX(target=instr.target, controls=(control, instr.control))
+                yield Instruction.CCX(target=instr.target, controls=(control, instr.control))
             case InstructionKind.RZ:
-                yield instruction.CRZ(control=control, target=instr.target, angle=instr.angle)
+                yield Instruction.CRZ(control=control, target=instr.target, angle=instr.angle)
             case InstructionKind.GPHASE:
-                gphase += instr.angle
+                yield Instruction.P(target=control, angle=instr.angle)
             case _:
                 assert_never(instr.kind)
-    yield Instruction.P(target=control, angle=gphase)
 
 
 def decompose_cj(instr: Instruction.CJ) -> Iterator[Instruction.RZ | Instruction.CNOT | Instruction.RY | Instruction.P]:
@@ -1412,13 +1382,13 @@ def decompose_cj(instr: Instruction.CJ) -> Iterator[Instruction.RZ | Instruction
     See :class:`~graphix.instruction.CJ` for more information.
     """
     delta = (instr.angle + ANGLE_PI) / 2
-    yield instruction.RZ(target=instr.target, angle=delta)
-    yield instruction.CNOT(control=instr.control, target=instr.target)
-    yield instruction.RZ(target=instr.target, angle=-delta)
-    yield instruction.RY(target=instr.target, angle=-ANGLE_PI / 4)
-    yield instruction.CNOT(control=instr.control, target=instr.target)
-    yield instruction.RY(target=instr.target, angle=ANGLE_PI / 4)
-    yield instruction.P(target=instr.control, angle=delta)
+    yield Instruction.RZ(target=instr.target, angle=delta)
+    yield Instruction.CNOT(control=instr.control, target=instr.target)
+    yield Instruction.RZ(target=instr.target, angle=-delta)
+    yield Instruction.RY(target=instr.target, angle=-ANGLE_PI / 4)
+    yield Instruction.CNOT(control=instr.control, target=instr.target)
+    yield Instruction.RY(target=instr.target, angle=ANGLE_PI / 4)
+    yield Instruction.P(target=instr.control, angle=delta)
 
 
 def decompose_p(instr: Instruction.P) -> Iterator[Instruction.RZ | Instruction.GPHASE]:
@@ -1515,8 +1485,7 @@ def instructions_to_jcz(
                     insert_control(instr.control, decompose_swap(Instruction.SWAP(instr.targets)))
                 )
             case InstructionKind.GPHASE:
-                # Global phase is currently ignored
-                pass
+                yield instr
             case _:
                 assert_never(instr.kind)
 
diff --git a/tests/test_qasm3_exporter_to_graphix_parser.py b/tests/test_qasm3_exporter_to_graphix_parser.py
index 727d3b24a..4b9f19649 100644
--- a/tests/test_qasm3_exporter_to_graphix_parser.py
+++ b/tests/test_qasm3_exporter_to_graphix_parser.py
@@ -35,7 +35,7 @@
 
 def check_round_trip(circuit: Circuit) -> None:
     qasm = circuit_to_qasm3(circuit)
-    check_circuit = circuit.transpile_j_to_rzh()
+    check_circuit = circuit.transpile_to_qasm_gates()
     parser = OpenQASMParser()
     parsed_circuit = parser.parse_str(qasm)
     for parsed_instr, instr in zip(parsed_circuit.instruction, check_circuit.instruction, strict=True):
@@ -69,7 +69,7 @@ def test_j_to_qasm3() -> None:
     qasm = circuit_to_qasm3(circuit)
     parser = OpenQASMParser()
     parsed_circuit = parser.parse_str(qasm)
-    assert parsed_circuit.instruction == circuit.transpile_j_to_rzh().instruction
+    assert parsed_circuit.instruction == circuit.transpile_to_qasm_gates().instruction
 
 
 def test_cj_to_qasm3() -> None:
@@ -77,7 +77,7 @@ def test_cj_to_qasm3() -> None:
     qasm = circuit_to_qasm3(circuit)
     parser = OpenQASMParser()
     parsed_circuit = parser.parse_str(qasm)
-    assert parsed_circuit.instruction == circuit.transpile_cj().instruction
+    assert parsed_circuit.instruction == circuit.transpile_to_qasm_gates().instruction
 
 
 def test_rzz_to_qasm3() -> None:
@@ -85,7 +85,7 @@ def test_rzz_to_qasm3() -> None:
     qasm = circuit_to_qasm3(circuit)
     parser = OpenQASMParser()
     parsed_circuit = parser.parse_str(qasm)
-    assert parsed_circuit.instruction == circuit.transpile_rzz().instruction
+    assert parsed_circuit.instruction == circuit.transpile_to_qasm_gates().instruction
 
 
 def test_gphase_to_qasm3() -> None:
diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py
index 696a19e71..83142f5ff 100644
--- a/tests/test_transpiler.py
+++ b/tests/test_transpiler.py
@@ -147,7 +147,7 @@ def test_transpile_measurements_to_z_axis(self, fx_bg: PCG64, jumps: int, axis:
         input_state = rand_state_vector(2, rng=rng)
         branch_selector = ConstBranchSelector(outcome)
         state = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).state
-        circuit_z = circuit.transpile_measurements_to_z_axis()
+        circuit_z = circuit.transpile_to_qasm_gates()
         assert all(instr.axis == Axis.Z for instr in circuit_z.instruction if instr.kind == InstructionKind.M)
         state_z = circuit.simulate(rng=rng, input_state=input_state, branch_selector=branch_selector).state
         assert state_z.isclose(state)
@@ -160,7 +160,7 @@ def test_transpile_j_to_rzh(self, fx_bg: PCG64, jumps: int) -> None:
         circuit = rand_circuit(nqubits, depth, rng, use_j=True, use_ccx=True, use_rzz=True)
         circuit.j(0, 0.5)  # Ensure that there is at least one J instruction
         assert any(instr.kind == InstructionKind.J for instr in circuit.instruction)
-        circuit2 = circuit.transpile_j_to_rzh()
+        circuit2 = circuit.transpile_to_qasm_gates()
         assert not any(instr.kind == InstructionKind.J for instr in circuit2.instruction)
         state = circuit.simulate(rng=rng).state
         state2 = circuit2.simulate(rng=rng).state
@@ -419,7 +419,7 @@ def test_transpile_cj(fx_rng: Generator) -> None:
     alpha = fx_rng.random()
     circuit = Circuit(2)
     circuit.cj(0, 1, alpha)
-    decomposed_circuit = circuit.transpile_cj()
+    decomposed_circuit = circuit.transpile_to_qasm_gates()
     assert check_circuit_equivalence(circuit, decomposed_circuit, rng=fx_rng)
 
 

From 70ed3ffa6b6d4659253f76aee1260ff7bcfae2a3 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 1 Sep 2026 07:26:31 +0200
Subject: [PATCH 13/18] Fix syntax and mypy

---
 graphix/parameter.py  | 39 +++++++++++++++++++++++++++++++++++++--
 graphix/transpiler.py |  2 +-
 2 files changed, 38 insertions(+), 3 deletions(-)

diff --git a/graphix/parameter.py b/graphix/parameter.py
index 891986f38..b813b6c37 100644
--- a/graphix/parameter.py
+++ b/graphix/parameter.py
@@ -21,12 +21,32 @@
     from collections.abc import Mapping
     from typing import Self
 
+    import numpy as np
+    from numpy.typing import NDArray
+
+    from graphix.fundamentals import Sign
+
 
 class Expression(ABC):
     """Expression with parameters."""
 
+    @overload
+    def __mul__(self, other: float) -> ExpressionOrFloat: ...
+
+    @overload
+    def __mul__(self, other: Sign) -> ExpressionOrFloat: ...
+
+    @overload
+    def __mul__(self, other: NDArray[np.object_]) -> NDArray[np.object_]: ...
+
+    @overload
+    def __mul__(self, other: NDArray[np.complex128]) -> NDArray[np.complex128]: ...
+
+    @overload
+    def __mul__(self, other: Expression) -> ExpressionOrFloat: ...
+
     @abstractmethod
-    def __mul__(self, other: object) -> ExpressionOrFloat:
+    def __mul__(self, other: object) -> ExpressionOrFloat | NDArray[np.object_] | NDArray[np.complex128]:
         """
         Return the product of this expression with another object.
 
@@ -178,8 +198,23 @@ def scale(self, k: float) -> ExpressionOrFloat:
             return 0
         return self.scale_non_null(k)
 
+    @overload
+    def __mul__(self, other: float) -> ExpressionOrFloat: ...
+
+    @overload
+    def __mul__(self, other: Sign) -> ExpressionOrFloat: ...
+
+    @overload
+    def __mul__(self, other: NDArray[np.object_]) -> NDArray[np.object_]: ...
+
+    @overload
+    def __mul__(self, other: NDArray[np.complex128]) -> NDArray[np.complex128]: ...
+
+    @overload
+    def __mul__(self, other: Expression) -> ExpressionOrFloat: ...
+
     @override
-    def __mul__(self, other: object) -> ExpressionOrFloat:
+    def __mul__(self, other: object) -> ExpressionOrFloat | NDArray[np.object_] | NDArray[np.complex128]:
         """Look to the documentation in the parent class."""
         if isinstance(other, SupportsFloat):
             return self.scale(float(other))
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 4785625b8..357eb4748 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -1296,7 +1296,7 @@ def decompose_rz(instr: instruction.RZ) -> Iterator[instruction.J | Instruction.
 
 
 def decompose_u(instr: instruction.U) -> Iterator[instruction.J | Instruction.GPHASE]:
-    """Yield a J decomposition of the U gate.
+    r"""Yield a J decomposition of the U gate.
 
     The :math:`U(\theta, \phi, \lambda)` gate is decomposed as :math:`e^{-i \theta/2} \cdot H \cdot J(\phi + \pi/2) \cdot J(\theta) \cdot J(\lambda - \pi/2)`, or equivalently, :math:`J(0) \cdot J(\phi + \pi/2) \cdot J(\theta) \cdot J(\lambda - \pi/2)`.
 

From 9e475238369751554ec17a1f751dc4ffab4238b5 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Tue, 1 Sep 2026 22:38:59 +0200
Subject: [PATCH 14/18] Document the basis and better typing

---
 graphix/instruction.py | 108 +++++++++++++++++++++++++++++++++++++++--
 graphix/ops.py         |  26 +++++-----
 2 files changed, 119 insertions(+), 15 deletions(-)

diff --git a/graphix/instruction.py b/graphix/instruction.py
index 1a0c68bf4..633bc0871 100644
--- a/graphix/instruction.py
+++ b/graphix/instruction.py
@@ -260,6 +260,17 @@ class CSWAP(_KindChecker, BaseInstruction):
         0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\
         0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``targets[0]``, ``targets[1]``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    targets : tuple[int, int]
+        Indices of the two target qubits.
     """
 
     control: int
@@ -550,14 +561,32 @@ class CU(_KindChecker, BaseInstruction):
           \cos\left(\frac{\theta}{2}\right)
       \end{matrix}\right]
 
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
     It can be decomposed as
 
     .. math::
-
       CU(\theta, \phi, \lambda, \gamma) =
-        \left(P\left(\frac{\gamma - \theta} 2\right) \otimes I)
+        \left(P\left(\frac{\gamma - \theta} 2\right) \otimes I\right)
         CJ(0) CJ\left(\phi + \frac \pi 2\right)
         CJ(\theta) CJ\left(\lambda - \frac \pi 2\right)
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    theta : ParameterizedAngle
+        Rotation angle around the Y axis.
+    phi : ParameterizedAngle
+        Rotation angle around the Z axis after the Y rotation.
+    lambda_ : ParameterizedAngle
+        Rotation angle around the Z axis before the Y rotation.
+    gamma : ParameterizedAngle
+        Global phase angle.
     """
 
     control: int
@@ -610,7 +639,7 @@ def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
 
 @dataclass(repr=False)
 class CP(_KindChecker, ControlledRotationInstruction):
-    r"""Controlled-P rotation circuit instruction.
+    r"""Controlled-P circuit instruction.
 
     The :math:`CP(\theta)` gate applies the matrix
 
@@ -622,6 +651,19 @@ class CP(_KindChecker, ControlledRotationInstruction):
         0 & 0 & 1 & 0\\
         0 & 0 & 0 & \mathrm e^{\mathrm i \theta}
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    angle : ParameterizedAngle
+        Phase angle.
     """
 
     kind: ClassVar[Literal[InstructionKind.CP]] = field(default=InstructionKind.CP, init=False)
@@ -641,6 +683,19 @@ class CRX(_KindChecker, ControlledRotationInstruction):
         0 & 0 & \cos \frac \theta 2 & -\mathrm i \sin \frac \theta 2\\
         0 & 0 & -\mathrm i \sin \frac \theta 2 & \cos \frac \theta 2
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    angle : ParameterizedAngle
+        Rotation angle.
     """
 
     kind: ClassVar[Literal[InstructionKind.CRX]] = field(default=InstructionKind.CRX, init=False)
@@ -658,6 +713,19 @@ class CRY(_KindChecker, ControlledRotationInstruction):
         0 & 0 & \cos \frac \theta 2 & - \sin \frac \theta 2\\
         0 & 0 & \sin \frac \theta 2 & \cos \frac \theta 2
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    angle : ParameterizedAngle
+        Rotation angle.
     """
 
     kind: ClassVar[Literal[InstructionKind.CRY]] = field(default=InstructionKind.CRY, init=False)
@@ -675,6 +743,19 @@ class CRZ(_KindChecker, ControlledRotationInstruction):
         0 & 0 & \mathrm e^{-\mathrm i \frac \theta 2} & 0\\
         0 & 0 & 0 & \mathrm e^{\mathrm i \frac \theta 2}
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    angle : ParameterizedAngle
+        Rotation angle.
     """
 
     kind: ClassVar[Literal[InstructionKind.CRZ]] = field(default=InstructionKind.CRZ, init=False)
@@ -695,6 +776,10 @@ class CJ(_KindChecker, ControlledRotationInstruction):
         0 & 0 & \frac 1 {\sqrt 2} & - \frac 1 {\sqrt 2} \mathrm e^{\mathrm i \alpha}
       \end{matrix}\right]
 
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
     Following Lemmas 4.3 and 5.1 of Barenco et al. (1995), we define:
 
     .. math::
@@ -716,6 +801,15 @@ class CJ(_KindChecker, ControlledRotationInstruction):
 
       CJ(\alpha) = (P(\delta) \otimes I) \, (I \otimes A) \, CX \, (I \otimes B) \, CX \, (I \otimes C)
 
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    angle : ParameterizedAngle
+        Rotation angle.
+
     References
     ----------
     Barenco, A., Bennett, C. H., Cleve, R., DiVincenzo, D. P., Margolus, N., Shor, P., Sleator, T., Smolin, J. A., & Weinfurter, H. (1995).
@@ -728,7 +822,13 @@ class CJ(_KindChecker, ControlledRotationInstruction):
 
 @dataclass(repr=False)
 class GPHASE(_KindChecker, BaseInstruction):
-    """GPHASE circuit instruction."""
+    """GPHASE circuit instruction.
+
+    Attributes
+    ----------
+    angle : ParameterizedAngle
+        Phase angle.
+    """
 
     angle: ParameterizedAngle = field(metadata={"repr": repr_angle})
     kind: ClassVar[Literal[InstructionKind.GPHASE]] = field(default=InstructionKind.GPHASE, init=False)
diff --git a/graphix/ops.py b/graphix/ops.py
index 9d4bfd864..5216a77ce 100644
--- a/graphix/ops.py
+++ b/graphix/ops.py
@@ -120,19 +120,21 @@ class Ops:
 
     @overload
     @staticmethod
-    def _cast_array(array: Iterable[Iterable[complex]], theta: Angle) -> npt.NDArray[np.complex128]: ...
+    def _cast_array(
+        array: Iterable[Iterable[complex]], parameters: tuple[Angle, ...]
+    ) -> npt.NDArray[np.complex128]: ...
 
     @overload
     @staticmethod
     def _cast_array(
-        array: Iterable[Iterable[ExpressionOrComplex]], theta: ParameterizedAngle
+        array: Iterable[Iterable[ExpressionOrComplex]], parameters: tuple[ParameterizedAngle, ...]
     ) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]: ...
 
     @staticmethod
     def _cast_array(
-        array: Iterable[Iterable[ExpressionOrComplex]], theta: ParameterizedAngle
+        array: Iterable[Iterable[ExpressionOrComplex]], parameters: tuple[ParameterizedAngle, ...]
     ) -> npt.NDArray[np.complex128] | npt.NDArray[np.object_]:
-        if isinstance(theta, Expression):
+        if all(isinstance(parameter, Expression) for parameter in parameters):
             return np.asarray(array, dtype=np.object_)
         return np.asarray(array, dtype=np.complex128)
 
@@ -159,7 +161,7 @@ def p(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.
         -------
         operator : 2*2 np.asarray
         """
-        return Ops._cast_array([[1, 0], [0, exp(1j * angle_to_rad(theta))]], theta)
+        return Ops._cast_array([[1, 0], [0, exp(1j * angle_to_rad(theta))]], (theta,))
 
     @overload
     @staticmethod
@@ -183,7 +185,7 @@ def rx(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np
         operator : 2*2 np.asarray
         """
         cos, sin = cos_sin(angle_to_rad(theta) / 2)
-        return Ops._cast_array([[cos, -1j * sin], [-1j * sin, cos]], theta)
+        return Ops._cast_array([[cos, -1j * sin], [-1j * sin, cos]], (theta,))
 
     @overload
     @staticmethod
@@ -207,7 +209,7 @@ def ry(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np
         operator : 2*2 np.asarray
         """
         cos, sin = cos_sin(angle_to_rad(theta) / 2)
-        return Ops._cast_array([[cos, -sin], [sin, cos]], theta)
+        return Ops._cast_array([[cos, -sin], [sin, cos]], (theta,))
 
     @overload
     @staticmethod
@@ -230,7 +232,9 @@ def rz(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np
         -------
         operator : 2*2 np.asarray
         """
-        return Ops._cast_array([[exp(-1j * angle_to_rad(theta) / 2), 0], [0, exp(1j * angle_to_rad(theta) / 2)]], theta)
+        return Ops._cast_array(
+            [[exp(-1j * angle_to_rad(theta) / 2), 0], [0, exp(1j * angle_to_rad(theta) / 2)]], (theta,)
+        )
 
     @staticmethod
     def u(
@@ -256,7 +260,7 @@ def u(
         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,
+            (theta, phi, lambda_),
         )
 
     @staticmethod
@@ -424,7 +428,7 @@ def j(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[np.
                 [1 / np.sqrt(2), (1 / np.sqrt(2)) * exp(1j * angle_to_rad(theta))],
                 [1 / np.sqrt(2), (-1 / np.sqrt(2)) * exp(1j * angle_to_rad(theta))],
             ],
-            theta,
+            (theta,),
         )
 
     @overload
@@ -453,7 +457,7 @@ def rzz(theta: ParameterizedAngle) -> npt.NDArray[np.complex128] | npt.NDArray[n
         -------
         operator : 4*4 np.asarray
         """
-        return Ops._cast_array(Ops.CNOT @ np.kron(Ops.I, Ops.rz(theta)) @ Ops.CNOT, theta)
+        return Ops._cast_array(Ops.CNOT @ np.kron(Ops.I, Ops.rz(theta)) @ Ops.CNOT, (theta,))
 
     @staticmethod
     def build_tensor_pauli_ops(n_qubits: int) -> npt.NDArray[np.complex128]:

From 24c5b5e75aa49abbc3640025e4cc6e5ff8a73ee0 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Wed, 2 Sep 2026 06:18:21 +0200
Subject: [PATCH 15/18] Add explanation in test_visit_qubit

---
 tests/test_instruction.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/tests/test_instruction.py b/tests/test_instruction.py
index 32b18f2bc..0354b8813 100644
--- a/tests/test_instruction.py
+++ b/tests/test_instruction.py
@@ -105,6 +105,9 @@ def test_visit_qubit(fx_rng: Generator, test_case: InstructionTestCase) -> None:
     visitor = VisitQubit()
     instr_visited = instr.visit(visitor, copy=True)
     assert instr == instr_copy
+    # instr_visited differs from the original instr_copy iff a qubit
+    # attribute has been modified by the visitor, and GPHASE gate has
+    # no qubit attribute.
     if test_case.name != "GPHASE":
         assert instr_visited != instr_copy
     instr.visit(visitor, copy=False)

From 986e790c81fa056d8a686413acb12e082b3d7226 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 7 Sep 2026 01:01:30 +0200
Subject: [PATCH 16/18] Better documentation

---
 graphix/instruction.py | 78 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 77 insertions(+), 1 deletion(-)

diff --git a/graphix/instruction.py b/graphix/instruction.py
index e502b6527..e40563059 100644
--- a/graphix/instruction.py
+++ b/graphix/instruction.py
@@ -150,6 +150,17 @@ class CCX(_KindChecker, BaseInstruction):
       0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
       0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``controls[0]``, ``controls[1]``, ``target``.
+
+    Attributes
+    ----------
+    controls : tuple[int, int]
+        Index of the control qubits.
+    target : int
+        Index of the target qubit.
     """
 
     target: int
@@ -183,7 +194,20 @@ class RZZ(_KindChecker, BaseInstruction):
         0 & 0 & 0 & \mathrm e^{-\mathrm i \frac \theta 2}
       \end{matrix}\right]
 
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
     We have :math:`\mathrm{RZZ}(\theta) = \mathrm{CNOT} (I \otimes \mathrm{RZ}(\theta)) \mathrm{CNOT}`.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    angle : ParameterizedAngle
+        Rotation angle.
     """
 
     target: int
@@ -224,7 +248,30 @@ def visit(self, visitor: InstructionVisitor, *, copy: bool = False) -> Self:
 
 @dataclass(repr=False)
 class CY(_KindChecker, ControlledSingleTargetInstruction):
-    """CY circuit instruction."""
+    r"""CY circuit instruction.
+
+    The CY gate applies the matrix
+
+    .. math::
+
+      \left[\begin{matrix}
+        1 & 0 & 0 & 0\\
+        0 & 1 & 0 & 0\\
+        0 & 0 & 0 & -\mathrm i\\
+        0 & 0 & \mathrm i & 0
+      \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
+    """
 
     kind: ClassVar[Literal[InstructionKind.CY]] = field(default=InstructionKind.CY, init=False)
 
@@ -243,6 +290,17 @@ class CNOT(_KindChecker, ControlledSingleTargetInstruction):
         0 & 0 & 0 & 1\\
         0 & 0 & 1 & 0
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``control``, ``target``.
+
+    Attributes
+    ----------
+    control : int
+        Index of the control qubit.
+    target : int
+        Index of the target qubit.
     """
 
     kind: ClassVar[Literal[InstructionKind.CNOT]] = field(default=InstructionKind.CNOT, init=False)
@@ -264,6 +322,15 @@ class CZ(_KindChecker, BaseInstruction):
         0 & 0 & 1 & 0\\
         0 & 0 & 0 & -1
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``targets[0]``, ``targets[1]``.
+
+    Attributes
+    ----------
+    targets : tuple[int, int]
+        Index of the target qubits.
     """
 
     targets: tuple[int, int]
@@ -293,6 +360,15 @@ class SWAP(_KindChecker, BaseInstruction):
         0 & 1 & 0 & 0\\
         0 & 0 & 0 & 1
       \end{matrix}\right]
+
+    in the computational basis. The basis states use big-endian
+    ordering, with the most significant qubit first. The qubits are
+    numbered in the order ``targets[0]``, ``targets[1]``.
+
+    Attributes
+    ----------
+    targets : tuple[int, int]
+        Index of the target qubits.
     """
 
     targets: tuple[int, int]

From 008e03fd4e0c660e16eff73463d80a05189e4d67 Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Mon, 7 Sep 2026 01:17:27 +0200
Subject: [PATCH 17/18] Fix documentation

---
 graphix/instruction.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/graphix/instruction.py b/graphix/instruction.py
index e40563059..641bd5880 100644
--- a/graphix/instruction.py
+++ b/graphix/instruction.py
@@ -26,7 +26,7 @@ def repr_angle(angle: ParameterizedAngle) -> str:
     """
     Return the representation string of an angle in radians.
 
-    This is used for pretty-printing instructions with `angle` parameters.
+    This is used for pretty-printing instructions with ``angle`` parameters.
     Delegates to :func:`pretty_print.angle_to_str`.
     """
     # Non-float-supporting objects are returned as-is
@@ -991,7 +991,6 @@ class CJ(_KindChecker, ControlledRotationInstruction):
     Following Lemmas 4.3 and 5.1 of Barenco et al. (1995), we define:
 
     .. math::
-
       \begin{aligned}
         A &= R_Y\left(\frac \pi 4\right),\\
         B &= R_Y\left(- \frac \pi 4\right) R_Z(- \delta),\\
@@ -1000,8 +999,7 @@ class CJ(_KindChecker, ControlledRotationInstruction):
       \end{aligned}
 
     These operators satisfy :math:`ABC = I` and
-    :math:`AXBXC = \mathrm e^{-\mathrm i \delta} J(\alpha)` with
-    :math:``.
+    :math:`AXBXC = \mathrm e^{-\mathrm i \delta} J(\alpha)`.
 
     Consequently, :math:`CJ(\alpha)` can be decomposed as:
 

From 368210902deffd3c9c2fc5726df17de4bb696bbb Mon Sep 17 00:00:00 2001
From: Thierry Martinez <thierry.martinez@inria.fr>
Date: Sun, 13 Sep 2026 17:15:46 +0200
Subject: [PATCH 18/18] Use TeamGraphix/graphix-qasm-parser#15

---
 .github/qasm-parser-requirements.txt | 2 +-
 noxfile.py                           | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/qasm-parser-requirements.txt b/.github/qasm-parser-requirements.txt
index 5c922358d..bfb39feb4 100644
--- a/.github/qasm-parser-requirements.txt
+++ b/.github/qasm-parser-requirements.txt
@@ -1 +1 @@
-graphix-qasm-parser @ git+https://github.com/TeamGraphix/graphix-qasm-parser@refs/pull/16/head
+graphix-qasm-parser @ git+https://github.com/TeamGraphix/graphix-qasm-parser@refs/pull/15/head
diff --git a/noxfile.py b/noxfile.py
index 32d853def..b47bd604c 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -97,7 +97,7 @@ class ReverseDependency:
     [
         ReverseDependency("https://github.com/thierry-martinez/graphix-symbolic", branch="in-place_methods"),
         ReverseDependency("https://github.com/thierry-martinez/graphix-stim-backend", branch="rename-simulate"),
-        ReverseDependency("https://github.com/TeamGraphix/graphix-qasm-parser", branch="refs/pull/16/head"),
+        ReverseDependency("https://github.com/TeamGraphix/graphix-qasm-parser", branch="refs/pull/15/head"),
         ReverseDependency(
             "https://github.com/thierry-martinez/graphix-ibmq", doctest_modules=False, branch="rename-simulate"
         ),