feat: add parametric Pauli rotation gates

src/tsim/circuit.py

@@ -183,6 +183,49 @@ 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"):
+                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."
+                    )
+                t_list: list[int] = (
+                    [targets] if isinstance(targets, int) else list(targets)  # type: ignore[arg-type]
+                )
+                if len(t_list) != 2:
+                    raise ValueError(
+                        f"For {name} gates, exactly two target qubits must be provided."
+                    )
+                q0, q1 = (
+                    t.value if isinstance(t, stim.GateTarget) else t for t in t_list
+                )
+                if q0 == q1:
+                    raise ValueError(
+                        f"Duplicate target qubits in {name}: both targets are qubit {q0}."
+                    )
+                pauli = name[2]  # 'X', 'Y', or 'Z'
+                target_fn = {
+                    "X": stim.target_x,
+                    "Y": stim.target_y,
+                    "Z": stim.target_z,
+                }[pauli]
+                tag = f"{name}(theta={args[0]}*pi)"
+                targets = [target_fn(q0), stim.target_combiner(), target_fn(q1)]
+                name = "SPP"
+                arg = None
+            elif name == "R_PAULI":
+                if arg is None:
+                    raise ValueError("For R_PAULI gates, an angle must be provided.")
+                args = list(arg) if isinstance(arg, Iterable) else [arg]
+                if len(args) != 1:
+                    raise ValueError(
+                        "For R_PAULI gates, a single angle must be provided."
+                    )
+                tag = f"R_PAULI(theta={args[0]}*pi)"
+                name = "SPP"
+                arg = None
 
             self._stim_circ.append(name=name, targets=targets, arg=arg, tag=tag)  # type: ignore
         else:
@@ -806,6 +849,25 @@ def fix_tags(circuit: stim.Circuit) -> stim.Circuit:
                         result.append("I", targets, args, tag=new_tag)
                         continue
 
+                # Stim flips SPP↔SPP_DAG during inverse.  For parametric
+                # Pauli rotations the angle already encodes direction, so
+                # undo the flip and negate θ instead.
+                if name in ("SPP", "SPP_DAG") and tag and tag != "T":
+                    parsed = parse_parametric_tag(instr)
+                    if parsed is not None:
+                        gate_name, params = parsed
+                        if name == "SPP_DAG":
+                            # Stim flipped SPP → SPP_DAG; undo flip by negating θ.
+                            theta = float(-params["theta"])
+                            new_tag = f"{gate_name}(theta={theta}*pi)"
+                            result.append(
+                                "SPP", instr.targets_copy(), args, tag=new_tag
+                            )
+                        else:
+                            # Stim flipped SPP_DAG → SPP; tag already correct.
+                            result.append(instr)
+                        continue
+
                 result.append(instr)
             return result
 

src/tsim/core/instructions.py

@@ -963,6 +963,31 @@ def _pauli_product_phase(
             s(b, qubit)
 
 
+def r_pauli(
+    b: GraphRepresentation,
+    paulis: list[tuple[Literal["X", "Y", "Z"], int]],
+    phase: Fraction,
+) -> None:
+    """Apply exp(-i * phase * pi/2 * P) for an arbitrary Pauli product P.
+
+    Generalises R_X / R_Y / R_Z to multi-qubit Pauli strings.
+    ``R_XX``, ``R_YY``, ``R_ZZ`` are two-qubit special cases.
+
+    Args:
+        b: The graph representation to modify.
+        paulis: List of (pauli_type, qubit) pairs defining the Pauli product P.
+        phase: Rotation angle in units of π  (α means exp(-i α π/2 P)).
+
+    """
+    _pauli_product_phase(
+        b,
+        paulis,
+        lambda b_, q: r_z(b_, q, phase),
+        lambda b_, q: r_z(b_, q, -phase),
+        dagger=False,
+    )
+
+
 def spp(
     b: GraphRepresentation,
     paulis: list[tuple[Literal["X", "Y", "Z"], int]],

src/tsim/core/parse.py

@@ -16,6 +16,7 @@
     mpad,
     mpp,
     observable_include,
+    r_pauli,
     r_x,
     r_y,
     r_z,
@@ -29,6 +30,10 @@
     "R_X": frozenset({"theta"}),
     "R_Y": frozenset({"theta"}),
     "R_Z": frozenset({"theta"}),
+    "R_XX": frozenset({"theta"}),
+    "R_YY": frozenset({"theta"}),
+    "R_ZZ": frozenset({"theta"}),
+    "R_PAULI": frozenset({"theta"}),
     "U3": frozenset({"theta", "phi", "lambda"}),
 }
 
@@ -38,7 +43,7 @@ def parse_parametric_tag(
 ) -> tuple[str, dict[str, Fraction]] | None:
     """Parse the parametric tag on an instruction (e.g. ``I[R_Z(theta=0.3*pi)]``).
 
-    Supports gates: R_Z, R_X, R_Y, U3.
+    Supports gates: R_Z, R_X, R_Y, R_XX, R_YY, R_ZZ, R_PAULI, U3.
 
     Args:
         instruction: The stim instruction whose tag will be parsed.
@@ -225,6 +230,17 @@ def parse_stim_circuit(
             for paulis, invert in _iter_pauli_products(instruction):
                 mpp(b, paulis, invert, p=p)
             continue
+        if name in ("SPP", "SPP_DAG") and instruction.tag and instruction.tag != "T":
+            result = parse_parametric_tag(instruction)
+            if result is not None:
+                gate_name, params = result
+                is_dag = name == "SPP_DAG"
+                for paulis, invert in _iter_pauli_products(instruction):
+                    phase = params["theta"]
+                    if is_dag ^ invert:
+                        phase = -phase
+                    r_pauli(b, paulis, phase)
+                continue
         if name in ("SPP", "SPP_DAG") and instruction.tag == "T":
             is_dag = name == "SPP_DAG"
             for paulis, invert in _iter_pauli_products(instruction):

src/tsim/utils/clifford.py

@@ -6,7 +6,7 @@
 
 import stim
 
-from tsim.core.parse import parse_parametric_tag
+from tsim.core.parse import _iter_pauli_products, parse_parametric_tag
 
 # Clifford decompositions for U3(θ, φ, λ) = R_Z(φ) · R_Y(θ) · R_Z(λ).
 # Keys: (θ_idx, φ_idx, λ_idx) where each index ∈ {0,1,2,3} is the angle in half-pi units.
@@ -84,6 +84,20 @@ def parametric_to_clifford_gates(
         table = {"R_Z": RZ_CLIFFORD, "R_X": RX_CLIFFORD, "R_Y": RY_CLIFFORD}[gate_name]
         return [table[idx]]
 
+    if gate_name in ("R_XX", "R_YY", "R_ZZ", "R_PAULI"):
+        idx = _to_half_pi_index(params["theta"])
+        if idx is None:
+            return None
+        # Clifford-angle Pauli rotations map to SPP / SPP_DAG / identity
+        # idx=0 → I, idx=1 → SPP, idx=2 → SPP·SPP, idx=3 → SPP_DAG
+        _SPP_CLIFFORD: dict[int, list[str]] = {
+            0: [],
+            1: ["SPP"],
+            2: ["SPP", "SPP"],
+            3: ["SPP_DAG"],
+        }
+        return _SPP_CLIFFORD[idx]
+
     if gate_name == "U3":
         theta_idx = _to_half_pi_index(params["theta"])
         phi_idx = _to_half_pi_index(params["phi"])
@@ -119,6 +133,13 @@ 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 and instr.tag != "T":
+            result = parse_parametric_tag(instr)
+            if result is not None:
+                _, params = result
+                if not is_half_pi_multiple(params["theta"]):
+                    return False
+
         if instr.name == "I" and instr.tag:
             result = parse_parametric_tag(instr)
             if result is None:
@@ -154,6 +175,11 @@ def expand_clifford_rotations(source: stim.Circuit) -> stim.Circuit:
                 )
             )
             continue
+        spp_exp = _try_spp_clifford_expansion(instr)
+        if spp_exp is not None:
+            for gate_name, gate_targets in spp_exp:
+                out.append(gate_name, gate_targets, [])
+            continue
         expansion = _try_clifford_expansion(instr)
         if expansion is not None:
             gates, targets = expansion
@@ -164,6 +190,38 @@ def expand_clifford_rotations(source: stim.Circuit) -> stim.Circuit:
     return out
 
 
+def _try_spp_clifford_expansion(
+    instr: stim.CircuitInstruction,
+) -> list[tuple[str, list[stim.GateTarget]]] | None:
+    """Try to expand a tagged ``SPP`` instruction into Clifford SPP/SPP_DAG.
+
+    Returns:
+        List of ``(gate_name, targets)`` pairs, or ``None`` if the instruction
+        is not an expandable parametric Pauli rotation.
+
+    """
+    if instr.name not in ("SPP", "SPP_DAG") or not instr.tag or instr.tag == "T":
+        return None
+
+    parsed = parse_parametric_tag(instr)
+    if parsed is None:
+        return None
+
+    _, params = parsed
+    effective_params = dict(params)
+    is_dag = instr.name == "SPP_DAG"
+    for _, invert in _iter_pauli_products(instr):
+        if is_dag ^ invert:
+            effective_params["theta"] = -effective_params["theta"]
+        break
+    gates = parametric_to_clifford_gates(parsed[0], effective_params)
+    if gates is None:
+        return None
+
+    targets = instr.targets_copy()
+    return [(g, targets) for g in gates] if gates else []
+
+
 def _try_clifford_expansion(
     instr: stim.CircuitInstruction,
 ) -> tuple[list[str], list[int]] | None:

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_(?:XX|YY|ZZ|PAULI|[XYZ])\([^)]*\)|R_(?:XX|YY|ZZ|PAULI|[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.25) 0 1      → SPP[R_XX(theta=0.25*pi)] X0*X1
+        R_YY(0.25) 0 1      → SPP[R_YY(theta=0.25*pi)] Y0*Y1
+        R_ZZ(0.25) 0 1      → SPP[R_ZZ(theta=0.25*pi)] Z0*Z1
+        R_PAULI(0.25) X0*Y1 → SPP[R_PAULI(theta=0.25*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,36 @@ 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)
 
+    # R_XX/R_YY/R_ZZ must come before R_X/R_Y/R_Z to avoid partial matches
+    def replace_r_pp(m: re.Match) -> str:
+        pauli = m.group(1)  # "XX", "YY", or "ZZ"
+        angle = float(m.group(2))
+        q0 = m.group(3)
+        q1 = m.group(4)
+        if q0 == q1:
+            raise ValueError(
+                f"Duplicate target qubits in R_{pauli}: both targets are qubit {q0}."
+            )
+        p = pauli[0]
+        return f"SPP[R_{pauli}(theta={angle}*pi)] {p}{q0}*{p}{q1}"
+
+    text = re.sub(
+        rf"\bR_(XX|YY|ZZ)\(({FLOAT_RE})\)\s+(\d+)\s+(\d+)",
+        replace_r_pp,
+        text,
+    )
+
+    def replace_r_pauli(m: re.Match) -> str:
+        angle = float(m.group(1))
+        targets = m.group(2)
+        return f"SPP[R_PAULI(theta={angle}*pi)] {targets}"
+
+    text = re.sub(
+        rf"\bR_PAULI\(({FLOAT_RE})\)\s+((?:[XYZ]\d+)(?:\*[XYZ]\d+)*)",
+        replace_r_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 +122,10 @@ 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_XX(theta=α*pi)] X0*X1 → R_XX(α) 0 1
+    - SPP[R_YY(theta=α*pi)] Y0*Y1 → R_YY(α) 0 1
+    - SPP[R_ZZ(theta=α*pi)] Z0*Z1 → R_ZZ(α) 0 1
+    - SPP[R_PAULI(theta=α*pi)] X0*Y1*Z2 → R_PAULI(α) X0*Y1*Z2
     - SPP[T] → TPP
     - SPP_DAG[T] → TPP_DAG
     - S[T] → T
@@ -115,6 +155,32 @@ def replace_rotation(m: re.Match) -> str:
         text,
     )
 
+    # Replace SPP[R_XX/R_YY/R_ZZ(...)] with R_XX/R_YY/R_ZZ(...)
+    def replace_spp_r_pp(m: re.Match) -> str:
+        pauli = m.group(1)
+        angle = m.group(2)
+        q0 = m.group(3)
+        q1 = m.group(4)
+        return f"R_{pauli}({angle}) {q0} {q1}"
+
+    text = re.sub(
+        rf"\bSPP\[R_(XX|YY|ZZ)\(theta=({FLOAT_RE})\*pi\)\]\s+[XYZ](\d+)\*[XYZ](\d+)",
+        replace_spp_r_pp,
+        text,
+    )
+
+    # Replace SPP[R_PAULI(...)] with R_PAULI(...)
+    def replace_spp_r_pauli(m: re.Match) -> str:
+        angle = m.group(1)
+        targets = m.group(2)
+        return f"R_PAULI({angle}) {targets}"
+
+    text = re.sub(
+        rf"\bSPP\[R_PAULI\(theta=({FLOAT_RE})\*pi\)\]\s+((?:[XYZ]\d+)(?:\*[XYZ]\d+)*)",
+        replace_spp_r_pauli,
+        text,
+    )
+
     # Replace SPP[T] and SPP_DAG[T] with TPP and TPP_DAG
     # Must come before S[T]/S_DAG[T] to avoid partial matches
     text = re.sub(r"(?<!\w)SPP_DAG\[T\](?!\w)", "TPP_DAG", text)

test/helpers/gate_matrices.py

@@ -114,6 +114,31 @@
     "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)
+    "R_XX": lambda frac: np.array(
+        [
+            [np.cos(frac * np.pi / 2), 0, 0, -1j * np.sin(frac * np.pi / 2)],
+            [0, np.cos(frac * np.pi / 2), -1j * np.sin(frac * np.pi / 2), 0],
+            [0, -1j * np.sin(frac * np.pi / 2), np.cos(frac * np.pi / 2), 0],
+            [-1j * np.sin(frac * np.pi / 2), 0, 0, np.cos(frac * np.pi / 2)],
+        ]
+    ),
+    "R_YY": lambda frac: np.array(
+        [
+            [np.cos(frac * np.pi / 2), 0, 0, 1j * np.sin(frac * np.pi / 2)],
+            [0, np.cos(frac * np.pi / 2), -1j * np.sin(frac * np.pi / 2), 0],
+            [0, -1j * np.sin(frac * np.pi / 2), np.cos(frac * np.pi / 2), 0],
+            [1j * np.sin(frac * np.pi / 2), 0, 0, np.cos(frac * np.pi / 2)],
+        ]
+    ),
+    "R_ZZ": lambda frac: np.array(
+        [
+            [np.exp(-1j * frac * np.pi / 2), 0, 0, 0],
+            [0, np.exp(1j * frac * np.pi / 2), 0, 0],
+            [0, 0, np.exp(1j * frac * np.pi / 2), 0],
+            [0, 0, 0, np.exp(-1j * frac * np.pi / 2)],
+        ]
+    ),
     "U3": lambda frac_theta, frac_phi, frac_lambda: np.array(
         [
             [