paulistrings/readout/
product_state.rs

1//! Product-state read-outs: [`ProductBasis`] and the [`PauliSum`] expectation value in one.
2
3use num_complex::Complex64;
4
5use crate::pauli_sum::PauliSum;
6
7/// One of the three single-qubit Pauli axes: the axis a [`ProductBasis`] qubit is an eigenstate of.
8#[derive(Clone, Copy, Debug, PartialEq, Eq)]
9pub enum PauliAxis {
10    /// `X`, with eigenstates `|+⟩` (`+1`) and `|-⟩` (`-1`).
11    X,
12    /// `Y`, with eigenstates `|+i⟩` (`+1`) and `|-i⟩` (`-1`).
13    Y,
14    /// `Z`, with eigenstates `|0⟩` (`+1`) and `|1⟩` (`-1`).
15    Z,
16}
17
18impl PauliAxis {
19    /// The axis Pauli's symplectic bits `(x, z)`.
20    const fn bits(self) -> (bool, bool) {
21        match self {
22            PauliAxis::X => (true, false),
23            PauliAxis::Y => (true, true),
24            PauliAxis::Z => (false, true),
25        }
26    }
27}
28
29/// A product of single-qubit stabilizer states, one per qubit: on qubit `q` an axis `A_q ∈ {X, Y, Z}` with a sign `s_q ∈ {+1, -1}`, i.e. `|ψ⟩ = ⊗_q |A_q, s_q⟩`.
30///
31/// Stored as bit masks in the symplectic key layout; build one with [`ProductBasis::uniform`] or [`ProductBasis::from_axes`].
32/// A term `P` has `⟨ψ|P|ψ⟩ = Π_q ⟨P_q⟩`, with `⟨P_q⟩ = 1` for `I`, `s_q` for `A_q`, and `0` for any other Pauli.
33#[derive(Clone, Debug, PartialEq, Eq)]
34pub struct ProductBasis<const W: usize> {
35    /// `x`-bit of each qubit's axis Pauli.
36    axis_x: [u64; W],
37    /// `z`-bit of each qubit's axis Pauli.
38    axis_z: [u64; W],
39    /// Sign bit per qubit: `1` selects the `-1` eigenstate.
40    negative: [u64; W],
41}
42
43impl<const W: usize> ProductBasis<W> {
44    /// The `+1` eigenstate of `axis` on every qubit: `|+…+⟩`, `|+i…+i⟩` or `|0…0⟩`.
45    pub fn uniform(axis: PauliAxis) -> Self {
46        let (axis_x, axis_z) = axis.bits();
47        let all = |b: bool| if b { [!0u64; W] } else { [0u64; W] };
48        Self { axis_x: all(axis_x), axis_z: all(axis_z), negative: [0u64; W] }
49    }
50
51    /// A basis from per-qubit `(axis, minus)` pairs: item `i` describes qubit `i`, and `minus = true` selects that axis's `-1` eigenstate.
52    ///
53    /// Qubits past the end of the iterator match only an identity factor, so supply one pair per qubit the sum uses.
54    ///
55    /// # Panics
56    ///
57    /// Panics if more than `64 * W` pairs are supplied.
58    pub fn from_axes<I>(axes: I) -> Self
59    where
60        I: IntoIterator<Item = (PauliAxis, bool)>,
61    {
62        let mut out = Self { axis_x: [0u64; W], axis_z: [0u64; W], negative: [0u64; W] };
63        for (q, (axis, minus)) in axes.into_iter().enumerate() {
64            assert!(q < 64 * W, "ProductBasis::from_axes: qubit {q} exceeds the {W}-word width",);
65            let word = q / 64;
66            let bit = 1u64 << (q % 64);
67            let (axis_x, axis_z) = axis.bits();
68            if axis_x {
69                out.axis_x[word] |= bit;
70            }
71            if axis_z {
72                out.axis_z[word] |= bit;
73            }
74            if minus {
75                out.negative[word] |= bit;
76            }
77        }
78        out
79    }
80}
81
82impl<const W: usize> PauliSum<W> {
83    /// Expectation value `⟨ψ|O|ψ⟩` in the product state `basis`; complex since `O` need not be Hermitian.
84    pub fn expectation_product(&self, basis: &ProductBasis<W>) -> Complex64 {
85        self.sum_over_terms(|x, z, coeff| {
86            // An equality on both key halves, not a subset test: an X factor on a Y-axis qubit must not match.
87            let mut mismatch = 0u64;
88            let mut sign_bits = 0u32;
89            for w in 0..W {
90                let support = x[w] | z[w];
91                mismatch |=
92                    (x[w] ^ (support & basis.axis_x[w])) | (z[w] ^ (support & basis.axis_z[w]));
93                sign_bits += (support & basis.negative[w]).count_ones();
94            }
95            match (mismatch == 0, sign_bits & 1 == 0) {
96                (false, _) => Complex64::new(0.0, 0.0),
97                (true, true) => coeff,
98                (true, false) => -coeff,
99            }
100        })
101    }
102}
103
104#[cfg(test)]
105mod tests;