ppvm_pauli_sum/sum/
rot1.rs

1// SPDX-FileCopyrightText: 2026 The PPVM Authors
2// SPDX-License-Identifier: Apache-2.0
3
4use crate::sum::PauliSum;
5use ppvm_traits::char::Pauli;
6use ppvm_traits::config::Config;
7use ppvm_traits::traits::*;
8
9pub(crate) fn rotate_1_map_insert_closure<T: Config>(
10    k: &T::PauliWordType,
11    v: &mut T::Coeff,
12    axis: Pauli,
13    addr0: usize,
14    sin: &T::Coeff,
15    cos: &T::Coeff,
16) -> Option<(T::PauliWordType, T::Coeff)> {
17    if axis == Pauli::L {
18        panic!("Rotation axis cannot be L");
19    }
20    // Check loss first — avoids reading xbit/zbit on lost qubits
21    if k.get_lbit(addr0) {
22        return None;
23    }
24    let xbit = k.get_xbit(addr0);
25    let zbit = k.get_zbit(addr0);
26    // Reconstruct 2-bit Pauli code (I=00, X=01, Z=10, Y=11) from bits
27    let p_g = (zbit as u8) << 1 | (xbit as u8);
28    let (eps, p_q) = levi_civita(p_g, axis as u8);
29    if eps == 0 {
30        None
31    } else {
32        let mut coeff = v.clone();
33        *v *= cos.clone();
34        let mut new_word = k.clone();
35        new_word.set_xbit(addr0, p_q & 0b01 != 0);
36        new_word.set_zbit(addr0, p_q & 0b10 != 0);
37        new_word.rehash();
38
39        coeff *= sin.mul_sign(eps);
40        Some((new_word, coeff))
41    }
42}
43
44impl<T: Config> RotationOne<T> for PauliSum<T>
45where
46    T::Coeff: std::ops::MulAssign,
47    T::Map: ACMapInsert<T::Storage, T::Coeff, T::BuildHasher, T::PauliWordType> + ACMapConsume,
48{
49    fn rotate_1(&mut self, axis: Pauli, addr0: usize, theta: <T as Config>::Coeff) {
50        let (sin, cos) = theta.sin_cos();
51        self.map_insert(|k, v| rotate_1_map_insert_closure::<T>(k, v, axis, addr0, &sin, &cos));
52    }
53
54    #[inline]
55    fn rx(&mut self, addr0: usize, theta: impl Into<T::Coeff>) {
56        // Axis = X (xbit=1, zbit=0). Commutes when zbit==false (I or X).
57        // Anticommuting: Z(xbit=0)→new Y, eps=+1; Y(xbit=1)→new Z, eps=-1
58        let (sin, cos) = theta.into().sin_cos();
59        self.map_insert(|k, v| {
60            if k.get_lbit(addr0) {
61                return None;
62            }
63            let zbit = k.get_zbit(addr0);
64            if !zbit {
65                return None;
66            }
67            let xbit = k.get_xbit(addr0);
68            let mut coeff = v.clone();
69            *v *= cos.clone();
70            let mut new_word = k.clone();
71            new_word.set_xbit(addr0, !xbit);
72            new_word.rehash();
73            let eps: i8 = if xbit { -1 } else { 1 };
74            coeff *= sin.mul_sign(eps);
75            Some((new_word, coeff))
76        });
77    }
78
79    #[inline]
80    fn ry(&mut self, addr0: usize, theta: impl Into<T::Coeff>) {
81        // Axis = Y (xbit=1, zbit=1). Commutes when xbit==zbit (I or Y).
82        // Anticommuting: X(zbit=0)→new Z, eps=+1; Z(zbit=1)→new X, eps=-1
83        let (sin, cos) = theta.into().sin_cos();
84        self.map_insert(|k, v| {
85            if k.get_lbit(addr0) {
86                return None;
87            }
88            let xbit = k.get_xbit(addr0);
89            let zbit = k.get_zbit(addr0);
90            if xbit == zbit {
91                return None;
92            }
93            let mut coeff = v.clone();
94            *v *= cos.clone();
95            let mut new_word = k.clone();
96            new_word.set_xbit(addr0, !xbit);
97            new_word.set_zbit(addr0, !zbit);
98            new_word.rehash();
99            let eps: i8 = if zbit { -1 } else { 1 };
100            coeff *= sin.mul_sign(eps);
101            Some((new_word, coeff))
102        });
103    }
104
105    #[inline]
106    fn rz(&mut self, addr0: usize, theta: impl Into<T::Coeff>) {
107        // Axis = Z (xbit=0, zbit=1). Commutes when xbit==false (I or Z).
108        // Anticommuting: X(zbit=0)→new Y, eps=-1; Y(zbit=1)→new X, eps=+1
109        let (sin, cos) = theta.into().sin_cos();
110        self.map_insert(|k, v| {
111            if k.get_lbit(addr0) {
112                return None;
113            }
114            let xbit = k.get_xbit(addr0);
115            if !xbit {
116                return None;
117            }
118            let zbit = k.get_zbit(addr0);
119            let mut coeff = v.clone();
120            *v *= cos.clone();
121            let mut new_word = k.clone();
122            new_word.set_zbit(addr0, !zbit);
123            new_word.rehash();
124            let eps: i8 = if zbit { 1 } else { -1 };
125            coeff *= sin.mul_sign(eps);
126            Some((new_word, coeff))
127        });
128    }
129}
130
131impl<T: Config> RotXY<T> for PauliSum<T>
132where
133    PauliSum<T>: RotationOne<T>,
134{
135    fn r(&mut self, addr0: usize, axis_angle: <T as Config>::Coeff, theta: <T as Config>::Coeff) {
136        // R(axis_angle, θ) = RZ(axis_angle)·RX(θ)·RZ(−axis_angle). PauliSum runs
137        // in the Heisenberg picture (observables propagate backward), so the
138        // sub-rotations are emitted in reverse of the tableau's forward order.
139        self.rz(addr0, axis_angle.clone());
140        self.rx(addr0, theta);
141        self.rz(addr0, -axis_angle);
142    }
143}
144
145/// 2-bit Pauli code: 00 I, 01 X, 10 Z, 11 Y
146/// Returns \(ε, k\) so that  –i \[P_i, P_j\]/2 = ε · P_k.
147/// For every commuting pair it yields (0, 0).
148#[inline]
149pub fn levi_civita(i: u8, j: u8) -> (i8, u8) {
150    // --------------------------------------------------- third Pauli by XOR
151    let k = i ^ j; // 0 when i == j
152
153    // ----------- commute ⇔ i==0  OR  j==0  OR  k==0  (no false positives)
154    let commute = ((i == 0) | (j == 0) | (k == 0)) as u8; // 1 = commute
155
156    // ------------------------------------------------------ sign ε_{ijk}
157    #[inline]
158    fn rank(p: u8) -> u8 {
159        let b1 = p >> 1; // MSB
160        (b1 << 1).wrapping_sub(b1 & (p & 1)) // 0,1,2 for X,Y,Z
161    }
162
163    let ri = rank(i);
164    let rj = rank(j);
165
166    // diff = (rj - ri) mod 3   without an actual modulus
167    let mut diff = rj.wrapping_sub(ri).wrapping_add(3); // 0…5
168    diff -= 3 & (0u8.wrapping_sub(diff >> 2)); // if ≥3 subtract 3
169
170    // +1 when diff == 1,  –1 when diff == 2
171    let eps_raw = 1i8 - 2 * ((diff >> 1) as i8);
172
173    // --------------------------------------------------- zero when commute
174    let eps = eps_raw * (1 - commute as i8); // 0 if commute
175    let k = k * (1 - commute); // 0 if commute
176
177    (eps, k)
178}
179
180#[cfg(test)]
181mod tests {
182    use crate::config::fxhash::ByteF64;
183
184    use super::*;
185
186    #[test]
187    fn test_rx() {
188        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
189        answer += ("X", 1.0);
190        answer.rx(0, 2.1);
191        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
192        expect += ("X", 1.0);
193        assert_eq!(answer, expect);
194
195        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
196        answer += ("Y", 1.0);
197        answer.rx(0, 2.1);
198        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
199        expect += ("Y", 2.1_f64.cos());
200        expect += ("Z", -2.1_f64.sin());
201        assert_eq!(answer, expect);
202
203        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
204        answer += ("Z", 1.0);
205        answer.rx(0, 2.1);
206        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
207        expect += ("Z", 2.1_f64.cos());
208        expect += ("Y", 2.1_f64.sin());
209        assert_eq!(answer, expect);
210
211        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
212        answer += ("I", 1.0);
213        answer.rx(0, 2.1);
214        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
215        expect += ("I", 1.0);
216        assert_eq!(answer, expect);
217    }
218
219    #[test]
220    fn test_ry() {
221        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
222        answer += ("X", 1.0);
223        answer.ry(0, 2.1);
224
225        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
226        expect += ("X", 2.1_f64.cos());
227        expect += ("Z", 2.1_f64.sin());
228        assert_eq!(answer, expect);
229
230        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
231        answer += ("Y", 1.0);
232        answer.ry(0, 2.1);
233        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
234        expect += ("Y", 1.0);
235        assert_eq!(answer, expect);
236
237        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
238        answer += ("Z", 1.0);
239        answer.ry(0, 2.1);
240        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
241        expect += ("Z", 2.1_f64.cos());
242        expect += ("X", -2.1_f64.sin());
243        assert_eq!(answer, expect);
244
245        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
246        answer += ("I", 1.0);
247        answer.ry(0, 2.1);
248        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
249        expect += ("I", 1.0);
250        assert_eq!(answer, expect);
251    }
252
253    #[test]
254    fn test_r() {
255        use std::f64::consts::FRAC_PI_2;
256        let theta = 2.1;
257
258        // r(axis_angle=0, θ) == rx(θ).
259        let mut via_r: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
260        via_r += ("Z", 1.0);
261        via_r.r(0, 0.0, theta);
262        let mut via_rx: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
263        via_rx += ("Z", 1.0);
264        via_rx.rx(0, theta);
265        assert!((via_r.overlap(&via_rx) - 1.0).abs() < 1e-9);
266
267        // r(axis_angle=π/2, θ) must equal ry(θ) — NOT ry(−θ). This is the case
268        // that distinguishes the Heisenberg (backward) order from the
269        // Schrödinger one: a forward-ordered impl would yield ry(−θ) here.
270        let mut via_r: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
271        via_r += ("Z", 1.0);
272        via_r.r(0, FRAC_PI_2, theta);
273        let mut via_ry: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
274        via_ry += ("Z", 1.0);
275        via_ry.ry(0, theta);
276        assert!((via_r.overlap(&via_ry) - 1.0).abs() < 1e-9);
277    }
278
279    #[test]
280    fn test_rz() {
281        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
282        answer += ("X", 1.0);
283        answer.rz(0, 2.1);
284        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
285        expect += ("X", 2.1_f64.cos());
286        expect += ("Y", -2.1_f64.sin());
287        assert_eq!(answer, expect);
288
289        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
290        answer += ("Y", 1.0);
291        answer.rz(0, 2.1);
292        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
293        expect += ("Y", 2.1_f64.cos());
294        expect += ("X", 2.1_f64.sin());
295        assert_eq!(answer, expect);
296
297        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
298        answer += ("Z", 1.0);
299        answer.rz(0, 2.1);
300        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
301        expect += ("Z", 1.0);
302        assert_eq!(answer, expect);
303
304        let mut answer: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
305        answer += ("I", 1.0);
306        answer.rz(0, 2.1);
307        let mut expect: PauliSum<ByteF64<2>> = PauliSum::builder().n_qubits(1).build();
308        expect += ("I", 1.0);
309        assert_eq!(answer, expect);
310    }
311}
ppvm_pauli_sum/sum/rot1.rs

ppvm_pauli_sum/sum/
rot1.rs
