paulistrings/readout/
stabilizer.rs

1//! [`StabilizerState<W>`], a reduced tableau for reading expectation values `⟨ψ|P|ψ⟩` in a stabilizer state; it is a read-out target, never evolved.
2
3use std::fmt;
4
5use num_complex::Complex64;
6
7use crate::pauli_string::PauliString;
8use crate::pauli_sum::PauliSum;
9use crate::phase::Phase;
10
11/// Why a set of generators does not define a stabilizer state.
12#[derive(Clone, Copy, Debug, PartialEq, Eq)]
13pub enum StabilizerError {
14    /// The generator count is not `num_qubits`.
15    GeneratorCount {
16        /// The required count, i.e. `num_qubits`.
17        expected: usize,
18        /// The number of generators supplied.
19        found: usize,
20    },
21    /// Generator `generator` acts on a qubit index at or beyond `num_qubits`.
22    QubitOutOfRange {
23        /// Index of the offending generator in the input slice.
24        generator: usize,
25        /// The state's qubit count.
26        num_qubits: usize,
27    },
28    /// Generators `first` and `second` anticommute, so they have no common eigenvector.
29    NotCommuting {
30        /// Index of the first of the two anticommuting generators.
31        first: usize,
32        /// Index of the second.
33        second: usize,
34    },
35    /// Generator `generator` is a product of the others up to sign, which also covers `-I ∈ S`.
36    Dependent {
37        /// Index of a generator that reduced to the identity key.
38        generator: usize,
39    },
40}
41
42impl fmt::Display for StabilizerError {
43    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
44        match self {
45            Self::GeneratorCount { expected, found } => write!(
46                f,
47                "a stabilizer state on {expected} qubits needs exactly {expected} generators, \
48                 got {found}",
49            ),
50            Self::QubitOutOfRange { generator, num_qubits } => {
51                write!(f, "generator {generator} acts on a qubit at or beyond index {num_qubits}",)
52            }
53            Self::NotCommuting { first, second } => write!(
54                f,
55                "generators {first} and {second} anticommute; stabilizer generators must \
56                 pairwise commute",
57            ),
58            Self::Dependent { generator } => write!(
59                f,
60                "generator {generator} is a product of the others (or the generators imply \
61                 -I ∈ S); stabilizer generators must be independent over GF(2)",
62            ),
63        }
64    }
65}
66
67impl std::error::Error for StabilizerError {}
68
69/// One GF(2) coordinate of a symplectic key: the `x`- or `z`-bit of one qubit.
70#[derive(Clone, Copy, Debug, PartialEq, Eq)]
71struct Column {
72    word: usize,
73    mask: u64,
74    is_z: bool,
75}
76
77impl Column {
78    fn is_set<const W: usize>(self, p: &PauliString<W>) -> bool {
79        let half = if self.is_z { &p.z } else { &p.x };
80        half[self.word] & self.mask != 0
81    }
82}
83
84/// Every symplectic coordinate of a `num_qubits`-qubit key, in elimination order.
85fn columns(num_qubits: usize) -> Vec<Column> {
86    let mut columns = Vec::with_capacity(2 * num_qubits);
87    for is_z in [false, true] {
88        for q in 0..num_qubits {
89            columns.push(Column { word: q / 64, mask: 1u64 << (q % 64), is_z });
90        }
91    }
92    columns
93}
94
95/// One row of the reduced tableau: the signed group element `(-1)^negative · key`, plus the column it pivots on.
96#[derive(Clone, Copy, Debug)]
97struct Row<const W: usize> {
98    key: PauliString<W>,
99    negative: bool,
100    pivot: Column,
101}
102
103/// A stabilizer state on `num_qubits` qubits, held as a reduced tableau of signed generators.
104///
105/// Read it out with [`crate::PauliSum::expectation_stabilizer`]; a product state reads faster through [`crate::ProductBasis`].
106#[derive(Clone, Debug)]
107pub struct StabilizerState<const W: usize> {
108    num_qubits: usize,
109    /// Echelon rows ascending in pivot column, each a signed element of the stabilizer group.
110    rows: Vec<Row<W>>,
111}
112
113impl<const W: usize> StabilizerState<W> {
114    /// Build the state stabilized by `generators`, each `(key, minus)` standing for `(-1)^minus · key`.
115    ///
116    /// Generators must number `num_qubits`, lie within `num_qubits`, pairwise commute and be independent, or a [`StabilizerError`] says which failed.
117    ///
118    /// # Panics
119    ///
120    /// Panics if `num_qubits > 64 · W`.
121    pub fn from_generators(
122        num_qubits: usize,
123        generators: &[(PauliString<W>, bool)],
124    ) -> Result<Self, StabilizerError> {
125        assert!(
126            num_qubits <= 64 * W,
127            "StabilizerState::from_generators: num_qubits {num_qubits} exceeds the {W}-word width",
128        );
129        if generators.len() != num_qubits {
130            return Err(StabilizerError::GeneratorCount {
131                expected: num_qubits,
132                found: generators.len(),
133            });
134        }
135        for (i, (key, _)) in generators.iter().enumerate() {
136            if !key.is_within(num_qubits) {
137                return Err(StabilizerError::QubitOutOfRange { generator: i, num_qubits });
138            }
139        }
140        for i in 0..generators.len() {
141            for j in (i + 1)..generators.len() {
142                if !generators[i].0.commutes_with(&generators[j].0) {
143                    return Err(StabilizerError::NotCommuting { first: i, second: j });
144                }
145            }
146        }
147
148        // `work[r] = (key, negative, origin)` is always the group element `(-1)^negative · key`; `origin` indexes the input for errors.
149        let mut work: Vec<(PauliString<W>, bool, usize)> = generators
150            .iter()
151            .enumerate()
152            .map(|(i, (key, negative))| (*key, *negative, i))
153            .collect();
154        let mut rows: Vec<Row<W>> = Vec::with_capacity(num_qubits);
155
156        for column in columns(num_qubits) {
157            let rank = rows.len();
158            let Some(p) = (rank..work.len()).find(|&r| column.is_set(&work[r].0)) else {
159                continue;
160            };
161            work.swap(rank, p);
162            let (key, negative, _) = work[rank];
163            for (r, row) in work.iter_mut().enumerate() {
164                if r == rank || !column.is_set(&row.0) {
165                    continue;
166                }
167                let phase = row.0.mul_assign(&key);
168                // Both factors are Hermitian and commute (checked above), so their product is too.
169                assert_eq!(
170                    phase.exponent() & 1,
171                    0,
172                    "a product of commuting Hermitian Paulis came out imaginary"
173                );
174                row.1 ^= negative ^ (phase == Phase::MINUS_ONE);
175            }
176            rows.push(Row { key, negative, pivot: column });
177            if rows.len() == num_qubits {
178                break;
179            }
180        }
181
182        if rows.len() < num_qubits {
183            return Err(StabilizerError::Dependent { generator: work[rows.len()].2 });
184        }
185        Ok(Self { num_qubits, rows })
186    }
187
188    /// Number of qubits the state is defined on.
189    pub fn num_qubits(&self) -> usize {
190        self.num_qubits
191    }
192
193    /// The stabilizer sign of `key`: `None` when `±key` is outside the group (expectation `0`), else `Some(true)` iff `⟨ψ|key|ψ⟩ = -1`.
194    pub(crate) fn sign_of(&self, key: &PauliString<W>) -> Option<bool> {
195        let mut reduced = *key;
196        let mut phase = Phase::ONE;
197        let mut negative = false;
198        for row in &self.rows {
199            if row.pivot.is_set(&reduced) {
200                phase += reduced.mul_assign(&row.key);
201                negative ^= row.negative;
202            }
203        }
204        if reduced != PauliString::<W>::identity() {
205            return None;
206        }
207        // Intermediate phases can be imaginary, but `key · ∏K_j = i^phase · I` with both factors Hermitian forces `i^phase = ±1`.
208        debug_assert_eq!(
209            phase.exponent() & 1,
210            0,
211            "stabilizer membership produced an imaginary phase i^{}",
212            phase.exponent(),
213        );
214        Some(negative ^ (phase == Phase::MINUS_ONE))
215    }
216}
217
218impl<const W: usize> PauliSum<W> {
219    /// Expectation value `⟨ψ|O|ψ⟩` in a stabilizer state; complex since `O` need not be Hermitian.
220    ///
221    /// # Panics
222    ///
223    /// Panics if `state.num_qubits()` differs from [`Self::num_qubits`].
224    pub fn expectation_stabilizer(&self, state: &StabilizerState<W>) -> Complex64 {
225        assert_eq!(
226            self.num_qubits(),
227            state.num_qubits(),
228            "PauliSum::expectation_stabilizer: num_qubits mismatch ({} vs {})",
229            self.num_qubits(),
230            state.num_qubits(),
231        );
232        self.sum_over_terms(|x, z, coeff| match state.sign_of(&PauliString::<W> { x: *x, z: *z }) {
233            None => Complex64::new(0.0, 0.0),
234            Some(false) => coeff,
235            Some(true) => -coeff,
236        })
237    }
238}
239
240#[cfg(test)]
241mod tests;
paulistrings/readout/stabilizer.rs

paulistrings/readout/
stabilizer.rs
