From 2cf56871e744f495d24129f623fa0d3de44eb0dd Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Tue, 11 Aug 2026 20:37:05 -0400 Subject: [PATCH] fix(tableau): correlated-loss convention, duplicate batch targets, CZ-block overlap MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three defects in the crates the `ppvm` wheel actually ships, found by auditing the Lean formalization against the Rust and adjudicated with independent re-derivation. Each is user-visible today through the Python API. **1. Correlated loss disagreed between backends by a factor of two.** `p[1]` is the probability that a *named* one of the pair is lost, so P(exactly one lost) is `2·p[1]` and the survivor scales by `1 − 2·p[1] − p[0]`. That is what the paper specifies, and what `ppvm-pauli-sum` has always computed — its own test comment said `(1 - 2*p[1] - p[0])` when the channel landed. `ppvm-tableau`'s trajectory and `ppvm-tableau-sum`'s mixture instead read the trait's ambiguous "losing either one qubit" as the *total*, giving `p[1]`. So a user got one answer from `LossyPauliSum` and half of it from `GeneralizedTableau`: p LossyPauliSum GeneralizedTableau 2*p[1] [0.0, 0.3, 0.0] 0.600000 0.298450 0.60 before [0.2, 0.4, 0.0] 0.800000 0.399300 0.80 before All three backends now agree with `2·p[1]` across the admissible region, including the saturating boundary `p[0] + 2·p[1] == 1`. The wording is the actual root cause, so it is now stated normatively once, on `ppvm-traits`' `CorrelatedLossChannel`, and every other site cites it — the two tableau backends, `ppvm-pauli-sum`, `mixins.py`, `paulisum.py` (whose "losing a single qubit" was the ambiguity that let the split through), and the usage skill, which mislabelled the triple as a Pauli-error vector `[p_x, p_y, p_z]` rather than `[p_LL, p_LQ, p_LN]`. Adds `debug_assert`s for the admissible region `p[0], p[1] >= 0`, `p[0] + 2·p[1] <= 1`, `p[2] ∈ [0, 1]` (with 1e-9 slack so a saturated triple like `[1/3, 1/3, _]` is not rejected by rounding). Tests, benches and Python tests that passed out-of-domain triples are corrected to admissible ones without changing a single assertion — e.g. `[0.0, 1.0, 0.0]`, which describes a map that is not completely positive, becomes `[0.0, 0.5, 0.0]`, the same "exactly one lost every shot" witness inside the domain. **2. Duplicate qubit indices silently collapsed in batched Cliffords.** `build_masks` ORs one bit per target, so a repeated index applied the gate once instead of `k` times. `X 0 0` is legal Stim meaning apply-per-target, so `run_string("X 0 0\nM 0")` returned `Some(true)` where the truth is `Some(false)`. Detected now via a popcount against the index count, falling back to the per-index loop, which conjugates by `G^k` correctly for every family — XOR-cancelling the mask would be wrong for `s`/`sqrt_x`/`sqrt_y`, where `S² = Z ≠ I`. The fallback bodies are `#[cold] #[inline(never)]` so the fused sweeps keep their register budget: the ten bit-plane batch rows measure 0.990–1.008× against `origin/main`. **3. The fused CZ block corrupted state when pairs overlapped.** `cz_block` and `cz_block_pairs` assumed disjoint support, which the Lean proves is necessary but nothing enforced. On `X₀X₁X₂`, `cz_block(0, 1, 2)` returned `+Y₀Y₁Y₂` where the per-pair loop gives `−Y₀X₁Y₂` — wrong in both bit planes and in the sign. Since `cz_block(0, 1, n)` is adjacent-pair brickwork, the natural call was the broken one. Now falls back per pair when `offset < count`; disjoint calls stay bit-for-bit on the fused kernel. `cargo test --workspace` green (65 targets), `cargo fmt` clean, and `pytest ppvm-python/test/` 217 passed against a rebuilt wheel. Benchmarked against a pristine `origin/main` checkout across 35 rows including the untouched two-qubit and scalar gates: 0.985–1.021×, all inside the noise floor. Known follow-up, deliberately not in this PR: `ppvm-pauli-sum`'s impl is coefficient-generic, and `Coefficient` carries no ordering, so it documents the admissible region but does not assert it — meaning it still produces negative coefficients where the tableau backends now raise. The guards are also `debug_assert`s, so release wheels gain no protection; surfacing this as a `ValueError` at the binding layer is a separate change. Co-Authored-By: Claude Opus 5 --- crates/ppvm-pauli-sum/src/sum/noise.rs | 16 +- crates/ppvm-tableau-sum/src/noise.rs | 171 ++++++++- .../ppvm-tableau-sum/tests/sampler_vs_pure.rs | 32 +- crates/ppvm-tableau/benches/micro.rs | 2 +- crates/ppvm-tableau/src/data.rs | 210 ++++++++++- crates/ppvm-tableau/src/gates/clifford.rs | 326 +++++++++++++++--- crates/ppvm-tableau/src/noise.rs | 191 +++++++++- crates/ppvm-traits/src/traits/noise.rs | 23 +- ppvm-python/src/ppvm/mixins.py | 12 +- ppvm-python/src/ppvm/paulisum.py | 8 +- .../test/generalized_tableau/test_loss.py | 56 ++- skills/ppvm-usage/SKILL.md | 4 +- 12 files changed, 945 insertions(+), 106 deletions(-) diff --git a/crates/ppvm-pauli-sum/src/sum/noise.rs b/crates/ppvm-pauli-sum/src/sum/noise.rs index 9b2859a60..d6dc3aec1 100644 --- a/crates/ppvm-pauli-sum/src/sum/noise.rs +++ b/crates/ppvm-pauli-sum/src/sum/noise.rs @@ -182,13 +182,23 @@ impl LossChannel for PauliSum { impl CorrelatedLossChannel for PauliSum { /// Apply a correlated loss channel to qubits at `addr0` and `addr1`. /// - /// The three probabilities are: + /// The parameters are exactly those of + /// [`CorrelatedLossChannel`](ppvm_traits::traits::CorrelatedLossChannel), + /// which is the normative statement; this impl only restates them: /// * `p[0]`: The probability of losing both qubits simultaneously when /// both of them are in the qubit subspace. - /// * `p[1]`: The probability of losing either one qubit when both of them are - /// in the qubit subspace. + /// * `p[1]`: The probability of losing a **named** one of the two when both + /// of them are in the qubit subspace — so the probability of losing + /// *exactly one* is `2·p[1]`, split evenly between the two qubits, and the + /// survivor is scaled by `1 − 2·p[1] − p[0]`. /// * `p[2]`: The probability of losing one qubit when the other one has already /// been lost prior to the channel. + /// + /// Admissible region: `p[0], p[1] >= 0`, `p[0] + 2·p[1] <= 1`, + /// `p[2] ∈ [0, 1]`. Outside it the map is not completely positive and this + /// impl will produce negative coefficients rather than raising — unlike the + /// tableau backends, which `debug_assert`. Guarding a generic + /// `Coefficient` (which carries no ordering) is tracked separately. fn correlated_loss_channel(&mut self, addr0: usize, addr1: usize, p: [T::Coeff; 3]) { self.map_insert_multiple(|k, v| { match (k.get(addr0), k.get(addr1)) { diff --git a/crates/ppvm-tableau-sum/src/noise.rs b/crates/ppvm-tableau-sum/src/noise.rs index eedbde003..756da6a0d 100644 --- a/crates/ppvm-tableau-sum/src/noise.rs +++ b/crates/ppvm-tableau-sum/src/noise.rs @@ -10,7 +10,8 @@ use num::{ }; use ppvm_pauli_word::pattern::NotIdentity; use ppvm_tableau::{ - data::GeneralizedTableau, sparsevec::SparseVector, tableau_index::TableauIndex, + data::GeneralizedTableau, noise::is_admissible_correlated_loss, sparsevec::SparseVector, + tableau_index::TableauIndex, }; use ppvm_traits::config::Config; use ppvm_traits::traits::{ @@ -450,6 +451,10 @@ where addr1: usize, p: [::Coeff; 3], ) { + debug_assert!( + is_admissible_correlated_loss(&p), + "correlated loss needs p0, p1 >= 0, p0 + 2*p1 <= 1, p2 in [0, 1]; got {p:?}" + ); let mut branches = Vec::<(GeneralizedTableau, T::Coeff, u64, u64)>::with_capacity( 3 * self.entries.len(), ); @@ -481,7 +486,11 @@ where } // if both are present, then we create 3 new branches: - // losing both (p[0]), one, or the other qubit (p[1]) + // losing both (p[0]), one, or the other qubit. `p[1]` is the + // probability that a *named* one of the pair is lost, so each + // single-loss branch carries `p[1]` and the survivor keeps + // `1 − p[0] − 2·p[1]`. See + // `ppvm_traits::traits::CorrelatedLossChannel`. let tab_seed_both = self.rng.random::(); let mut tab_lose_both = tab.fork(Some(tab_seed_both)); @@ -503,7 +512,7 @@ where tab_lose_0.is_lost[addr0] = true; branches.push(( tab_lose_0, - p_sum.clone() * (p[1].clone() / 2.0.into()), + p_sum.clone() * p[1].clone(), word_fp, phase_loss ^ loss_mask(addr0), )); @@ -513,12 +522,12 @@ where tab_lose_1.is_lost[addr1] = true; branches.push(( tab_lose_1, - p_sum.clone() * (p[1].clone() / 2.0.into()), + p_sum.clone() * p[1].clone(), word_fp, phase_loss ^ loss_mask(addr1), )); - let p_total = p[0].clone() + p[1].clone(); + let p_total = p[0].clone() + p[1].clone() + p[1].clone(); *p_sum *= T::Coeff::one() - p_total; }); @@ -573,3 +582,155 @@ where .insert_or_merge_batch(branches, &self.sum_cutoff); } } + +#[cfg(test)] +mod tests { + // === G-040 — the correlated-loss `p[1]` convention === + // + // The paper (`ppvm-paper/main.tex:462`, `:523`, `:845`) is the definition of + // record: `p[1]` is `p_LQ`, the probability that a **named** one of the pair + // is lost, so each single-loss branch carries `p[1]`, the total weight on + // "exactly one lost" is `2·p[1]`, and the survivor keeps `1 − p[0] − 2·p[1]`. + // The same number is observable three ways — as a Heisenberg coefficient + // (`ppvm-pauli-sum`), as a branch weight (this mixture) and as a sampling + // frequency (`ppvm-tableau`'s trajectory) — and the three must agree, since + // the cross-backend disagreement is what a Python user actually hits. + + use ppvm_pauli_sum::config::fxhash::ByteF64; + use ppvm_pauli_sum::prelude::*; + use ppvm_tableau::prelude::*; + use ppvm_traits::traits::{CorrelatedLossChannel, NoStrategy}; + + use crate::data::GeneralizedTableauSum; + use crate::storage::EntryStore; + + type Cfg = ByteF64<1>; + type TabSum = GeneralizedTableauSum; + type Tab = GeneralizedTableau; + type LossyTestSum = PauliSum< + ppvm_pauli_sum::config::fxhash::Byte< + 1, + f64, + NoStrategy, + LossyPauliWord<[u8; 1], fxhash::FxBuildHasher>, + >, + >; + + /// Total mixture weight on branches with exactly one lost qubit. + /// `sum_cutoff = 0.0`, so no branch is truncated away and no + /// renormalization can hide a mis-weighted survivor. + fn mixture_single_loss_weight(p: [f64; 3]) -> f64 { + let mut sum: TabSum = GeneralizedTableauSum::new_with_seed(2, 1e-12, 0.0, 7); + sum.correlated_loss_channel(0, 1, p); + sum.entries + .iter() + .filter(|(tab, _)| tab.is_lost[0] ^ tab.is_lost[1]) + .map(|(_, probability)| *probability) + .sum() + } + + /// Total mixture weight on branches where both qubits are still present. + fn mixture_survivor_weight(p: [f64; 3]) -> f64 { + let mut sum: TabSum = GeneralizedTableauSum::new_with_seed(2, 1e-12, 0.0, 7); + sum.correlated_loss_channel(0, 1, p); + sum.entries + .iter() + .filter(|(tab, _)| !tab.is_lost[0] && !tab.is_lost[1]) + .map(|(_, probability)| *probability) + .sum() + } + + /// The Heisenberg scale factor `ppvm-pauli-sum` applies to a fully + /// in-subspace observable, i.e. `1 − p[0] − P(exactly one lost)`. + fn pauli_sum_survivor(p: [f64; 3]) -> f64 { + let mut sum = LossyTestSum::builder().n_qubits(2).build(); + sum += ("ZZ", 1.0); + sum.correlated_loss_channel(0, 1, p); + let zz: LossyPauliWord<[u8; 1], fxhash::FxBuildHasher> = "ZZ".into(); + *sum.data() + .get(&zz) + .expect("the all-present term survives a pure rescale") + } + + /// The trajectory's sampled fraction of runs that lose exactly one qubit. + fn trajectory_single_loss_fraction(p: [f64; 3], trials: u64) -> f64 { + let mut hits = 0u64; + for seed in 0..trials { + let mut tab: Tab = GeneralizedTableau::new_with_seed(2, 1e-12, seed); + tab.correlated_loss_channel(0, 1, p); + if tab.is_lost[0] ^ tab.is_lost[1] { + hits += 1; + } + } + hits as f64 / trials as f64 + } + + #[test] + fn correlated_loss_exactly_one_lost_is_two_p1_on_every_backend() { + let p1 = 0.3_f64; + let p = [0.0, p1, 0.0]; + let expected = 2.0 * p1; + + let mixture = mixture_single_loss_weight(p); + let pauli_sum = 1.0 - pauli_sum_survivor(p); + let trajectory = trajectory_single_loss_fraction(p, 20_000); + + // Report every dissenting backend, so a failure names the split rather + // than only its first symptom. + let mut wrong = Vec::new(); + for (backend, value, tolerance) in [ + ("mixture", mixture, 1e-12), + ("ppvm-pauli-sum", pauli_sum, 1e-12), + ("trajectory (20k seeds)", trajectory, 0.02), + ] { + if (value - expected).abs() >= tolerance { + wrong.push(format!(" {backend}: {value}")); + } + } + assert!( + wrong.is_empty(), + "p = {p:?}: P(exactly one lost) must be 2*p[1] = {expected} on every \ + backend, got\n{}", + wrong.join("\n") + ); + } + + #[test] + fn correlated_loss_survivor_weight_is_one_minus_p0_minus_two_p1() { + let p = [0.2_f64, 0.3, 0.4]; + let expected = 1.0 - p[0] - 2.0 * p[1]; + let mixture = mixture_survivor_weight(p); + assert!( + (mixture - expected).abs() < 1e-12, + "mixture survivor weight {mixture}, want 1 - p0 - 2*p1 = {expected}" + ); + let pauli_sum = pauli_sum_survivor(p); + assert!( + (pauli_sum - expected).abs() < 1e-12, + "ppvm-pauli-sum survivor {pauli_sum}, want {expected}" + ); + } + + // G-043 — the admissible region `p0, p1 >= 0`, `p0 + 2·p1 <= 1`, + // `p2 ∈ [0, 1]` is what makes the channel completely positive. Outside it + // the mixture truncates a negative survivor weight and renormalizes, + // silently. + + #[test] + #[cfg(debug_assertions)] + #[should_panic(expected = "p0 + 2*p1 <= 1")] + fn correlated_loss_rejects_inadmissible_probabilities() { + let mut sum: TabSum = GeneralizedTableauSum::new_with_seed(2, 1e-12, 0.0, 7); + sum.correlated_loss_channel(0, 1, [0.6, 0.6, 0.0]); + } + + /// The saturated boundary `p0 + 2·p1 == 1` is admissible and must not trip + /// the guard. + #[test] + fn correlated_loss_saturated_boundary_is_admissible() { + let mut sum: TabSum = GeneralizedTableauSum::new_with_seed(2, 1e-12, 0.0, 7); + sum.correlated_loss_channel(0, 1, [0.2, 0.4, 1.0]); + let survivor = mixture_survivor_weight([0.2, 0.4, 1.0]); + assert!(survivor.abs() < 1e-15, "survivor {survivor} should vanish"); + } +} diff --git a/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs b/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs index 2595dd339..aa7dce372 100644 --- a/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs +++ b/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs @@ -1190,10 +1190,14 @@ fn bell_pair_with_two_qubit_pauli_error_nonuniform() { // Correlated two-qubit loss channel // --------------------------------------------------------------------------- // -// Probability layout (matches the trait spec and GeneralizedTableau): +// Probability layout (normative statement: `ppvm_traits::traits::noise::CorrelatedLossChannel`): // p[0] = P(lose both | both present) -// p[1] = P(lose either one | both present) → split 50/50 across q0/q1 +// p[1] = P(lose a *named* one | both present) → P(exactly one) = 2·p[1], +// split 50/50 across q0/q1 // p[2] = P(lose remaining | the other was lost prior) +// Admissible region: p[0], p[1] >= 0, p[0] + 2·p[1] <= 1, p[2] in [0, 1]. Every +// triple below stays inside it — a triple that violates it names a channel that +// is not completely positive, and the backends now `debug_assert` against it. #[test] fn correlated_loss_channel_zero_prob_is_noop() { @@ -1229,10 +1233,10 @@ fn correlated_loss_channel_single_loss_certain_is_5050_between_qubits() { // Never both, never neither. let shots = 8000; let sum = run_sum(2, shots, 1e-12, |t| { - t.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]); + t.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]); }); let pure = run_pure(2, shots, |t| { - t.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]); + t.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]); }); // Sanity: no shot has both lost or neither lost. assert!( @@ -1268,16 +1272,16 @@ fn correlated_loss_channel_marginals_on_ground_state() { // probabilities and TVD against pure. // p[0] = 0.20, p[1] = 0.40, p[2] is unused (no qubit pre-lost). // Expected outcome probabilities: - // P(both lost) = 0.20 - // P(only q0 lost) = 0.20 (p[1]/2) - // P(only q1 lost) = 0.20 - // P(none lost) = 0.40 + // P(both lost) = 0.20 (p[0]) + // P(only q0 lost) = 0.20 (p[1]) + // P(only q1 lost) = 0.20 (p[1]) + // P(none lost) = 0.40 (1 − p[0] − 2·p[1]) let shots = 8000; let sum = run_sum(2, shots, 1e-12, |t| { - t.correlated_loss_channel(0, 1, [0.20, 0.40, 0.0]); + t.correlated_loss_channel(0, 1, [0.20, 0.20, 0.0]); }); let pure = run_pure(2, shots, |t| { - t.correlated_loss_channel(0, 1, [0.20, 0.40, 0.0]); + t.correlated_loss_channel(0, 1, [0.20, 0.20, 0.0]); }); for (label, data) in [("sum", &sum), ("pure", &pure)] { let both = data @@ -1326,11 +1330,11 @@ fn correlated_loss_channel_preexisting_loss_falls_back_to_p2() { let shots = 8000; let sum = run_sum(2, shots, 1e-12, |t| { t.loss_channel(0, 1.0); - t.correlated_loss_channel(0, 1, [0.5, 0.5, 0.3]); + t.correlated_loss_channel(0, 1, [0.5, 0.25, 0.3]); }); let pure = run_pure(2, shots, |t| { t.loss_channel(0, 1.0); - t.correlated_loss_channel(0, 1, [0.5, 0.5, 0.3]); + t.correlated_loss_channel(0, 1, [0.5, 0.25, 0.3]); }); // q0 must be lost every shot. q1 must be lost with frequency p[2]. assert!(sum.iter().all(|s| s[0].is_none())); @@ -1360,12 +1364,12 @@ fn correlated_loss_channel_both_preexisting_loss_is_noop() { let sum = run_sum(2, shots, 1e-12, |t| { t.loss_channel(0, 1.0); t.loss_channel(1, 1.0); - t.correlated_loss_channel(0, 1, [0.5, 0.5, 0.5]); + t.correlated_loss_channel(0, 1, [0.5, 0.25, 0.5]); }); let pure = run_pure(2, shots, |t| { t.loss_channel(0, 1.0); t.loss_channel(1, 1.0); - t.correlated_loss_channel(0, 1, [0.5, 0.5, 0.5]); + t.correlated_loss_channel(0, 1, [0.5, 0.25, 0.5]); }); assert!(sum.iter().all(|s| s == &vec![None, None])); assert!(pure.iter().all(|s| s == &vec![None, None])); diff --git a/crates/ppvm-tableau/benches/micro.rs b/crates/ppvm-tableau/benches/micro.rs index 0d598d7e2..85a50467f 100644 --- a/crates/ppvm-tableau/benches/micro.rs +++ b/crates/ppvm-tableau/benches/micro.rs @@ -277,7 +277,7 @@ fn bench_noise(c: &mut Criterion) { group.bench_function("correlated_loss_channel", |b| { b.iter_batched_ref( || tab.fork(None), - |t| t.correlated_loss_channel(0, 1, [0.5, 0.3, 0.2]), + |t| t.correlated_loss_channel(0, 1, [0.4, 0.3, 0.2]), criterion::BatchSize::SmallInput, ); }); diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index 7b8d008db..e4f7ab640 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -200,8 +200,23 @@ impl Tableau { /// All pairs must be in the same u64 word. This replaces N individual CZ calls /// with a single word-level shift+XOR operation per row. /// + /// # Preconditions + /// + /// `offset >= count`, i.e. the pairs have pairwise-disjoint supports (pairs + /// `i != j` collide iff `i - j == offset`, which needs `offset < count`, or + /// `offset == 0`). Both fused updates need it: the single + /// `count_ones() & 1` phase reads the pre-update `z` plane for every pair at + /// once (`Batch.lean::czSeq_phase` proves it equal to the sequential loop + /// only under pairwise disjointness, and `czSeq_phase_needs_disjoint` + /// exhibits a counterexample), and the `z` delta `OR`s the two shifted + /// planes, so a `z` bit written by two pairs is set once instead of + /// XOR-cancelled. Callers that cannot guarantee it must use + /// [`GeneralizedTableau::cz_block_pairs`] or + /// [`GeneralizedTableau::cz_block`], which fall back to the per-pair loop. + /// /// # Panics - /// Debug-asserts that all bits are within the same word. + /// Debug-asserts that all bits are within the same word and that the pairs + /// have pairwise-disjoint supports. #[inline] pub fn cz_block_pairs(&mut self, base: usize, offset: usize, count: usize) where @@ -220,6 +235,11 @@ impl Tableau { word_idx, "All CZ pairs must be in the same word" ); + debug_assert!( + offset >= count, + "cz_block pairs must have pairwise-disjoint supports \ + (Batch.lean::czSeq_phase_needs_disjoint)" + ); let one = ::Store::one(); let zero = ::Store::zero(); @@ -743,7 +763,10 @@ where } /// Apply CZ to N pairs with constant offset: (base+i, base+offset+i) for i in 0..count. - /// Falls back to individual CZ calls if any qubit in the range is lost. + /// Falls back to individual CZ calls if any qubit in the range is lost, or + /// if the pairs' supports overlap (`offset < count`, including the + /// degenerate `offset == 0`), which the fused kernel does not handle — see + /// [`Tableau::cz_block_pairs`]'s preconditions. pub fn cz_block_pairs(&mut self, base: usize, offset: usize, count: usize) where <::Store as TryFrom>::Error: Debug, @@ -752,16 +775,28 @@ where // Check if any qubit in the range is lost let any_lost = (0..count).any(|i| self.is_lost[base + i] || self.is_lost[base + offset + i]); - if !any_lost { + if !any_lost && offset >= count { self.tableau.cz_block_pairs(base, offset, count); } else { - // Fallback to individual CZ calls - for i in 0..count { - let c = base + i; - let t = base + offset + i; - if !self.is_lost[c] && !self.is_lost[t] { - Clifford::cz(&mut self.tableau, c, t); - } + self.cz_pairs_each(base, offset, count); + } + } + + /// The [`Self::cz_block_pairs`] fallback: one scalar `cz` per surviving + /// pair, in order. Outlined `#[cold]`/`#[inline(never)]` so this `count`-way + /// loop of full `2n`-row sweeps stays out of the fused kernel's register + /// budget and cache lines — the fused path keeps a single call. + #[cold] + #[inline(never)] + fn cz_pairs_each(&mut self, base: usize, offset: usize, count: usize) + where + ::Store: PrimInt, + { + for i in 0..count { + let c = base + i; + let t = base + offset + i; + if !self.is_lost[c] && !self.is_lost[t] { + Clifford::cz(&mut self.tableau, c, t); } } } @@ -807,6 +842,13 @@ where /// word) or [`Self::cz_block_pairs_cross_word`] (straddling two words), so /// callers never need to reason about the `u64` packing. CZ is symmetric, /// so the two bases may be passed in either order. + /// + /// Overlapping pairs — `|target_base - control_base| < count`, of which + /// adjacent-pair brickwork `cz_block(0, 1, n)` is the extreme case — are + /// handled: the same-word segments route through [`Self::cz_block_pairs`], + /// which falls back to the per-pair loop. A cross-word segment is always + /// disjoint (its `run` never exceeds the offset), so the fused kernel stays + /// live there. pub fn cz_block(&mut self, control_base: usize, target_base: usize, count: usize) where <::Store as TryFrom>::Error: Debug, @@ -1551,6 +1593,154 @@ mod tests { snapshot_tableau(&tab3.tableau) ); } + + /// Deterministic pseudo-random Clifford frame. `s`/`z` guarantee the rows + /// carry pre-existing `z` bits *and* non-zero phases, which is what makes + /// the fused block's `OR`ed `z` update and its single `count_ones() & 1` + /// parity observable. + fn prepare_frame(tab: &mut Tableau, n: usize, seed: u64) + where + ::Store: PrimInt, + { + let mut s = seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1; + let mut next = move || { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + (s >> 11) as usize + }; + for _ in 0..6 * n { + let q = next() % n; + match next() % 6 { + 0 => Clifford::h(tab, q), + 1 => Clifford::s(tab, q), + 2 => Clifford::z(tab, q), + 3 => Clifford::x(tab, q), + 4 => Clifford::cnot(tab, q, (q + 1 + next() % (n - 1)) % n), + _ => Clifford::cz(tab, q, (q + 1 + next() % (n - 1)) % n), + } + } + } + + /// The fused CZ block must agree with the sequential per-pair `cz` loop + /// even when the pairs' supports **overlap** (`offset < count`, or the + /// degenerate `offset == 0`) — bits *and* sign. + /// + /// `czSeq_phase` proves the fused `count_ones() & 1` phase equals the + /// sequential loop only under pairwise-disjoint supports, and + /// `czSeq_phase_needs_disjoint` exhibits the counterexample replayed below. + /// Adjacent-pair brickwork — `cz_block(0, 1, n)` — is precisely the + /// overlapping shape, so the most natural call is the broken one (G-060). + #[test] + fn test_cz_block_matches_per_pair_loop_on_overlapping_pairs() { + use ppvm_pauli_sum::config::fx64hash::Byte8F64; + type GTab = GeneralizedTableau>; + let n = 16; + // (control_base, target_base, count); offset = target - control. + let configs = [ + (0usize, 1usize, 2usize), // czSeq_phase_needs_disjoint's witness + (0, 1, 3), + (0, 1, 15), // adjacent-pair brickwork over the whole register + (3, 4, 8), + (0, 2, 5), // the ledger's second witness + (2, 5, 9), + (4, 4, 3), // degenerate: offset == 0, i.e. pairs (q, q) + (0, 1, 1), // a single pair is disjoint: must stay on the fused path + (0, 8, 8), // offset == count: the disjoint boundary + ]; + for seed in 0..8u64 { + let mut prepared: GTab = GeneralizedTableau::new(n, 1e-12); + prepare_frame(&mut prepared.tableau, n, seed); + // The frame must actually exercise the z plane and the phase. + assert!( + snapshot_tableau(&prepared.tableau) + .iter() + .any(|&(_, z, _)| z != [0]), + "frame (seed {seed}) has no pre-existing z bits" + ); + assert!( + snapshot_tableau(&prepared.tableau) + .iter() + .any(|&(.., phase)| phase != 0), + "frame (seed {seed}) has no non-zero phase" + ); + + for &(cb, tb, count) in &configs { + let mut seq = prepared.clone(); + let mut fused = prepared.clone(); + let mut fused_pairs = prepared.clone(); + + for i in 0..count { + Clifford::cz(&mut seq.tableau, cb + i, tb + i); + } + fused.cz_block(cb, tb, count); + fused_pairs.cz_block_pairs(cb, tb - cb, count); + + assert_eq!( + snapshot_tableau(&seq.tableau), + snapshot_tableau(&fused.tableau), + "cz_block({cb}, {tb}, {count}) disagrees with the per-pair loop (seed {seed})" + ); + assert_eq!( + snapshot_tableau(&seq.tableau), + snapshot_tableau(&fused_pairs.tableau), + "cz_block_pairs({cb}, {}, {count}) disagrees with the per-pair loop (seed {seed})", + tb - cb + ); + } + } + } + + /// `czSeq_phase_needs_disjoint`'s exact witness, with the absolute answer + /// derived from the Pauli algebra rather than from the loop: `CZ` + /// conjugates `X⊗X ↦ +Y⊗Y` and `Y⊗X ↦ −X⊗Y`, so `CZ₁₂·CZ₀₁` sends + /// `X₀X₁X₂` to `−Y₀X₁Y₂` — row coordinates `x = 111`, `z = 101`, phase `2`. + #[test] + fn test_cz_block_overlapping_pairs_hits_the_algebraic_answer() { + use ppvm_pauli_sum::config::fx64hash::Byte8F64; + type GTab = GeneralizedTableau>; + let mut fused: GTab = GeneralizedTableau::new(3, 1e-12); + // Row 0 is the destabilizer X₀; the two CNOTs spread it to X₀X₁X₂. + Clifford::cnot(&mut fused.tableau, 0, 1); + Clifford::cnot(&mut fused.tableau, 0, 2); + let mut seq = fused.clone(); + assert_eq!(snapshot_tableau(&fused.tableau)[0], ([0b111], [0b000], 0)); + + Clifford::cz(&mut seq.tableau, 0, 1); + Clifford::cz(&mut seq.tableau, 1, 2); + fused.cz_block(0, 1, 2); + + assert_eq!(snapshot_tableau(&seq.tableau)[0], ([0b111], [0b101], 2)); + assert_eq!(snapshot_tableau(&fused.tableau)[0], ([0b111], [0b101], 2)); + } + + /// Overlapping brickwork that straddles a storage-word boundary: the + /// segment split must not turn an overlapping run into a fused one. + #[test] + fn test_cz_block_overlapping_pairs_across_word_boundary() { + use ppvm_pauli_sum::config::fx64hash::Byte8F64; + type GTab = GeneralizedTableau>; + let n = 85; + let (control_base, target_base, count) = (60, 61, 10); + for seed in 0..8u64 { + let mut seq: GTab = GeneralizedTableau::new(n, 1e-12); + prepare_frame(&mut seq.tableau, n, seed); + let mut fused = seq.clone(); + + for i in 0..count { + Clifford::cz(&mut seq.tableau, control_base + i, target_base + i); + } + fused.cz_block(control_base, target_base, count); + + assert_eq!( + snapshot_tableau(&seq.tableau), + snapshot_tableau(&fused.tableau), + "cz_block({control_base}, {target_base}, {count}) disagrees \ + with the per-pair loop (seed {seed})" + ); + } + } + // ─── reset_all ──────────────────────────────────────────────────── /// `GeneralizedTableau::reset_all` restores the full state to a fresh diff --git a/crates/ppvm-tableau/src/gates/clifford.rs b/crates/ppvm-tableau/src/gates/clifford.rs index 60fe5345d..7420c3bd8 100644 --- a/crates/ppvm-tableau/src/gates/clifford.rs +++ b/crates/ppvm-tableau/src/gates/clifford.rs @@ -354,15 +354,102 @@ where impl_generalized_tableau_clifford_pair!(cy); } +/// How many sites a per-word mask covers — one bit per distinct index, so this +/// equals the index count exactly when the indices are distinct. +#[inline] +fn mask_bits(masks: &[S]) -> u64 { + masks.iter().map(|m| u64::from(m.count_ones())).sum() +} + +/// Take the fused per-word masks, or bail out of the caller: `return` when there +/// is nothing to do, and on a **repeated** index hand the whole batch to the +/// gate's outlined `$fallback`. +/// +/// A mask carries one bit per site, so a repeated index would make the fused +/// sweep apply the gate once where the per-index loop — the `CliffordBatch` +/// contract, and what a legal `X 0 0` in a `.stim` file means — conjugates by +/// `Gᵏ`. `Gᵏ` is the identity only for the involutory gates (`S² = Z`, +/// `(√X)² = X`, `(√Y)² = Y`), so neither one bit nor an XOR-cancelled bit is +/// right for every family; only the per-index loop is (G-061). +/// +/// The fast path costs exactly one branch on a flag `build_masks` has already +/// computed, and the slow path is a *call* — no scalar loop is inlined into the +/// fused body. Spelling the fallback out here instead (so each gate inlined a +/// whole `2n`-row sweep per index into its own preamble) cost 22–31% on the five +/// gates that map the x/z bit planes (`h`, `√X`, `√X†`, `√Y`, `√Y†`) when it was +/// measured on the `-2` crates: the bodies with real work to do lose registers +/// and I-cache to code that never runs. The phase-only gates (`x`/`y`/`z`/`s`) +/// never move, which is why A/B'ing only those rows hides it. +macro_rules! masks_or_scalar_fallback { + ($self:ident, $indices:ident, $fallback:ident) => { + match $self.build_masks($indices) { + // Distinct indices: the fused sweep is licensed. + Some((masks, n_words, false)) => (masks, n_words), + // No rows or no indices: the batch is a no-op. + None => return, + // A repeated index: owe `Gᵏ`, i.e. the per-index loop. + Some((.., true)) => return $self.$fallback($indices), + } + }; +} + +/// Define the outlined per-index loops the batched gates fall back to on a +/// repeated index (`$fallback` runs the canonical `$one` once per occurrence). +/// +/// `#[cold]` + `#[inline(never)]` is the point: each loop is a full `2n`-row +/// sweep per index, and keeping it out of line keeps it off the fused bodies' +/// register budget and out of their cache lines — the caller keeps a call. +macro_rules! repeated_index_loops { + ($($fallback:ident => $one:ident),+ $(,)?) => { + impl Tableau + where + ::Store: PrimInt, + { + $( + #[cold] + #[inline(never)] + fn $fallback(&mut self, indices: &[usize]) { + for &q in indices { + self.$one(q); + } + } + )+ + } + }; +} + +repeated_index_loops! { + x_each => x, + y_each => y, + z_each => z, + h_each => h, + s_each => s, + s_dag_each => s_dag, + sqrt_x_each => sqrt_x, + sqrt_x_dag_each => sqrt_x_dag, + sqrt_y_each => sqrt_y, + sqrt_y_dag_each => sqrt_y_dag, +} + impl Tableau where ::Store: PrimInt, { /// Build per-word bitmasks from a list of qubit indices. - /// Returns `(masks, n_words)`. Stack-allocates for up to 8 storage words; - /// spills to the heap beyond that, so there is no hard qubit cap. + /// Returns `(masks, n_words, repeated)`. Stack-allocates for up to 8 storage + /// words; spills to the heap beyond that, so there is no hard qubit cap. + /// + /// Distinctness is **checked**, not assumed, and checked here, in the third + /// tuple slot: the mask holds one bit per *distinct* index, so it is short of + /// `indices.len()` bits exactly when an index repeats. That is `n_words` + /// popcounts (typically one or two) *after* the index walk, which is left + /// byte-for-byte as it was — folding the test *into* the walk, as a + /// `repeated |= word & bit` accumulator, measured worse on the `-2` crates + /// (it adds a loop-carried dependency to the walk), and outlining this whole + /// function was worse again. `repeated` sends the caller, via + /// [`masks_or_scalar_fallback`], to its outlined per-index loop. #[inline] - fn build_masks(&self, indices: &[usize]) -> Option<(MaskBuf, usize)> { + fn build_masks(&self, indices: &[usize]) -> Option<(MaskBuf, usize, bool)> { if self.data.is_empty() || indices.is_empty() { return None; } @@ -375,7 +462,8 @@ where masks[addr0 / bits_per_word] = masks[addr0 / bits_per_word] | (one << (addr0 % bits_per_word)); } - Some((masks, n_words)) + let repeated = mask_bits(&masks) != indices.len() as u64; + Some((masks, n_words, repeated)) } } @@ -386,10 +474,7 @@ where /// `X` is bit-preserving: phase flips for each masked qubit where z=1. #[inline] fn x_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, x_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -409,10 +494,7 @@ where /// `Y` is bit-preserving: phase flips for each masked qubit where x⊕z=1. #[inline] fn y_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, y_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -433,10 +515,7 @@ where /// `Z` is bit-preserving: phase flips for each masked qubit where x=1. #[inline] fn z_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, z_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -456,10 +535,7 @@ where /// Forward `S`: phase flips where x&z=1, then z ^= x for masked qubits. #[inline] fn s_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, s_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -514,10 +590,7 @@ where /// phase += 2 when x=1 & z=1 (Y goes to -Y). #[inline] fn h_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, h_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -579,10 +652,7 @@ where /// Phase flips where x&!z=1, then z ^= x for masked qubits. #[inline] fn s_dag_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, s_dag_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -637,10 +707,7 @@ where /// reducing N individual operations to O(n_words) per row. #[inline] fn sqrt_y_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, sqrt_y_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -667,10 +734,7 @@ where /// Apply `(√Y)†` to multiple qubits using combined bitmask operations. #[inline] fn sqrt_y_dag_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, sqrt_y_dag_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -697,10 +761,7 @@ where /// Apply `√X` to multiple qubits using combined bitmask operations. #[inline] fn sqrt_x_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, sqrt_x_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -723,10 +784,7 @@ where /// Apply `(√X)†` to multiple qubits using combined bitmask operations. #[inline] fn sqrt_x_dag_many(&mut self, indices: &[usize]) { - let (masks, n_words) = match self.build_masks(indices) { - Some(m) => m, - None => return, - }; + let (masks, n_words) = masks_or_scalar_fallback!(self, indices, sqrt_x_dag_each); let zero = ::Store::zero(); self.data.iter_mut().for_each(|pw| { @@ -1647,4 +1705,184 @@ mod tests { assert_eq!(read(n + t), (false, true, false, true, 0)); } } + + // G-061: a repeated qubit index in a batched Clifford owes conjugation by + // `Gᵏ`, i.e. the per-index loop — `X 0 0` is legal Stim meaning + // apply-to-each-target-in-order, and the fused mask carries one bit per + // *distinct* site. `Gᵏ` is the identity only for the involutory gates + // (`S² = Z`, `(√X)² = X`, `(√Y)² = Y`), so neither one bit nor an + // XOR-cancelled bit is right for every family. + mod duplicate_index_tests { + use super::*; + use ppvm_pauli_sum::config::fxhash::ByteF64; + + // 2 u8 words → qubits 0..8 in word 0, qubits 8..16 in word 1. + type TC = ByteF64<2>; + type TTab = Tableau; + type GenTab = GeneralizedTableau; + + fn snapshot(tab: &TTab) -> Vec<(Vec, Vec, u8)> { + tab.data + .iter() + .map(|pw| { + ( + pw.word.xbits.data.as_raw_slice().to_vec(), + pw.word.zbits.data.as_raw_slice().to_vec(), + pw.phase, + ) + }) + .collect() + } + + /// A frame with all four bit patterns (`I`/`X`/`Y`/`Z`) live in the + /// rows, so a wrong number of gate applications cannot cancel out by + /// accident. + fn prepared(n: usize) -> TTab { + let mut tab = TTab::new(n); + for q in 0..n { + tab.h(q); + if q % 2 == 0 { + tab.s(q); + } + } + for q in 1..n { + tab.cnot(q - 1, q); + } + tab + } + + /// The batch on `$indices` must equal `$one` applied once per + /// occurrence. + macro_rules! check { + ($many:ident, $one:ident, $indices:expr) => {{ + let indices: &[usize] = $indices; + let mut batched = prepared(16); + let mut looped = prepared(16); + batched.$many(indices); + for &q in indices { + looped.$one(q); + } + assert_eq!( + snapshot(&batched), + snapshot(&looped), + "{} on repeated indices {:?} must equal the {} loop", + stringify!($many), + indices, + stringify!($one), + ); + }}; + } + + macro_rules! check_all { + ($indices:expr) => {{ + check!(x_many, x, $indices); + check!(y_many, y, $indices); + check!(z_many, z, $indices); + check!(h_many, h, $indices); + check!(s_many, s, $indices); + check!(s_dag_many, s_dag, $indices); + check!(sqrt_x_many, sqrt_x, $indices); + check!(sqrt_x_dag_many, sqrt_x_dag, $indices); + check!(sqrt_y_many, sqrt_y, $indices); + check!(sqrt_y_dag_many, sqrt_y_dag, $indices); + }}; + } + + #[test] + fn duplicate_indices_batch_equals_per_index_loop() { + check_all!(&[0, 0]); + check_all!(&[1, 1, 1]); + check_all!(&[0, 2, 2, 3]); + } + + /// Independent of the loop: `G²` is a Clifford in its own right, so + /// `g_many(&[q, q])` must equal conjugation by `G²`. + #[test] + fn duplicate_indices_equal_conjugation_by_g_squared() { + // Involutory: G² = I, so the batch must leave the frame untouched. + for (name, apply) in [ + ("x_many", TTab::x_many as fn(&mut TTab, &[usize])), + ("y_many", TTab::y_many), + ("z_many", TTab::z_many), + ("h_many", TTab::h_many), + ] { + let mut batched = prepared(16); + let untouched = prepared(16); + apply(&mut batched, &[1, 1]); + assert_eq!( + snapshot(&batched), + snapshot(&untouched), + "{name}(&[1, 1]) must be the identity (G² = I)" + ); + } + + // S² = Z, (S†)² = Z, (√X)² = X, (√X†)² = X, (√Y)² = Y, (√Y†)² = Y. + for (name, apply, square) in [ + ( + "s_many", + TTab::s_many as fn(&mut TTab, &[usize]), + TTab::z as fn(&mut TTab, usize), + ), + ("s_dag_many", TTab::s_dag_many, TTab::z), + ("sqrt_x_many", TTab::sqrt_x_many, TTab::x), + ("sqrt_x_dag_many", TTab::sqrt_x_dag_many, TTab::x), + ("sqrt_y_many", TTab::sqrt_y_many, TTab::y), + ("sqrt_y_dag_many", TTab::sqrt_y_dag_many, TTab::y), + ] { + let mut batched = prepared(16); + let mut squared = prepared(16); + apply(&mut batched, &[1, 1]); + square(&mut squared, 1); + assert_eq!( + snapshot(&batched), + snapshot(&squared), + "{name}(&[1, 1]) must equal conjugation by G²" + ); + } + } + + /// The same defect across storage words: the duplicate sits in word 1 + /// while another index sits in word 0. + #[test] + fn duplicate_indices_batch_equals_loop_multiword() { + check!(h_many, h, &[3, 10, 10]); + check!(s_many, s, &[3, 10, 10]); + check!(sqrt_y_many, sqrt_y, &[3, 10, 10]); + } + + /// The `GeneralizedTableau` wrappers filter lost qubits and forward the + /// rest, duplicates included, so they inherit the same contract. + #[test] + fn duplicate_indices_generalized_tableau_equals_loop() { + let indices: &[usize] = &[0, 2, 2]; + for lost in [None, Some(0usize)] { + let mut batched: GenTab = GeneralizedTableau::new(16, 1e-12); + let mut looped: GenTab = GeneralizedTableau::new(16, 1e-12); + for tab in [&mut batched, &mut looped] { + tab.tableau = prepared(16); + if let Some(q) = lost { + tab.is_lost[q] = true; + } + } + batched.sqrt_x_many(indices); + for &q in indices { + looped.sqrt_x(q); + } + assert_eq!( + snapshot(&batched.tableau), + snapshot(&looped.tableau), + "sqrt_x_many(&[0, 2, 2]) with lost = {lost:?}" + ); + } + } + + /// The user-visible symptom from G-061's `.stim` reproduction: + /// `X 0 0` on `|0…0⟩` is the identity, so the sampled bit must be 0. + #[test] + fn batched_x_twice_on_ground_state_measures_zero() { + let mut tab: GenTab = GeneralizedTableau::new(1, 1e-12); + tab.x_many(&[0, 0]); + assert_eq!(tab.measure(0), Some(false)); + } + } } diff --git a/crates/ppvm-tableau/src/noise.rs b/crates/ppvm-tableau/src/noise.rs index 8efb5e926..cd23017fd 100644 --- a/crates/ppvm-tableau/src/noise.rs +++ b/crates/ppvm-tableau/src/noise.rs @@ -56,6 +56,28 @@ where } } +/// `true` iff `p` is an admissible correlated-loss parameter triple. +/// +/// The channel is completely positive exactly when every event weight is +/// nonnegative, i.e. `p0, p1 >= 0`, `p0 + 2·p1 <= 1` (`p1` is the probability +/// that a *named* one of the pair is lost, so the exactly-one event carries +/// `2·p1`) and `p2 ∈ [0, 1]`. Outside that region the mixture truncates a +/// negative survivor weight and renormalizes, and the trajectory's cumulative +/// scan stops being a categorical sampler — both silently. +/// +/// Coefficients that do not convert to `f64` (symbolic ones) are not checked. +pub fn is_admissible_correlated_loss(p: &[C; 3]) -> bool { + // The slack matches `ppvm-tableau-2`'s guard: reject genuinely out-of-region + // parameters, not last-bit noise in a legitimately saturated triple. An exact + // `<= 1.0` would spuriously fire on e.g. `[1/3, 1/3, _]`, where + // `p0 + 2·p1` rounds to `1.0000000000000002`. + const SLACK: f64 = 1e-9; + let (Some(p0), Some(p1), Some(p2)) = (p[0].to_f64(), p[1].to_f64(), p[2].to_f64()) else { + return true; + }; + p0 >= 0.0 && p1 >= 0.0 && p0 + 2.0 * p1 <= 1.0 + SLACK && (0.0..=1.0).contains(&p2) +} + // === Noise trait impls === // // Orphan rules (E0210) forbid `impl> Depolarizing for X`, @@ -257,8 +279,12 @@ where /// The three probabilities are: /// * `p[0]`: The probability of losing both qubits simultaneously when /// both of them are in the qubit subspace. - /// * `p[1]`: The probability of losing either one qubit when both of them are - /// in the qubit subspace. + /// * `p[1]`: The probability of losing a **named** one of the two qubits when + /// both of them are in the qubit subspace, so the probability of losing + /// *exactly one* is `2·p[1]` and the both-present survivor keeps + /// `1 − 2·p[1] − p[0]` (which qubit is lost is 50/50). See + /// [`ppvm_traits::traits::CorrelatedLossChannel`] for the normative + /// statement. /// * `p[2]`: The probability of losing one qubit when the other one has already /// been lost prior to the channel. fn correlated_loss_channel( @@ -267,6 +293,10 @@ where addr1: usize, p: [::Coeff; 3], ) { + debug_assert!( + is_admissible_correlated_loss(&p), + "correlated loss needs p0, p1 >= 0, p0 + 2*p1 <= 1, p2 in [0, 1]; got {p:?}" + ); if self.is_lost[addr0] { self.loss_channel(addr1, p[2].clone()); return; @@ -278,7 +308,14 @@ where let r = self.tableau.rng.random::(); let mut cumulative = T::Coeff::zero(); for (i, p_i) in p[..2].iter().enumerate() { - cumulative += p_i.clone(); + // `p[1]` is the probability that a *named* one of the pair is lost, + // so the exactly-one event carries `2·p[1]` in this categorical + // scan; the fair coin below then picks which qubit. + cumulative += if i == 1 { + p_i.clone() + p_i.clone() + } else { + p_i.clone() + }; if cumulative > r { if i == 0 { // both lost @@ -731,12 +768,16 @@ mod tests { #[test] fn correlated_loss_p1_exactly_one_lost() { - // p[1]=1 → exactly one qubit lost each time. + // p[1]=0.5 → exactly one qubit lost each time, since `p[1]` is the + // probability that a *named* one of the pair is lost and so the + // exactly-one event carries 2·p[1] = 1. (`[0.0, 1.0, 0.0]`, which this + // test used before, is inadmissible under that convention: + // p0 + 2·p1 = 2 > 1.) let trials = 200; for seed in 0..trials { let mut t = tab(2); t.tableau.rng = rand::SeedableRng::seed_from_u64(seed); - t.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]); + t.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]); assert!( t.is_lost[0] ^ t.is_lost[1], "Expected exactly one lost qubit (seed {seed})" @@ -746,13 +787,14 @@ mod tests { #[test] fn correlated_loss_p1_both_qubits_chosen_equally() { - // With p[1]=1 the coin flip should lose addr0 and addr1 with equal frequency. + // With 2·p[1]=1 the coin flip should lose addr0 and addr1 with equal + // frequency. let trials = 1000u64; let mut addr0_lost = 0u64; for seed in 0..trials { let mut t = tab(2); t.tableau.rng = rand::SeedableRng::seed_from_u64(seed); - t.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]); + t.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]); if t.is_lost[0] { addr0_lost += 1; } @@ -790,7 +832,8 @@ mod tests { let mut t = tab(2); t.tableau.rng = rand::SeedableRng::seed_from_u64(seed); t.x(0); // put addr0 in |1⟩ - t.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]); + // 2·p[1] = 1: exactly one qubit is lost in every trial. + t.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]); if t.is_lost[0] { t.is_lost[0] = false; assert!(!t.measure(0).unwrap(), "Lost qubit should be reset to |0⟩"); @@ -853,8 +896,11 @@ mod tests { #[test] fn correlated_loss_statistics_single() { - // P(exactly one lost) should converge to p[1]. - let p_single = 0.4_f64; + // `p[1]` is the probability that a *named* one of the pair is lost, so + // P(exactly one lost) converges to 2·p[1]. (This test previously + // asserted `p[1]`, i.e. the rejected convention.) + let p_single = 0.2_f64; + let expected = 2.0 * p_single; let trials = 1000u64; let mut one_lost = 0u64; for seed in 0..trials { @@ -868,11 +914,132 @@ mod tests { let fraction = one_lost as f64 / trials as f64; // 5σ: σ = sqrt(0.4*0.6/1000) ≈ 0.015 assert!( - (fraction - p_single).abs() < 0.08, - "Expected ~{p_single:.2}, got {fraction:.3}" + (fraction - expected).abs() < 0.08, + "Expected ~{expected:.2}, got {fraction:.3}" ); } + // === G-040 — the correlated-loss `p[1]` convention === + // + // The paper (`ppvm-paper/main.tex:462`, `:523`, `:845`) is the definition of + // record: `p[1]` is `p_LQ`, the probability that a **named** one of the pair + // is lost, so `P(exactly one lost) = 2·p[1]` and the both-present survivor + // scales by `1 − 2·p[1] − p[0]`. `ppvm-pauli-sum` already reads it that way; + // these tests pin this trajectory sampler to the same number, because the + // cross-backend disagreement is the actual bug a Python user hits. + + type LossyTestSum = PauliSum< + ppvm_pauli_sum::config::fxhash::Byte< + 1, + f64, + NoStrategy, + LossyPauliWord<[u8; 1], fxhash::FxBuildHasher>, + >, + >; + + /// The Heisenberg scale factor `ppvm-pauli-sum` applies to a fully + /// in-subspace observable, i.e. `1 − p[0] − P(exactly one lost)`. + fn pauli_sum_survivor(p: [f64; 3]) -> f64 { + let mut sum = LossyTestSum::builder().n_qubits(2).build(); + sum += ("ZZ", 1.0); + sum.correlated_loss_channel(0, 1, p); + let zz: LossyPauliWord<[u8; 1], fxhash::FxBuildHasher> = "ZZ".into(); + *sum.data() + .get(&zz) + .expect("the all-present term survives a pure rescale") + } + + #[test] + fn correlated_loss_exactly_one_lost_is_two_p1_and_agrees_with_pauli_sum() { + let p1 = 0.3_f64; + let p = [0.0, p1, 0.0]; + let trials = 20_000u64; + let mut one_lost = 0u64; + for seed in 0..trials { + let mut t = tab(2); + t.tableau.rng = rand::SeedableRng::seed_from_u64(seed); + t.correlated_loss_channel(0, 1, p); + if t.is_lost[0] ^ t.is_lost[1] { + one_lost += 1; + } + } + let trajectory = one_lost as f64 / trials as f64; + // Same number, read off the (already correct) Heisenberg backend. + let pauli_sum = 1.0 - pauli_sum_survivor(p); + assert!( + (pauli_sum - 2.0 * p1).abs() < 1e-12, + "ppvm-pauli-sum P(exactly one) = {pauli_sum}, want 2*p[1] = {}", + 2.0 * p1 + ); + assert!( + (trajectory - 2.0 * p1).abs() < 0.02, + "trajectory P(exactly one lost) = {trajectory:.4}, want 2*p[1] = {} \ + (ppvm-pauli-sum says {pauli_sum})", + 2.0 * p1 + ); + } + + /// The paper's transport prediction: with `[p/3, p/3, p/3]` the per-event + /// both-present survival is `1 − p[0] − 2·p[1] = 1 − p`. + #[test] + fn correlated_loss_paper_transport_survival_is_one_minus_p() { + let p = 0.15_f64; + let third = p / 3.0; + let probabilities = [third, third, third]; + let trials = 20_000u64; + let mut none_lost = 0u64; + for seed in 0..trials { + let mut t = tab(2); + t.tableau.rng = rand::SeedableRng::seed_from_u64(seed); + t.correlated_loss_channel(0, 1, probabilities); + if !t.is_lost[0] && !t.is_lost[1] { + none_lost += 1; + } + } + let fraction = none_lost as f64 / trials as f64; + assert!( + (fraction - (1.0 - p)).abs() < 0.02, + "trajectory P(no loss) = {fraction:.4}, want 1 - {p} = {}", + 1.0 - p + ); + let survivor = pauli_sum_survivor(probabilities); + assert!( + (survivor - (1.0 - p)).abs() < 1e-12, + "ppvm-pauli-sum survivor {survivor}, want 1 - {p}" + ); + } + + // G-043 — the admissible region `p0, p1 >= 0`, `p0 + 2·p1 <= 1`, + // `p2 ∈ [0, 1]` is what makes the channel completely positive. Outside it + // the cumulative scan below silently stops being a categorical sampler. + + #[test] + #[cfg(debug_assertions)] + #[should_panic(expected = "p0 + 2*p1 <= 1")] + fn correlated_loss_rejects_inadmissible_probabilities() { + let mut t = tab(2); + t.correlated_loss_channel(0, 1, [0.6, 0.6, 0.0]); + } + + #[test] + #[cfg(debug_assertions)] + #[should_panic(expected = "p0 + 2*p1 <= 1")] + fn correlated_loss_rejects_negative_probabilities() { + let mut t = tab(2); + t.correlated_loss_channel(0, 1, [5.0, -3.0, 17.0]); + } + + /// The saturated boundary `p0 + 2·p1 == 1` is admissible and must not trip + /// the guard — including the `[1/3, 1/3, _]` triple whose sum rounds above + /// one in binary floating point. + #[test] + fn correlated_loss_saturated_boundary_is_admissible() { + let mut t = tab(2); + t.correlated_loss_channel(0, 1, [0.2, 0.4, 1.0]); + let mut t = tab(2); + t.correlated_loss_channel(0, 1, [1.0 / 3.0, 1.0 / 3.0, 0.5]); + } + // === z_expectation === #[test] diff --git a/crates/ppvm-traits/src/traits/noise.rs b/crates/ppvm-traits/src/traits/noise.rs index 80544ecb8..1dfaa6488 100644 --- a/crates/ppvm-traits/src/traits/noise.rs +++ b/crates/ppvm-traits/src/traits/noise.rs @@ -121,14 +121,33 @@ pub trait LossChannel { } /// Correlated two-qubit loss channel. +/// +/// # The `p[1]` convention (normative) +/// +/// This is the one place the parameterization is defined; every backend +/// (`ppvm-pauli-sum`, `ppvm-tableau`'s trajectory, `ppvm-tableau-sum`'s mixture) +/// and every binding cites it rather than restating it. In the paper's notation +/// `p = [p_LL, p_LQ, p_LN]`, and `p[1] = p_LQ` is the probability that a +/// **named** one of the two atoms is lost while the other survives. The two +/// single-loss events are disjoint, so +/// +/// * the probability of losing *exactly one* atom is `2·p[1]`, and +/// * the probability that both remain in the qubit subspace — the factor a +/// fully in-subspace observable is scaled by — is `1 − 2·p[1] − p[0]`. +/// +/// The channel is completely positive exactly on `p[0], p[1] >= 0`, +/// `p[0] + 2·p[1] <= 1`, `p[2] ∈ [0, 1]`; the tableau backends `debug_assert` +/// that region. pub trait CorrelatedLossChannel { /// Apply a correlated loss channel to qubits at `addr0` and `addr1`. /// /// The three probabilities are: /// * `p[0]`: The probability of losing both qubits simultaneously when /// both of them are in the qubit subspace. - /// * `p[1]`: The probability of losing either one qubit when both of them are - /// in the qubit subspace. + /// * `p[1]`: The probability of losing a **named** one of the two qubits when + /// both of them are in the qubit subspace, so losing *exactly one* has + /// probability `2·p[1]` and the both-present survivor is scaled by + /// `1 − 2·p[1] − p[0]` (which qubit is lost is 50/50). /// * `p[2]`: The probability of losing one qubit when the other one has already /// been lost prior to the channel. fn correlated_loss_channel(&mut self, addr0: usize, addr1: usize, p: [T::Coeff; 3]); diff --git a/ppvm-python/src/ppvm/mixins.py b/ppvm-python/src/ppvm/mixins.py index a181dc145..7080f9315 100644 --- a/ppvm-python/src/ppvm/mixins.py +++ b/ppvm-python/src/ppvm/mixins.py @@ -505,10 +505,18 @@ def correlated_loss_channel( - ``p[0]``: probability of losing both qubits simultaneously when both are in the qubit subspace. - - ``p[1]``: probability of losing exactly one qubit when both - are in the qubit subspace (which qubit is lost is 50/50 random). + - ``p[1]``: probability of losing a *named* one of the two + qubits when both are in the qubit subspace. The two + single-loss events are disjoint, so the probability of losing + **exactly one** qubit is ``2 * p[1]`` and the probability that + both remain is ``1 - 2 * p[1] - p[0]`` (which qubit is lost is + 50/50 random). - ``p[2]``: probability of losing the remaining active qubit when the other has already been lost prior to this channel. + + The channel is physical (completely positive) exactly on + ``p[0], p[1] >= 0``, ``p[0] + 2 * p[1] <= 1`` and + ``0 <= p[2] <= 1``. """ self._interface.correlated_loss_channel(addr0, addr1, p) diff --git a/ppvm-python/src/ppvm/paulisum.py b/ppvm-python/src/ppvm/paulisum.py index dca573d8d..67fe75a82 100644 --- a/ppvm-python/src/ppvm/paulisum.py +++ b/ppvm-python/src/ppvm/paulisum.py @@ -471,8 +471,12 @@ def correlated_loss_channel( The channel accepts 3 probabilities as argument: * `p[0]`: The probability of losing both qubits, when they are originally in the qubit subspace. - * `p[1]`: The probability of losing a single qubit, when both qubits - are originally in the qubit subspace. + * `p[1]`: The probability of losing a *named* one of the two qubits, + when both qubits are originally in the qubit subspace. The two + single-loss events are disjoint, so the probability of losing + **exactly one** qubit is `2 * p[1]`, and the probability that + both qubits remain in the qubit subspace — the factor a fully + in-subspace observable is scaled by — is `1 - 2 * p[1] - p[0]`. * `p[2]`: The probability of losing one qubit when the other one has already been lost prior to applying the channel. This is to account for the fact that when one qubit is missing during e.g. diff --git a/ppvm-python/test/generalized_tableau/test_loss.py b/ppvm-python/test/generalized_tableau/test_loss.py index c01f5de13..d8c5e3513 100644 --- a/ppvm-python/test/generalized_tableau/test_loss.py +++ b/ppvm-python/test/generalized_tableau/test_loss.py @@ -1,4 +1,4 @@ -from ppvm import GeneralizedTableau +from ppvm import GeneralizedTableau, LossyPauliSum from ppvm.generalized_tableau import MeasurementResult @@ -80,22 +80,25 @@ def test_correlated_loss_p0_both_lost(): def test_correlated_loss_p1_exactly_one_lost(): - # p[1]=1 → exactly one qubit lost in every trial. + # ``p[1]`` is the probability that a *named* one of the pair is lost, so + # P(exactly one lost) = 2 * p[1]. p[1] = 0.5 therefore loses exactly one + # qubit in every trial; p[1] = 1 would be inadmissible (p[0] + 2*p[1] > 1). for seed in range(200): tab = GeneralizedTableau(n_qubits=2, seed=seed) - tab.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]) + tab.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]) assert tab.is_lost(0) ^ tab.is_lost(1), f"Expected exactly one lost qubit (seed {seed})" def test_correlated_loss_p1_both_qubits_chosen_equally(): - # With p[1]=1 the 50/50 coin flip should lose addr0 and addr1 equally. + # With p[1]=0.5 (2*p[1] = 1) the 50/50 coin flip should lose addr0 and addr1 + # equally. trials = 1000 addr0_lost = sum( 1 for seed in range(trials) if ( tab := GeneralizedTableau(n_qubits=2, seed=seed), - tab.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]), + tab.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]), tab.is_lost(0), )[-1] ) @@ -118,11 +121,12 @@ def test_correlated_loss_both_lost_resets_to_zero(): def test_correlated_loss_single_lost_resets_to_zero(): - # The lost qubit should be reset to |0⟩. + # The lost qubit should be reset to |0⟩. p[1] = 0.5 (i.e. 2*p[1] = 1) is the + # admissible witness for "exactly one lost"; [0, 1, 0] is out of the region. for seed in range(1000): tab = GeneralizedTableau(n_qubits=2, seed=seed) tab.x(0) - tab.correlated_loss_channel(0, 1, [0.0, 1.0, 0.0]) + tab.correlated_loss_channel(0, 1, [0.0, 0.5, 0.0]) if tab.is_lost(0): tab.reset_loss_channel(0) assert tab.measure(0) == MeasurementResult.ZERO, "Lost qubit should be reset to |0⟩" @@ -171,8 +175,10 @@ def test_correlated_loss_statistics_both(): def test_correlated_loss_statistics_single(): - # P(exactly one lost) should converge to p[1]. + # ``p[1]`` is the probability that a *named* one of the pair is lost, so + # P(exactly one lost) should converge to 2 * p[1]. p_single = 0.4 + expected = 2 * p_single trials = 1000 one_lost = sum( 1 @@ -184,7 +190,39 @@ def test_correlated_loss_statistics_single(): )[-1] ) fraction = one_lost / trials - assert abs(fraction - p_single) < 0.08, f"Expected ~{p_single:.2f}, got {fraction:.3f}" + assert abs(fraction - expected) < 0.08, f"Expected ~{expected:.2f}, got {fraction:.3f}" + + +def test_correlated_loss_exactly_one_lost_agrees_with_lossy_pauli_sum(): + # ``p[1]`` is the probability that a *named* one of the pair is lost, so + # P(exactly one lost) = ``2 * p[1]`` on *every* backend. Reading the same + # number off ``LossyPauliSum`` and off ``GeneralizedTableau`` used to give + # answers that differed by a factor of two — that cross-backend split is the + # actual bug this test pins. + p1 = 0.3 + p = [0.0, p1, 0.0] + + # Heisenberg picture: a fully in-subspace observable is scaled by + # ``1 - 2 * p[1] - p[0]``, so its complement is P(exactly one lost). + ps = LossyPauliSum.new(n_qubits=2, terms=["ZZ"]) + ps.correlated_loss_channel(0, 1, p) + pauli_sum = 1.0 - ps.overlap_with_zero() + assert abs(pauli_sum - 2 * p1) < 1e-9, ( + f"LossyPauliSum P(exactly one) = {pauli_sum}, want 2 * p[1] = {2 * p1}" + ) + + # Trajectory picture: the sampled frequency of the exactly-one-lost outcome. + trials = 4000 + one_lost = 0 + for seed in range(trials): + tab = GeneralizedTableau(n_qubits=2, seed=seed) + tab.correlated_loss_channel(0, 1, p) + one_lost += tab.is_lost(0) ^ tab.is_lost(1) + trajectory = one_lost / trials + assert abs(trajectory - 2 * p1) < 0.04, ( + f"GeneralizedTableau P(exactly one lost) = {trajectory:.4f}, " + f"want 2 * p[1] = {2 * p1} (LossyPauliSum says {pauli_sum})" + ) # === AsymmetricLossChannel === diff --git a/skills/ppvm-usage/SKILL.md b/skills/ppvm-usage/SKILL.md index 9d7180edc..68e84edae 100644 --- a/skills/ppvm-usage/SKILL.md +++ b/skills/ppvm-usage/SKILL.md @@ -110,7 +110,7 @@ for _ in range(50): print(ps.overlap_with_zero()) ``` -Loss channels live on `LossyPauliSum` (same API, plus `loss_channel(q, p)` and `correlated_loss_channel(q0, q1, [p_x, p_y, p_z])`). +Loss channels live on `LossyPauliSum` (same API, plus `loss_channel(q, p)` and `correlated_loss_channel(q0, q1, [p_LL, p_LQ, p_LN])`). ### Generalized stabilizer tableau @@ -266,7 +266,7 @@ Important: the six off-diagonal two-qubit rotations (`rxy`, `rxz`, `ryx`, `ryz`, | `two_qubit_pauli_error(q0, q1, p[15])` | ✓ | ✓ | ✓ | | `amplitude_damping(q, gamma)` | ✓ | ✓ | — | | `loss_channel(q, p)` (Lossy types) | ✓ | ✓\* | ✓ | -| `correlated_loss_channel(q0, q1, [px,py,pz])` | ✓ | ✓\* | ✓ | +| `correlated_loss_channel(q0, q1, [p_LL, p_LQ, p_LN])` | ✓ | ✓\* | ✓ | | `reset_loss_channel(q)` | ✓ | ✓\* | ✓ | \* Python side: loss methods live on `LossyPauliSum`, not the plain `PauliSum`.