UnitaryHack 2026 - issue #142

CHANGELOG.md

@@ -8,6 +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)
 - 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.
 
 

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).

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

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:

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.
 

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,
@@ -104,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]]:
@@ -230,6 +264,15 @@ 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]
+                _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)
+                continue
         if name in ("SPP", "SPP_DAG"):
             is_dag = name == "SPP_DAG"
             for paulis, invert in _iter_pauli_products(instruction):

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:

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+)(?!\*)\b",
+        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)

test/helpers/gate_matrices.py

@@ -114,6 +114,23 @@
     "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"]),
+    # 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(
         [
             [