From a74dbe96bbec154f9c443b53ceb9e2d5fd5f1364 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 15:35:57 -0400
Subject: [PATCH 01/14] Add changelog entry for Pauli-product rotation gates
 (#142)

---
 CHANGELOG.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 125d4f9d..c5a6157e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+
+`R_XX`, `R_YY`, `R_ZZ`, and `R_PAULI` parametric Pauli-product rotation gates. `R_XX(alpha) q0 q1` applies `exp(-i * alpha * pi/2 * XX)` (and likewise for `YY`/`ZZ`), while `R_PAULI(alpha) X0*Y1*Z2` applies the rotation for an arbitrary Pauli product (up to 64 qubits per instruction). The `alpha` parameter is in half-turns, matching the existing `R_X`/`R_Y`/`R_Z` convention. Duplicate target qubits are rejected at parse time. (#142)
+
 - Fast path in the detector sampler for components whose output is deterministically given by a single error variable. These components now skip the JAX compilation and autoregressive sampling pipeline, significantly speeding up detector sampling for surface-code circuits at low physical error rates.
 
 

From 1a7da1490fc2cab53b1e9cbc8496f0a039f4314f Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 17:12:19 -0400
Subject: [PATCH 02/14] Add r_pauli Pauli-product rotation primitive (#142)

---
 src/tsim/core/instructions.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/tsim/core/instructions.py b/src/tsim/core/instructions.py
index 0e62d734..98236034 100644
--- a/src/tsim/core/instructions.py
+++ b/src/tsim/core/instructions.py
@@ -1001,6 +1001,35 @@ def tpp(
     _pauli_product_phase(b, paulis, t, t_dag, dagger)
 
 
+def r_pauli(
+    b: GraphRepresentation,
+    paulis: list[tuple[Literal["X", "Y", "Z"], int]],
+    theta: Fraction,
+    dagger: bool = False,
+) -> None:
+    """Apply exp(-i theta pi/2 P) for a Pauli product P, with theta in half-turns.
+
+    Generalizes ``SPP`` (theta=1/2) and ``TPP`` (theta=1/4) to an arbitrary angle.
+    For a single qubit, ``r_pauli`` reduces to the named single-qubit rotation in
+    that Pauli's basis (e.g. ``R_PAULI Z0`` is ``R_Z``). If `dagger` is True, apply
+    exp(+i theta pi/2 P) instead.
+
+    Args:
+        b: The graph representation to modify.
+        paulis: List of (pauli_type, qubit) pairs defining the Pauli product P.
+        theta: Rotation angle in half-turns (units of pi).
+        dagger: If True, negate the rotation angle.
+
+    """
+    _pauli_product_phase(
+        b,
+        paulis,
+        lambda b, qubit: r_z(b, qubit, theta),
+        lambda b, qubit: r_z(b, qubit, -theta),
+        dagger,
+    )
+
+
 def mpad(b: GraphRepresentation, value: int, p: float = 0) -> None:
     """Pad measurement record with a fixed bit value.
 

From 7dc2108e23fa60f369d15764587e768d8e7ad302 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 17:21:54 -0400
Subject: [PATCH 03/14] Parse SPP[R_PAULI(...)] tag into r_pauli rotations,
 enforce 64-qubit limit (#142)

---
 src/tsim/core/parse.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/src/tsim/core/parse.py b/src/tsim/core/parse.py
index 0e74c712..6dd4251b 100644
--- a/src/tsim/core/parse.py
+++ b/src/tsim/core/parse.py
@@ -16,6 +16,7 @@
     mpad,
     mpp,
     observable_include,
+    r_pauli,
     r_x,
     r_y,
     r_z,
@@ -29,9 +30,12 @@
     "R_X": frozenset({"theta"}),
     "R_Y": frozenset({"theta"}),
     "R_Z": frozenset({"theta"}),
+    "R_PAULI": frozenset({"theta"}),
     "U3": frozenset({"theta", "phi", "lambda"}),
 }
 
+R_PAULI_MAX_QUBITS = 64
+
 
 def parse_parametric_tag(
     instruction: stim.CircuitInstruction,
@@ -230,6 +234,22 @@ def parse_stim_circuit(
             for paulis, invert in _iter_pauli_products(instruction):
                 tpp(b, paulis, dagger=is_dag ^ invert)
             continue
+        if name in ("SPP", "SPP_DAG") and instruction.tag:
+            parsed = parse_parametric_tag(instruction)
+            if parsed is not None and parsed[0] == "R_PAULI":
+                params = parsed[1]
+                n_qubits = len(
+                    {t.value for t in instruction.targets_copy() if not t.is_combiner}
+                )
+                if n_qubits > R_PAULI_MAX_QUBITS:
+                    raise ValueError(
+                        f"R_PAULI supports at most {R_PAULI_MAX_QUBITS} qubits per "
+                        f"instruction, got {n_qubits}."
+                    )
+                is_dag = name == "SPP_DAG"
+                for paulis, invert in _iter_pauli_products(instruction):
+                    r_pauli(b, paulis, params["theta"], dagger=is_dag ^ invert)
+                continue
         if name in ("SPP", "SPP_DAG"):
             is_dag = name == "SPP_DAG"
             for paulis, invert in _iter_pauli_products(instruction):

From d9d01e38ccc65f66a04b714a85ece2431c8404ca Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 17:28:59 -0400
Subject: [PATCH 04/14] Add R_XX/R_YY/R_ZZ/R_PAULI dispatch to Circuit.append
 (#142)

---
 src/tsim/circuit.py | 55 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/src/tsim/circuit.py b/src/tsim/circuit.py
index 09210567..ca57c5a4 100644
--- a/src/tsim/circuit.py
+++ b/src/tsim/circuit.py
@@ -30,6 +30,48 @@
     stim_to_shorthand,
 )
 
+_PAULI_TARGET = {"X": stim.target_x, "Y": stim.target_y, "Z": stim.target_z}
+
+
+def _single_angle(name: str, arg: float | Iterable[float] | None) -> float:
+    """Extract the single rotation angle required by a parametric gate."""
+    if arg is None:
+        raise ValueError(f"For {name} gates, an angle must be provided.")
+    args = list(arg) if isinstance(arg, Iterable) else [arg]
+    if len(args) != 1:
+        raise ValueError(f"For {name} gates, a single angle must be provided.")
+    return args[0]
+
+
+def _two_distinct_qubits(
+    name: str,
+    targets: int | stim.GateTarget | stim.PauliString | Iterable,
+) -> tuple[int, int]:
+    """Validate and unpack the two distinct qubit targets of an R_XX/R_YY/R_ZZ gate."""
+    qubits = list(targets) if isinstance(targets, Iterable) else [targets]
+    if len(qubits) != 2:
+        raise ValueError(f"For {name} gates, exactly two qubit targets are required.")
+    q0, q1 = qubits
+    if not isinstance(q0, int) or not isinstance(q1, int):
+        raise ValueError(f"For {name} gates, both targets must be qubit indices.")
+    if q0 == q1:
+        raise ValueError(
+            f"For {name} gates, the two target qubits must be distinct, got {q0} {q1}."
+        )
+    return q0, q1
+
+
+def _pauli_product_targets(
+    paulis: list[tuple[str, int]],
+) -> list[stim.GateTarget]:
+    """Build combiner-joined Pauli targets (e.g. ``X0*X1``) for an SPP instruction."""
+    out: list[stim.GateTarget] = []
+    for pauli, qubit in paulis:
+        if out:
+            out.append(stim.target_combiner())
+        out.append(_PAULI_TARGET[pauli](qubit))
+    return out
+
 
 class Circuit:
     """Quantum circuit as a thin wrapper around stim.Circuit.
@@ -183,6 +225,19 @@ def append(
                 tag = f"U3(theta={theta}*pi, phi={phi}*pi, lambda={lam}*pi)"
                 name = "I"
                 arg = None
+            elif name in ("R_XX", "R_YY", "R_ZZ"):
+                alpha = _single_angle(name, arg)
+                pauli = name[2]
+                q0, q1 = _two_distinct_qubits(name, targets)
+                targets = _pauli_product_targets([(pauli, q0), (pauli, q1)])
+                tag = f"R_PAULI(theta={alpha}*pi)"
+                name = "SPP"
+                arg = None
+            elif name == "R_PAULI":
+                alpha = _single_angle(name, arg)
+                tag = f"R_PAULI(theta={alpha}*pi)"
+                name = "SPP"
+                arg = None
 
             self._stim_circ.append(name=name, targets=targets, arg=arg, tag=tag)  # type: ignore
         else:

From aceb143fc5be21b7b9e4bb813949f3da84b04072 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 17:37:34 -0400
Subject: [PATCH 05/14] Support R_XX/R_YY/R_ZZ/R_PAULI shorthand round-trip
 (#142)

---
 src/tsim/utils/program_text.py | 58 ++++++++++++++++++++++++++++++++--
 1 file changed, 56 insertions(+), 2 deletions(-)

diff --git a/src/tsim/utils/program_text.py b/src/tsim/utils/program_text.py
index ed8f6f41..47cecec1 100644
--- a/src/tsim/utils/program_text.py
+++ b/src/tsim/utils/program_text.py
@@ -5,9 +5,11 @@
 # Matches valid numeric literals including scientific notation (e.g. 0.5, 4e-4, 1.2e3)
 FLOAT_RE = r"[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?"
 
-_TSIM_GATES = {"R_X", "R_Y", "R_Z", "U3"}
+_TSIM_GATES = {"R_X", "R_Y", "R_Z", "R_XX", "R_YY", "R_ZZ", "R_PAULI", "U3"}
 _GATE_NOT_FOUND_RE = re.compile(r"Gate not found: '(\w+)'")
-_GATE_USAGE_RE = re.compile(r"(?<!\[)\b(R_[A-Z]\([^)]*\)|R_[XYZ]\b|U3\([^)]*\)|U3\b)")
+_GATE_USAGE_RE = re.compile(
+    r"(?<!\[)\b(R_PAULI\([^)]*\)|R_[XYZ]{1,2}\([^)]*\)|R_[XYZ]\b|U3\([^)]*\)|U3\b)"
+)
 
 
 def enriched_stim_error(exc: ValueError, converted_text: str) -> ValueError:
@@ -38,6 +40,10 @@ def shorthand_to_stim(text: str) -> str:
         R_Z(0.3) 0          → I[R_Z(theta=0.3*pi)] 0
         R_X(0.25) 0         → I[R_X(theta=0.25*pi)] 0
         R_Y(-0.5) 0         → I[R_Y(theta=-0.5*pi)] 0
+        R_XX(0.5) 0 1       → SPP[R_PAULI(theta=0.5*pi)] X0*X1
+        R_YY(0.5) 0 1       → SPP[R_PAULI(theta=0.5*pi)] Y0*Y1
+        R_ZZ(0.5) 0 1       → SPP[R_PAULI(theta=0.5*pi)] Z0*Z1
+        R_PAULI(0.3) X0*Y1  → SPP[R_PAULI(theta=0.3*pi)] X0*Y1
         U3(0.3, 0.24, 0.49) 0 → I[U3(theta=0.3*pi, phi=0.24*pi, lambda=0.49*pi)] 0
 
     """
@@ -48,6 +54,30 @@ def shorthand_to_stim(text: str) -> str:
     text = re.sub(r"(?<!\[)\bT_DAG\b(?!\[)", "S_DAG[T]", text)
     text = re.sub(r"(?<!\[)\bT\b(?!\[)", "S[T]", text)
 
+    def replace_pauli_pair(m: re.Match) -> str:
+        pauli = m.group(1)
+        alpha = float(m.group(2))
+        q0, q1 = m.group(3), m.group(4)
+        if q0 == q1:
+            raise ValueError(
+                f"R_{pauli}{pauli} target qubits must be distinct, got {q0} {q1}."
+            )
+        return f"SPP[R_PAULI(theta={alpha}*pi)] {pauli}{q0}*{pauli}{q1}"
+
+    text = re.sub(
+        rf"\bR_([XYZ])\1\(({FLOAT_RE})\)\s+(\d+)\s+(\d+)", replace_pauli_pair, text
+    )
+
+    def replace_pauli(m: re.Match) -> str:
+        alpha = float(m.group(1))
+        return f"SPP[R_PAULI(theta={alpha}*pi)] {m.group(2)}"
+
+    text = re.sub(
+        rf"\bR_PAULI\(({FLOAT_RE})\)\s+((?:[XYZ]\d+)(?:\*[XYZ]\d+)*)",
+        replace_pauli,
+        text,
+    )
+
     def replace_rotation(m: re.Match) -> str:
         axis = m.group(1)
         return f"I[R_{axis}(theta={float(m.group(2))}*pi)]"
@@ -86,6 +116,8 @@ def stim_to_shorthand(text: str) -> str:
     Rewrites:
     - I[U3(theta=θ*pi, phi=φ*pi, lambda=λ*pi)] → U3(θ, φ, λ)
     - I[R_X(theta=α*pi)] / I[R_Y(...)] / I[R_Z(...)] → R_X(α) / R_Y(α) / R_Z(α)
+    - SPP[R_PAULI(theta=α*pi)] P0*P1 → R_PP(α) 0 1 (P ∈ {X, Y, Z})
+    - SPP[R_PAULI(theta=α*pi)] X0*Y1 → R_PAULI(α) X0*Y1
     - SPP[T] → TPP
     - SPP_DAG[T] → TPP_DAG
     - S[T] → T
@@ -103,6 +135,28 @@ def replace_u3(m: re.Match) -> str:
         text,
     )
 
+    # SPP[R_PAULI(...)] with a same-axis two-qubit product → R_XX/R_YY/R_ZZ shorthand.
+    # Must precede the general R_PAULI rewrite below.
+    def replace_pauli_pair(m: re.Match) -> str:
+        alpha, pauli, q0, q1 = m.group(1), m.group(2), m.group(3), m.group(4)
+        return f"R_{pauli}{pauli}({alpha}) {q0} {q1}"
+
+    text = re.sub(
+        rf"\bSPP\[R_PAULI\(theta=({FLOAT_RE})\*pi\)\] ([XYZ])(\d+)\*\2(\d+)",
+        replace_pauli_pair,
+        text,
+    )
+
+    def replace_pauli(m: re.Match) -> str:
+        alpha, product = m.group(1), m.group(2)
+        return f"R_PAULI({alpha}) {product}"
+
+    text = re.sub(
+        rf"\bSPP\[R_PAULI\(theta=({FLOAT_RE})\*pi\)\] ((?:[XYZ]\d+)(?:\*[XYZ]\d+)*)",
+        replace_pauli,
+        text,
+    )
+
     # Replace I[R_X(...)] / I[R_Y(...)] / I[R_Z(...)] with R_X(α) / R_Y(α) / R_Z(α)
     def replace_rotation(m: re.Match) -> str:
         axis = m.group(1)

From 6cc9b22e8089969084061653d4e42c57f3873691 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 17:45:17 -0400
Subject: [PATCH 06/14] Classify and expand Clifford-angle Pauli rotations
 (#142)

---
 src/tsim/utils/clifford.py | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/src/tsim/utils/clifford.py b/src/tsim/utils/clifford.py
index f27cc06f..9f78cece 100644
--- a/src/tsim/utils/clifford.py
+++ b/src/tsim/utils/clifford.py
@@ -119,6 +119,12 @@ def is_half_pi_multiple(phase: Fraction) -> bool:
         if instr.name in ["S", "S_DAG", "SPP", "SPP_DAG"] and instr.tag == "T":
             return False
 
+        if instr.name in ["SPP", "SPP_DAG"] and instr.tag:
+            result = parse_parametric_tag(instr)
+            if result is not None and not is_half_pi_multiple(result[1]["theta"]):
+                return False
+            continue
+
         if instr.name == "I" and instr.tag:
             result = parse_parametric_tag(instr)
             if result is None:
@@ -154,6 +160,8 @@ def expand_clifford_rotations(source: stim.Circuit) -> stim.Circuit:
                 )
             )
             continue
+        if _expand_clifford_spp(instr, out):
+            continue
         expansion = _try_clifford_expansion(instr)
         if expansion is not None:
             gates, targets = expansion
@@ -164,6 +172,33 @@ def expand_clifford_rotations(source: stim.Circuit) -> stim.Circuit:
     return out
 
 
+def _expand_clifford_spp(instr: stim.CircuitInstruction, out: stim.Circuit) -> bool:
+    """Expand a Clifford-angle ``SPP[R_PAULI(...)]`` into plain SPP gates.
+
+    Appends the equivalent tag-free instructions to ``out`` and returns True when
+    the instruction was an expandable Clifford-angle Pauli rotation; returns False
+    (appending nothing) otherwise, so the caller falls through to other handling.
+    """
+    if instr.name not in ("SPP", "SPP_DAG") or not instr.tag:
+        return False
+
+    parsed = parse_parametric_tag(instr)
+    if parsed is None or parsed[0] != "R_PAULI":
+        return False
+
+    idx = _to_half_pi_index(parsed[1]["theta"])
+    if idx is None:
+        return False
+
+    targets = instr.targets_copy()
+    if instr.name == "SPP_DAG":
+        idx = (4 - idx) % 4
+    repeats = {0: [], 1: ["SPP"], 2: ["SPP", "SPP"], 3: ["SPP_DAG"]}[idx]
+    for gate in repeats:
+        out.append(gate, targets, [])
+    return True
+
+
 def _try_clifford_expansion(
     instr: stim.CircuitInstruction,
 ) -> tuple[list[str], list[int]] | None:

From a1c1d3f5cf839200877893f431ea6cc248206079 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 19:59:22 -0400
Subject: [PATCH 07/14] Add analytic R_XX/R_YY/R_ZZ rotation matrices for tests
 (#142)

---
 test/helpers/gate_matrices.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/test/helpers/gate_matrices.py b/test/helpers/gate_matrices.py
index f669f196..9182fa0d 100644
--- a/test/helpers/gate_matrices.py
+++ b/test/helpers/gate_matrices.py
@@ -114,6 +114,19 @@
     "R_Z": lambda frac: np.array(
         [[np.exp(-1j * np.pi / 2 * frac), 0], [0, np.exp(1j * np.pi / 2 * frac)]]
     ),
+    # Two-qubit Pauli rotations: exp(-i frac pi/2 PP) = cos(frac pi/2) I - i sin(frac pi/2) PP
+    "R_XX": lambda frac: np.cos(frac * np.pi / 2) * np.eye(4)
+    - 1j
+    * np.sin(frac * np.pi / 2)
+    * np.kron(SINGLE_QUBIT_GATE_MATRICES["X"], SINGLE_QUBIT_GATE_MATRICES["X"]),
+    "R_YY": lambda frac: np.cos(frac * np.pi / 2) * np.eye(4)
+    - 1j
+    * np.sin(frac * np.pi / 2)
+    * np.kron(SINGLE_QUBIT_GATE_MATRICES["Y"], SINGLE_QUBIT_GATE_MATRICES["Y"]),
+    "R_ZZ": lambda frac: np.cos(frac * np.pi / 2) * np.eye(4)
+    - 1j
+    * np.sin(frac * np.pi / 2)
+    * np.kron(SINGLE_QUBIT_GATE_MATRICES["Z"], SINGLE_QUBIT_GATE_MATRICES["Z"]),
     "U3": lambda frac_theta, frac_phi, frac_lambda: np.array(
         [
             [

From f799131c651c29116587d3973767b8eb63786d04 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 20:28:18 -0400
Subject: [PATCH 08/14] Test Pauli-rotation arbitrary-angle unitaries (#142)

---
 test/integration/test_gate_unitaries.py | 30 +++++++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/test/integration/test_gate_unitaries.py b/test/integration/test_gate_unitaries.py
index 2e14fe54..4274331a 100644
--- a/test/integration/test_gate_unitaries.py
+++ b/test/integration/test_gate_unitaries.py
@@ -336,3 +336,33 @@ def test_heralded_erase_bell_state():
     mat = get_heralded_matrix(sampler, herald_value=1, batch_size=10000)
     expected = np.full((2, 2), 0.5)
     assert np.allclose(mat, expected, atol=0.05)
+
+
+@pytest.mark.parametrize("gate", ["R_XX", "R_YY", "R_ZZ"])
+def test_pauli_rotation_arbitrary_angle(gate: str):
+    c = Circuit(f"""
+        R 0 1 2 3
+        H 0 1
+        CNOT 0 2 1 3
+        {gate}(0.345) 0 1
+        M 0 1 2 3
+        """)
+    sampler = CompiledStateProbs(c)
+    mat = get_matrix(sampler)
+    expected = np.abs(ROT_GATE_MATRICES[gate](0.345)) ** 2
+    assert np.allclose(mat, expected)
+
+
+def test_r_pauli_arbitrary_angle():
+    c = Circuit("""
+        R 0 1 2 3
+        H 0 1
+        CNOT 0 2 1 3
+        R_PAULI(0.345) X0*Y1
+        M 0 1 2 3
+        """)
+    sampler = CompiledStateProbs(c)
+    mat = get_matrix(sampler)
+    product = np.kron(SINGLE_QUBIT_GATE_MATRICES["X"], SINGLE_QUBIT_GATE_MATRICES["Y"])
+    expected = np.abs(ROT_GATE_MATRICES["R_PAULI"](0.345, product)) ** 2
+    assert np.allclose(mat, expected)

From 19c34d36edd7b9756ba2ce98440ebd50a570b869 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 20:30:59 -0400
Subject: [PATCH 09/14] Add analytic R_PAULI rotation matrix for tests (#142)

---
 test/helpers/gate_matrices.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/helpers/gate_matrices.py b/test/helpers/gate_matrices.py
index 9182fa0d..eaca39b3 100644
--- a/test/helpers/gate_matrices.py
+++ b/test/helpers/gate_matrices.py
@@ -127,6 +127,10 @@
     - 1j
     * np.sin(frac * np.pi / 2)
     * np.kron(SINGLE_QUBIT_GATE_MATRICES["Z"], SINGLE_QUBIT_GATE_MATRICES["Z"]),
+    # General Pauli-product rotation: exp(-i frac pi/2 P) for a Pauli product matrix P
+    "R_PAULI": lambda frac, pauli_product: np.cos(frac * np.pi / 2)
+    * np.eye(pauli_product.shape[0])
+    - 1j * np.sin(frac * np.pi / 2) * pauli_product,
     "U3": lambda frac_theta, frac_phi, frac_lambda: np.array(
         [
             [

From efceda5a3f1f5b8585a7157740281c0387a3d3ed Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 20:32:22 -0400
Subject: [PATCH 10/14] Test Pauli-rotation Clifford stim-parity and inverse
 round-trip (#142)

---
 test/unit/test_circuit.py | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/test/unit/test_circuit.py b/test/unit/test_circuit.py
index 061e3563..4d5dd829 100644
--- a/test/unit/test_circuit.py
+++ b/test/unit/test_circuit.py
@@ -183,6 +183,23 @@ def test_two_qubit_gate(stim_gate: str):
     assert unitaries_equal_up_to_global_phase(c.to_matrix(), stim_c_matrix)
 
 
+@pytest.mark.parametrize("gate", ["R_XX", "R_YY", "R_ZZ"])
+@pytest.mark.parametrize("alpha", [0.0, 1.0, 2.0, 3.0])
+def test_pauli_rotation_clifford_matches_stim(gate: str, alpha: float):
+    """At Clifford angles, R_XX/R_YY/R_ZZ match stim's reference Clifford gate.
+
+    R_PP(alpha) = exp(-i alpha pi/2 PP) is Clifford for integer alpha: identity for
+    even alpha, and the Pauli pair PP for odd alpha (both up to global phase).
+    """
+    pauli = gate[2]
+    stim_program = "I 0\nI 1" if alpha % 2 == 0 else f"{pauli} 0\n{pauli} 1"
+    c = Circuit(f"{gate}({alpha}) 0 1")
+    stim_matrix = (
+        stim.Circuit(stim_program).to_tableau().to_unitary_matrix(endian="big")
+    )
+    assert unitaries_equal_up_to_global_phase(c.to_matrix(), stim_matrix)
+
+
 def test_num_measurements():
     c = Circuit()
     assert c.num_measurements == 0
@@ -754,6 +771,20 @@ def test_inverse_tpp_dag():
     assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
 
 
+def test_inverse_r_xx():
+    c = Circuit("R_XX(0.345) 0 1")
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
+def test_inverse_r_pauli():
+    c = Circuit("R_PAULI(0.345) X0*Y1*Z2")
+    c_inv = c.inverse()
+    combined = (c + c_inv).to_matrix()
+    assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
+
+
 def test_inverse_mixed_circuit():
     c = Circuit("H 0\nT 0\nR_Z(0.22) 0\nCNOT 0 1")
     c_inv = c.inverse()

From ba4b1fa3930148c2d5d76f01bacae33b99937180 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 20:36:06 -0400
Subject: [PATCH 11/14] Document Pauli-product rotation gates  (#142)

---
 CHANGELOG.md  |  4 +---
 README.md     | 13 ++++++++++++-
 docs/index.md |  2 +-
 3 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c5a6157e..bd07380a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,9 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
-
-`R_XX`, `R_YY`, `R_ZZ`, and `R_PAULI` parametric Pauli-product rotation gates. `R_XX(alpha) q0 q1` applies `exp(-i * alpha * pi/2 * XX)` (and likewise for `YY`/`ZZ`), while `R_PAULI(alpha) X0*Y1*Z2` applies the rotation for an arbitrary Pauli product (up to 64 qubits per instruction). The `alpha` parameter is in half-turns, matching the existing `R_X`/`R_Y`/`R_Z` convention. Duplicate target qubits are rejected at parse time. (#142)
-
+- `R_XX`, `R_YY`, `R_ZZ`, and `R_PAULI` parametric Pauli-product rotation gates. `R_XX(alpha) q0 q1` applies `exp(-i * alpha * pi/2 * XX)` (and likewise for `YY`/`ZZ`), while `R_PAULI(alpha) X0*Y1*Z2` applies the rotation for an arbitrary Pauli product (up to 64 qubits per instruction). The `alpha` parameter is in half-turns, matching the existing `R_X`/`R_Y`/`R_Z` convention. Duplicate target qubits are rejected at parse time. (#142)
 - Fast path in the detector sampler for components whose output is deterministically given by a single error variable. These components now skip the JAX compilation and autoregressive sampling pipeline, significantly speeding up detector sampling for surface-code circuits at low physical error rates.
 
 
diff --git a/README.md b/README.md
index e57cb47e..8568635c 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,7 @@ An introductory tutorial is available [here](https://tsim.mintlify.app/tutorials
 
 For many existing scripts, replacing `stim` with `tsim` should just work. Tsim mirrors the Stim API and currently supports all [Stim instructions](https://github.com/quantumlib/Stim/wiki/Stim-v1.13-Gate-Reference).
 
-Additionally, Tsim supports the instructions `T`, `T_DAG`, `R_Z`, `R_X`, `R_Y`, `U3`, `TPP`, and `TPP_DAG`.
+Additionally, Tsim supports the instructions `T`, `T_DAG`, `R_Z`, `R_X`, `R_Y`, `U3`, `TPP`, `TPP_DAG`, `R_XX`, `R_YY`, `R_ZZ`, and `R_PAULI`.
 ```python
 import tsim
 
@@ -144,6 +144,17 @@ TPP X0*Y1      # Apply exp(-i π/8 · X0⊗Y1) (up to global phase)
 TPP_DAG Z0     # Apply exp(+i π/8 · Z) = T_DAG (up to global phase)
 ```
 
+### Pauli Rotation Gates: `R_XX`, `R_YY`, `R_ZZ`, `R_PAULI`
+
+Parametric Pauli-product rotations apply exp(−i α·π/2 · P) for a Pauli product P, where α is specified as the parameter (in half-turns, matching `R_X`/`R_Y`/`R_Z`). `R_XX`, `R_YY` and `R_ZZ` are the two-qubit specializations; `R_PAULI` takes an arbitrary Pauli product (up to 64 qubits) in Stim's `X0*Y1*Z2` syntax. The two target qubits of `R_XX`/`R_YY`/`R_ZZ` must be distinct.
+
+```
+R_XX(0.5) 0 1          # Apply exp(-i π/4 · X0⊗X1)
+R_YY(0.25) 0 1         # Apply exp(-i π/8 · Y0⊗Y1)
+R_ZZ(0.5) 0 2          # Apply exp(-i π/4 · Z0⊗Z2)
+R_PAULI(0.3) X0*Y1*Z2  # Apply exp(-i 0.3·π/2 · X0⊗Y1⊗Z2)
+```
+
 ## Publications Using Tsim
 
 - [Simulating magic state cultivation with few Clifford terms](https://arxiv.org/abs/2509.08658) by Kwok Ho Wan and Zhenghao Zhong (2025).
diff --git a/docs/index.md b/docs/index.md
index f648655b..1db07e6e 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -9,7 +9,7 @@ A detailed description of Tsim is given in [arXiv:2604.01059](https://arxiv.org/
 
 An introductory tutorial is available [here](demos/encoding_demo.ipynb). For many existing scripts, replacing `stim` with `tsim` should just work. Tsim mirrors the Stim API and supports all [Stim instructions](https://github.com/quantumlib/Stim/wiki/Stim-v1.13-Gate-Reference).
 
-Additionally, Tsim supports the instructions `T`, `T_DAG`, `R_Z`, `R_X`, `R_Y`, and `U3`.
+Additionally, Tsim supports the instructions `T`, `T_DAG`, `R_Z`, `R_X`, `R_Y`, `U3`, `TPP`, `TPP_DAG`, `R_XX`, `R_YY`, `R_ZZ`, and `R_PAULI`.
 
 ```python
 import tsim

From cc7933573eb8ea518b8ee23425b1d74d1058e659 Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 21:27:27 -0400
Subject: [PATCH 12/14] Reject duplicate R_PAULI targets per product (#142)

---
 src/tsim/core/parse.py | 39 +++++++++++++++++++++++++++++++--------
 1 file changed, 31 insertions(+), 8 deletions(-)

diff --git a/src/tsim/core/parse.py b/src/tsim/core/parse.py
index 6dd4251b..f3632684 100644
--- a/src/tsim/core/parse.py
+++ b/src/tsim/core/parse.py
@@ -108,6 +108,36 @@ def parse_parametric_tag(
 }
 
 
+def _validate_r_pauli_targets(instruction: stim.CircuitInstruction) -> None:
+    """Validate an ``R_PAULI``-tagged SPP instruction's raw targets.
+
+    Rejects duplicate target qubits within a single Pauli product (before any
+    algebraic simplification) and enforces the per-instruction qubit limit, both
+    checked on the raw (un-reduced) target list.
+    """
+    targets = instruction.targets_copy()
+    total_qubits = sum(1 for t in targets if not t.is_combiner)
+    if total_qubits > R_PAULI_MAX_QUBITS:
+        raise ValueError(
+            f"R_PAULI supports at most {R_PAULI_MAX_QUBITS} qubits per instruction, "
+            f"got {total_qubits}."
+        )
+
+    seen: set[int] = set()
+    for i, target in enumerate(targets):
+        if target.is_combiner:
+            continue
+        if target.value in seen:
+            raise ValueError(
+                f"R_PAULI target qubits must be distinct within a product, "
+                f"got repeated qubit {target.value} in {str(instruction)!r}."
+            )
+        seen.add(target.value)
+        next_idx = i + 1
+        if next_idx >= len(targets) or not targets[next_idx].is_combiner:
+            seen = set()
+
+
 def _iter_pauli_products(
     instruction: stim.CircuitInstruction,
 ) -> Iterator[tuple[list[tuple[Literal["X", "Y", "Z"], int]], bool]]:
@@ -238,14 +268,7 @@ def parse_stim_circuit(
             parsed = parse_parametric_tag(instruction)
             if parsed is not None and parsed[0] == "R_PAULI":
                 params = parsed[1]
-                n_qubits = len(
-                    {t.value for t in instruction.targets_copy() if not t.is_combiner}
-                )
-                if n_qubits > R_PAULI_MAX_QUBITS:
-                    raise ValueError(
-                        f"R_PAULI supports at most {R_PAULI_MAX_QUBITS} qubits per "
-                        f"instruction, got {n_qubits}."
-                    )
+                _validate_r_pauli_targets(instruction)
                 is_dag = name == "SPP_DAG"
                 for paulis, invert in _iter_pauli_products(instruction):
                     r_pauli(b, paulis, params["theta"], dagger=is_dag ^ invert)

From 1938dc8074d4d26a5b65cdd9a7998dd14bee7eba Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 21:29:03 -0400
Subject: [PATCH 13/14] Round-trip multi-factor R_PAULI products (#142)

---
 src/tsim/utils/program_text.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/tsim/utils/program_text.py b/src/tsim/utils/program_text.py
index 47cecec1..a994de55 100644
--- a/src/tsim/utils/program_text.py
+++ b/src/tsim/utils/program_text.py
@@ -142,7 +142,7 @@ def replace_pauli_pair(m: re.Match) -> str:
         return f"R_{pauli}{pauli}({alpha}) {q0} {q1}"
 
     text = re.sub(
-        rf"\bSPP\[R_PAULI\(theta=({FLOAT_RE})\*pi\)\] ([XYZ])(\d+)\*\2(\d+)",
+        rf"\bSPP\[R_PAULI\(theta=({FLOAT_RE})\*pi\)\] ([XYZ])(\d+)\*\2(\d+)(?!\*)\b",
         replace_pauli_pair,
         text,
     )

From 1e1b74fec0aefebb707d73f556516e7e5a7295bd Mon Sep 17 00:00:00 2001
From: Mostafa Atallah <m.atallah.2014@gmail.com>
Date: Wed, 3 Jun 2026 21:30:19 -0400
Subject: [PATCH 14/14] Test R_PAULI round-trip and duplicate rejection (#142)

---
 test/unit/test_circuit.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/test/unit/test_circuit.py b/test/unit/test_circuit.py
index 4d5dd829..0ed39ad9 100644
--- a/test/unit/test_circuit.py
+++ b/test/unit/test_circuit.py
@@ -785,6 +785,19 @@ def test_inverse_r_pauli():
     assert unitaries_equal_up_to_global_phase(combined, np.eye(combined.shape[0]))
 
 
+def test_r_pauli_duplicate_target_in_product_rejected():
+    """Repeated qubits within one R_PAULI product are rejected before simplification."""
+    with pytest.raises(ValueError, match="distinct"):
+        Circuit("R_PAULI(0.25) X0*X0").get_graph()
+
+
+def test_r_pauli_long_product_roundtrip():
+    """A same-axis product with >2 factors round-trips as R_PAULI, not a mangled R_XX."""
+    c = Circuit("R_PAULI(0.3) X0*X1*X2")
+    assert str(c) == "R_PAULI(0.3) X0*X1*X2"
+    assert Circuit(str(c)) == c
+
+
 def test_inverse_mixed_circuit():
     c = Circuit("H 0\nT 0\nR_Z(0.22) 0\nCNOT 0 1")
     c_inv = c.inverse()
