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;