feat(gates): add RotXY in-plane single-qubit rotation (R gate)

crates/ppvm-pauli-sum/src/sum/rot1.rs

@@ -128,6 +128,20 @@ where
     }
 }
 
+impl<T: Config> RotXY<T> for PauliSum<T>
+where
+    PauliSum<T>: RotationOne<T>,
+{
+    fn r(&mut self, addr0: usize, axis_angle: <T as Config>::Coeff, theta: <T as Config>::Coeff) {
+        // R(axis_angle, θ) = RZ(axis_angle)·RX(θ)·RZ(−axis_angle). PauliSum runs
+        // in the Heisenberg picture (observables propagate backward), so the
+        // sub-rotations are emitted in reverse of the tableau's forward order.
+        self.rz(addr0, axis_angle.clone());
+        self.rx(addr0, theta);
+        self.rz(addr0, -axis_angle);
+    }
+}
+
 /// 2-bit Pauli code: 00 I, 01 X, 10 Z, 11 Y
 /// Returns \(ε, k\) so that  –i \[P_i, P_j\]/2 = ε · P_k.
 /// For every commuting pair it yields (0, 0).
@@ -236,6 +250,32 @@ mod tests {
         assert_eq!(answer, expect);
     }
 
+    #[test]
+    fn test_r() {
+        use std::f64::consts::FRAC_PI_2;
+        let theta = 2.1;
+
+        // r(axis_angle=0, θ) == rx(θ).
+        let mut via_r: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
+        via_r += ("Z", 1.0);
+        via_r.r(0, 0.0, theta);
+        let mut via_rx: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
+        via_rx += ("Z", 1.0);
+        via_rx.rx(0, theta);
+        assert!((via_r.overlap(&via_rx) - 1.0).abs() < 1e-9);
+
+        // r(axis_angle=π/2, θ) must equal ry(θ) — NOT ry(−θ). This is the case
+        // that distinguishes the Heisenberg (backward) order from the
+        // Schrödinger one: a forward-ordered impl would yield ry(−θ) here.
+        let mut via_r: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
+        via_r += ("Z", 1.0);
+        via_r.r(0, FRAC_PI_2, theta);
+        let mut via_ry: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
+        via_ry += ("Z", 1.0);
+        via_ry.ry(0, theta);
+        assert!((via_r.overlap(&via_ry) - 1.0).abs() < 1e-9);
+    }
+
     #[test]
     fn test_rz() {
         let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();

crates/ppvm-python-native/src/interface.rs

@@ -278,6 +278,12 @@ macro_rules! create_interface {
                 if truncate { self.inner.truncate(); }
             }
 
+            #[pyo3(signature = (addr0, axis_angle, theta, truncate = true))]
+            pub fn r(&mut self, addr0: usize, axis_angle: f64, theta: f64, truncate: bool) {
+                self.inner.r(addr0, axis_angle, theta);
+                if truncate { self.inner.truncate(); }
+            }
+
             // rot2
             #[pyo3(signature = (targets, theta, truncate = true))]
             pub fn rxx(&mut self, targets: Vec<usize>, theta: f64, truncate: bool) -> PyResult<()> {

crates/ppvm-python-native/src/interface_tableau.rs

@@ -186,6 +186,10 @@ macro_rules! create_interface {
                 self.inner.u3(addr0, theta, phi, lam);
             }
 
+            pub fn r(&mut self, addr0: usize, axis_angle: f64, theta: f64) {
+                self.inner.r(addr0, axis_angle, theta);
+            }
+
             // rot2
             pub fn rxx(&mut self, targets: Vec<usize>, theta: f64) -> PyResult<()> {
                 let pairs = crate::flat_pairs(&targets)?;

crates/ppvm-tableau/src/gates/rot1.rs

@@ -42,6 +42,20 @@ where
     }
 }
 
+impl<T: Config, I, C: SparseVector<Complex<T::Coeff>, I>> RotXY<T> for GeneralizedTableau<T, I, C>
+where
+    GeneralizedTableau<T, I, C>: RotationOne<T>,
+{
+    fn r(&mut self, addr0: usize, axis_angle: T::Coeff, theta: T::Coeff) {
+        // R(axis_angle, θ) = RZ(axis_angle)·RX(θ)·RZ(−axis_angle). The tableau
+        // runs in the Schrödinger picture, so the sub-rotations are applied in
+        // forward order: RZ(−axis_angle) first, then RX(θ), then RZ(axis_angle).
+        self.rz(addr0, -axis_angle.clone());
+        self.rx(addr0, theta);
+        self.rz(addr0, axis_angle);
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -112,6 +126,62 @@ mod tests {
         assert!(!tab.measure(0).unwrap());
     }
 
+    /// R(axis_angle=0, θ=0) = identity: |0⟩ stays |0⟩, no branching.
+    #[test]
+    fn test_r_identity() {
+        let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12);
+        tab.r(0, 0.0, 0.0);
+        assert_eq!(tab.coefficients.len(), 1);
+        assert!(!tab.measure(0).unwrap());
+    }
+
+    /// R(axis_angle=0, θ=π) = RX(π): flips |0⟩ → |1⟩ with no branching.
+    #[test]
+    fn test_r_axis_zero_is_rx() {
+        let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12);
+        tab.r(0, 0.0, PI);
+        assert_eq!(tab.coefficients.len(), 1);
+        assert!(tab.measure(0).unwrap());
+    }
+
+    /// R(axis_angle=π/2, θ=π) = RY(π): flips |0⟩ → |1⟩ with no branching.
+    #[test]
+    fn test_r_axis_half_pi_is_ry() {
+        let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12);
+        tab.r(0, FRAC_PI_2, PI);
+        assert_eq!(tab.coefficients.len(), 1);
+        assert!(tab.measure(0).unwrap());
+    }
+
+    /// A partial rotation about an in-plane axis branches into two terms.
+    #[test]
+    fn test_r_branches() {
+        let mut tab: TestTableau = GeneralizedTableau::new(1, 1e-12);
+        tab.r(0, 0.37 * PI, FRAC_PI_2);
+        assert_eq!(tab.coefficients.len(), 2);
+    }
+
+    /// R(axis_angle, θ) must give identical per-seed measurement statistics
+    /// to the manual decomposition RZ(axis_angle)·RX(θ)·RZ(−axis_angle).
+    #[test]
+    fn test_r_matches_rz_rx_rz() {
+        let (axis_angle, theta) = (0.21 * PI, 0.34 * PI);
+
+        let mut tab_r: TestTableau = GeneralizedTableau::new_with_seed(1, 1e-12, 0);
+        tab_r.r(0, axis_angle, theta);
+
+        let mut tab_manual: TestTableau = GeneralizedTableau::new_with_seed(1, 1e-12, 0);
+        tab_manual.rz(0, -axis_angle);
+        tab_manual.rx(0, theta);
+        tab_manual.rz(0, axis_angle);
+
+        for seed in 0..200 {
+            let result_r = tab_r.fork(Some(seed)).measure(0).unwrap();
+            let result_manual = tab_manual.fork(Some(seed)).measure(0).unwrap();
+            assert_eq!(result_r, result_manual, "mismatch at seed {}", seed);
+        }
+    }
+
     #[test]
     fn test_two_qubit_case() {
         let mut tab: TestTableau = GeneralizedTableau::new(2, 1e-10);

crates/ppvm-traits/src/traits/branch/mod.rs

@@ -10,7 +10,7 @@ mod u3;
 
 pub use crx::CRx;
 pub use proj::Projection;
-pub use rot1::RotationOne;
+pub use rot1::{RotXY, RotationOne};
 pub use rot2::RotationTwo;
 pub use tgate::TGate;
 pub use u3::U3Gate;

crates/ppvm-traits/src/traits/branch/rot1.rs

@@ -43,3 +43,13 @@ pub trait RotationOne<T: Config> {
         }
     }
 }
+
+/// Rotation about an axis in the x/y plane:
+/// `R(axis_angle, θ) = exp(-i θ/2 · (cos(axis_angle)·X + sin(axis_angle)·Y))`.
+///
+/// The in-plane axis is `X` rotated about `Z` by `axis_angle`, so
+/// `R(axis_angle, θ) = RZ(axis_angle)·RX(θ)·RZ(−axis_angle)`.
+pub trait RotXY<T: Config> {
+    /// `R(axis_angle, θ)` on qubit `addr0`.
+    fn r(&mut self, addr0: usize, axis_angle: T::Coeff, theta: T::Coeff);
+}

crates/ppvm-traits/src/traits/mod.rs

@@ -15,7 +15,7 @@ mod strategy;
 mod trace;
 mod word_trait;
 
-pub use branch::{CRx, Projection, RotationOne, RotationTwo, TGate, U3Gate};
+pub use branch::{CRx, Projection, RotXY, RotationOne, RotationTwo, TGate, U3Gate};
 pub use clifford::{Clifford, CliffordBatch, CliffordExtensions, CliffordExtensionsBatch};
 pub use coefficient::{Coefficient, ComplexCoefficient};
 pub use hash::HashFinalize;

ppvm-python/src/ppvm/mixins.py

@@ -221,6 +221,22 @@ def rz(self, *args: Any, theta: float | None = None) -> None:
         targets, theta = _split_targets_parameter(args, theta, "theta")
         self._interface.rz(targets, theta)
 
+    def r(self, addr0: int, axis_angle: float, theta: float) -> None:
+        """Apply a rotation about an axis in the X-Y plane to the specified qubit.
+
+        ```math
+        R(\\phi, \\theta) = e^{-i \\frac{\\theta}{2} (\\cos\\phi\\, X + \\sin\\phi\\, Y)}
+            = R_Z(\\phi) R_X(\\theta) R_Z(-\\phi)
+        ```
+
+        Args:
+            addr0: The index of the target qubit.
+            axis_angle: The angle ``φ`` (in radians) of the rotation axis
+                within the X-Y plane, measured from the X-axis.
+            theta: The rotation angle in radians.
+        """
+        self._interface.r(addr0, axis_angle, theta)
+
     # Two qubit rotations
     def rxx(self, *args: Any, theta: float | None = None) -> None:
         """Apply RXX (Ising XX) rotation gates over consecutive qubit pairs.
@@ -601,6 +617,21 @@ def rz(self, *args: Any, theta: float | None = None, truncate: bool = True) -> N
         targets, theta, truncate = _split_targets_parameter_truncate(args, theta, "theta", truncate)
         self._interface.rz(targets, theta, truncate)
 
+    def r(self, addr0: int, axis_angle: float, theta: float, *, truncate: bool = True) -> None:
+        """Apply a rotation about an axis in the X-Y plane to the specified qubit.
+
+        See `RotationsMixin.r` for the gate definition.
+
+        Args:
+            addr0: The index of the target qubit.
+            axis_angle: The angle ``φ`` (in radians) of the rotation axis
+                within the X-Y plane, measured from the X-axis.
+            theta: The rotation angle in radians.
+            truncate: See `rx`.
+        """
+        self._interface.r(addr0, axis_angle, theta, truncate=truncate)
+
+    # Two qubit rotations
     def rxx(self, *args: Any, theta: float | None = None, truncate: bool = True) -> None:
         """Apply RXX (Ising XX) rotation gates over consecutive qubit pairs.
 

ppvm-python/test/generalized_tableau/test_basics.py

@@ -155,6 +155,20 @@ def test_ry_pi_equals_y():
     assert tab.measure(0)
 
 
+def test_r_axis_zero_equals_rx():
+    # r(axis_angle=0, θ=π) = rx(π), so |0> → |1>
+    tab = GeneralizedTableau(2)
+    tab.r(0, 0.0, math.pi)
+    assert tab.measure(0)
+
+
+def test_r_axis_half_pi_equals_ry():
+    # r(axis_angle=π/2, θ=π) = ry(π), so |0> → |1>
+    tab = GeneralizedTableau(2)
+    tab.r(0, math.pi / 2, math.pi)
+    assert tab.measure(0)
+
+
 def test_clifford_extensions_sqrt_x():
     # sqrt_x followed by sqrt_x_dag is identity
     tab = GeneralizedTableau(2)

ppvm-python/test/test_basics.py

@@ -191,6 +191,16 @@ def t(state):
     s.rz(0, theta=-PI / 2)
     assert pytest.approx(t(s).get("YI", 0.0)) == 1.0
 
+    # r(axis_angle=0, θ=π/2) = rx(π/2): ZI → YI
+    s = PauliSum(n_qubits=2, initial_terms=["ZI"], coefficients=[1.0])
+    s.r(0, 0.0, PI / 2)
+    assert pytest.approx(t(s).get("YI", 0.0)) == 1.0
+
+    # r(axis_angle=π/2, θ=π/2) = ry(π/2): ZI → −XI
+    s = PauliSum(n_qubits=2, initial_terms=["ZI"], coefficients=[1.0])
+    s.r(0, PI / 2, PI / 2)
+    assert pytest.approx(t(s).get("XI", 0.0)) == -1.0
+
     # rxx(π/2): ZI → YX  [cos·ZI + sin·YX at θ=π/2 → YX]
     s = PauliSum(n_qubits=2, initial_terms=["ZI"], coefficients=[1.0])
     s.rxx(0, 1, theta=PI / 2)