From 305da3c170cf453d44d18debc9a12df97557ba60 Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sat, 11 Apr 2026 23:52:22 -0400 Subject: [PATCH 1/8] chore(autotune): set up rayon threading experiment with baselines Co-Authored-By: Claude Opus 4.6 (1M context) --- crates/ppvm-tableau/examples/profile_msd.rs | 149 ++++++++++++++++++ .../ppvm-tableau/examples/profile_scaling.rs | 109 +++++++++++++ .../2026-04-11-tableau-rayon-threading/log.md | 31 ++++ .../metric.toml | 9 ++ 4 files changed, 298 insertions(+) create mode 100644 crates/ppvm-tableau/examples/profile_msd.rs create mode 100644 crates/ppvm-tableau/examples/profile_scaling.rs create mode 100644 docs/autotune/2026-04-11-tableau-rayon-threading/log.md create mode 100644 docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml diff --git a/crates/ppvm-tableau/examples/profile_msd.rs b/crates/ppvm-tableau/examples/profile_msd.rs new file mode 100644 index 000000000..6a17201d8 --- /dev/null +++ b/crates/ppvm-tableau/examples/profile_msd.rs @@ -0,0 +1,149 @@ +use std::time::Instant; + +use ppvm_runtime::config::fx64hash::Byte8F64; +use ppvm_tableau::prelude::*; + +type Tab = GeneralizedTableau, u128>; + +fn encode(tab: &mut Tab, qubits: &[usize]) { + if qubits.len() == 17 { + for i in [0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16] { + tab.sqrt_y(qubits[i]); + } + for [i, j] in [[1, 3], [7, 10], [12, 14], [13, 16]] { + tab.cz(qubits[i], qubits[j]); + } + for i in [7, 16] { + tab.sqrt_y_adj(qubits[i]); + } + for [i, j] in [[4, 7], [8, 10], [11, 14], [15, 16]] { + tab.cz(qubits[i], qubits[j]); + } + for i in [4, 10, 14, 16] { + tab.sqrt_y_adj(qubits[i]); + } + for [i, j] in [[2, 4], [6, 8], [7, 9], [10, 13], [14, 16]] { + tab.cz(qubits[i], qubits[j]); + } + for i in [3, 6, 9, 10, 12, 13] { + tab.sqrt_y(qubits[i]); + } + for [i, j] in [[0, 2], [3, 6], [5, 8], [10, 12], [11, 13]] { + tab.cz(qubits[i], qubits[j]); + } + for i in [1, 2, 3, 4, 6, 7, 8, 9, 11, 12, 14] { + tab.sqrt_y(qubits[i]); + } + for [i, j] in [[0, 1], [2, 3], [4, 5], [6, 7], [8, 9], [12, 15]] { + tab.cz(qubits[i], qubits[j]); + } + for i in [0, 2, 5, 6, 8, 10, 12] { + tab.sqrt_y_adj(qubits[i]); + } + } +} + +fn main() { + let n_qubits = 85; + let n_runs = 20; + + let mut encoding_total = 0u128; + let mut middle_gates_total = 0u128; + let mut measure_total = 0u128; + let mut full_total = 0u128; + + for _ in 0..n_runs { + let full_start = Instant::now(); + + let mut tab: Tab = GeneralizedTableau::new(n_qubits, 1e-10); + let qubit_addrs: Vec = (0..n_qubits).collect(); + let ql: Vec<&[usize]> = qubit_addrs.chunks_exact(17).collect(); + + // Phase 1: Encoding (H + T + encode per block) + let t0 = Instant::now(); + for q in ql.iter() { + let encoding_qubit = q[7]; + tab.h(encoding_qubit); + tab.t(encoding_qubit); + encode(&mut tab, q); + } + encoding_total += t0.elapsed().as_nanos(); + + // Phase 2: Middle gates (sqrt_x, cz, sqrt_y, sqrt_x_adj layers) + let t1 = Instant::now(); + for i in [0, 1, 4] { + for q in ql[i] { + tab.sqrt_x(*q); + } + } + for (control, target) in ql[0].iter().zip(ql[1]) { + tab.cz(*control, *target); + } + for (control, target) in ql[2].iter().zip(ql[3]) { + tab.cz(*control, *target); + } + for q in ql[0] { + tab.sqrt_y(*q); + } + for q in ql[3] { + tab.sqrt_y(*q); + } + for (control, target) in ql[0].iter().zip(ql[2]) { + tab.cz(*control, *target); + } + for (control, target) in ql[3].iter().zip(ql[4]) { + tab.cz(*control, *target); + } + for q in ql[0] { + tab.sqrt_x_adj(*q); + } + for (control, target) in ql[0].iter().zip(ql[4]) { + tab.cz(*control, *target); + } + for (control, target) in ql[1].iter().zip(ql[3]) { + tab.cz(*control, *target); + } + for i in 0..5 { + for q in ql[i] { + tab.sqrt_x_adj(*q); + } + } + middle_gates_total += t1.elapsed().as_nanos(); + + // Phase 3: Measurement + let t2 = Instant::now(); + for i in 0..n_qubits { + tab.measure(i); + } + measure_total += t2.elapsed().as_nanos(); + + full_total += full_start.elapsed().as_nanos(); + } + + let div = n_runs as f64; + println!("=== MSD 85-qubit profile ({n_runs} runs) ==="); + println!( + "Encoding (H+T+encode): {:.1} µs", + encoding_total as f64 / div / 1000.0 + ); + println!( + "Middle gates: {:.1} µs", + middle_gates_total as f64 / div / 1000.0 + ); + println!( + "Measurement: {:.1} µs", + measure_total as f64 / div / 1000.0 + ); + println!( + "Full run: {:.1} µs", + full_total as f64 / div / 1000.0 + ); + println!(); + let total = encoding_total as f64; + let mid = middle_gates_total as f64; + let meas = measure_total as f64; + let full = full_total as f64; + println!("Encoding: {:.1}%", total / full * 100.0); + println!("Middle: {:.1}%", mid / full * 100.0); + println!("Measurement: {:.1}%", meas / full * 100.0); +} diff --git a/crates/ppvm-tableau/examples/profile_scaling.rs b/crates/ppvm-tableau/examples/profile_scaling.rs new file mode 100644 index 000000000..3b8beb91a --- /dev/null +++ b/crates/ppvm-tableau/examples/profile_scaling.rs @@ -0,0 +1,109 @@ +use std::time::Instant; + +use ppvm_runtime::config::fx64hash::Byte8F64; +use ppvm_tableau::prelude::*; + +type Tab = GeneralizedTableau, usize>; + +fn main() { + let n_runs = 10; + + // Profile 1: Scaling with qubit count (Clifford-heavy) + println!("=== Clifford gate scaling (CNOT chain) ==="); + for n_qubits in [32, 64, 96, 128] { + let mut total_ns = 0u128; + for _ in 0..n_runs { + let mut tab: Tab = GeneralizedTableau::new(n_qubits, 1e-10); + tab.h(0); + let t0 = Instant::now(); + for i in 0..n_qubits - 1 { + tab.cnot(i, i + 1); + } + total_ns += t0.elapsed().as_nanos(); + } + println!( + " n={:3}: {:.1} µs ({} cnots, {} rows)", + n_qubits, + total_ns as f64 / n_runs as f64 / 1000.0, + n_qubits - 1, + 2 * n_qubits + ); + } + + // Profile 2: T-gate coefficient branching scaling + println!("\n=== T-gate branching scaling (85 qubits) ==="); + for n_tgates in [1, 4, 8, 12, 16] { + let mut tab: Tab = GeneralizedTableau::new(85, 1e-10); + for i in 0..n_tgates.min(85) { + tab.h(i); // ensure non-deterministic for branching + } + + let mut total_ns = 0u128; + let mut coeff_count = 0; + for _ in 0..n_runs { + let mut t = tab.fork(Some(42)); + let t0 = Instant::now(); + for i in 0..n_tgates { + t.t(i % 85); + } + total_ns += t0.elapsed().as_nanos(); + coeff_count = t.coefficients.len(); + } + println!( + " t={:2}: {:.1} µs (coeffs={})", + n_tgates, + total_ns as f64 / n_runs as f64 / 1000.0, + coeff_count, + ); + } + + // Profile 3: Measurement after T gates + println!("\n=== Measurement scaling after T gates (85 qubits) ==="); + for n_tgates in [1, 4, 8, 12] { + let mut tab: Tab = GeneralizedTableau::new(85, 1e-10); + for i in 0..n_tgates.min(85) { + tab.h(i); + } + for i in 0..n_tgates { + tab.t(i % 85); + } + + let mut total_ns = 0u128; + for _ in 0..n_runs { + let mut t = tab.fork(Some(42)); + let t0 = Instant::now(); + for i in 0..85 { + t.measure(i); + } + total_ns += t0.elapsed().as_nanos(); + } + println!( + " t={:2}: {:.1} µs (85 measures, initial_coeffs={})", + n_tgates, + total_ns as f64 / n_runs as f64 / 1000.0, + tab.coefficients.len(), + ); + } + + // Profile 4: Individual gate cost at different qubit counts + println!("\n=== Single gate cost (per-call overhead test) ==="); + for n_qubits in [32, 64, 128] { + let tab: Tab = GeneralizedTableau::new(n_qubits, 1e-10); + let n_calls = 1000; + let mut total_ns = 0u128; + for _ in 0..n_runs { + let mut t = tab.fork(Some(42)); + let t0 = Instant::now(); + for _ in 0..n_calls { + t.h(0); + } + total_ns += t0.elapsed().as_nanos(); + } + println!( + " n={:3}: {:.0} ns/call (H gate, {} rows)", + n_qubits, + total_ns as f64 / (n_runs * n_calls) as f64, + 2 * n_qubits + ); + } +} diff --git a/docs/autotune/2026-04-11-tableau-rayon-threading/log.md b/docs/autotune/2026-04-11-tableau-rayon-threading/log.md new file mode 100644 index 000000000..95557ecf0 --- /dev/null +++ b/docs/autotune/2026-04-11-tableau-rayon-threading/log.md @@ -0,0 +1,31 @@ +# Log for 2026-04-11-tableau-rayon-threading + +## 2026-04-11 + +### Architecture Notes + +**Profiling breakdown (85-qubit MSD circuit, ~135µs total):** +- Encoding (H+T+encode): 55.5µs (41.1%) — includes 5 T gates that branch coefficients +- Middle gates (Clifford): 38.7µs (28.7%) — all Clifford operations +- Measurement: 39.3µs (29.1%) — 85 sequential measurements + +**T-gate branching scaling (85 qubits):** +- 1 T gate: 0.2µs (2 coefficients) +- 4 T gates: 2.1µs (16 coefficients) +- 8 T gates: 10.1µs (256 coefficients) +- 12 T gates: 105.1µs (4096 coefficients) +- 16 T gates: 1569.5µs (65536 coefficients) + +**Single Clifford gate cost:** +- 32 qubits (64 rows): 43 ns/call +- 64 qubits (128 rows): 82 ns/call +- 128 qubits (256 rows): 161 ns/call + +**Key insight:** Parallelizing individual Clifford gate row iteration is NOT viable — per-call overhead (~1µs rayon) far exceeds the gate cost (43-161ns). The primary parallelism target is `branch_with_coefficients` which dominates at high coefficient counts (exponential growth with T gates). + +**Parallelism targets (by expected impact):** +1. `branch_with_coefficients` (data.rs) — coefficient loop, O(M) where M grows exponentially with T gates. Each coefficient independently computes phase+branch, but accumulates into shared HashMap. Strategy: parallel fold with local HashMaps, then merge. +2. Measurement overlap computation (measure.rs) — O(M) iteration computing z_overlap. +3. Measurement coefficient update (measure.rs) — processing b_keys entries. + +**Design constraint:** All parallelism behind `features = ["rayon"]`. When off, zero regression — identical codepath. diff --git a/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml b/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml new file mode 100644 index 000000000..51072dd39 --- /dev/null +++ b/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml @@ -0,0 +1,9 @@ +[[metric]] +"commit" = "5da53ba128809b6598a477b0cd8207b72216014b" +"status" = "keep" +"description" = "baseline (no rayon)" +"t_gate_1_ns" = 192.0 +"t_gate_4_us" = 1.07 +"t_gate_8_us" = 8.25 +"t_gate_12_us" = 50.87 +"msd_us" = 124.97 From b574c45828cc59daa0c00a7084f911536a31013a Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 14:52:01 -0400 Subject: [PATCH 2/8] feat(tableau): add rayon feature for parallel coefficient branching Adds an optional `rayon` feature to ppvm-tableau that parallelizes the coefficient accumulation in `branch_with_coefficients` and `compute_coefficients_after_pauli_apply` using rayon's parallel fold/reduce. Key design decisions: - Feature-gated: all parallel code behind `features = ["rayon"]` - Threshold-based: only uses rayon when coefficient count >= 4096, below that sequential is always faster - Zero regression: when feature is off, code is identical to before (refactored to use extracted free functions but same logic) - Same API: no interface changes, just enable the feature Also adds Send + Sync + Copy + num::Num bounds to impl blocks across the crate (tgate, rot1, rot2, measure, stim, noise, reset) - all concrete types (f64, usize, u128, BUint) satisfy these trivially. Co-Authored-By: Claude Opus 4.6 (1M context) --- Cargo.lock | 1 + crates/ppvm-tableau/Cargo.toml | 5 + crates/ppvm-tableau/src/data.rs | 328 ++++++++++++++++++++----- crates/ppvm-tableau/src/gates/reset.rs | 9 +- crates/ppvm-tableau/src/gates/rot1.rs | 7 +- crates/ppvm-tableau/src/gates/rot2.rs | 7 +- crates/ppvm-tableau/src/gates/tgate.rs | 7 +- crates/ppvm-tableau/src/measure.rs | 9 +- crates/ppvm-tableau/src/noise.rs | 32 ++- crates/ppvm-tableau/src/stim.rs | 9 +- 10 files changed, 328 insertions(+), 86 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index ec43e798c..b1a7d9e8c 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -806,6 +806,7 @@ dependencies = [ "num", "ppvm-runtime", "rand", + "rayon", ] [[package]] diff --git a/crates/ppvm-tableau/Cargo.toml b/crates/ppvm-tableau/Cargo.toml index f1499e937..d2bd1ad4c 100644 --- a/crates/ppvm-tableau/Cargo.toml +++ b/crates/ppvm-tableau/Cargo.toml @@ -3,6 +3,10 @@ name = "ppvm-tableau" version = "0.1.0" edition = "2024" +[features] +default = [] +rayon = ["dep:rayon"] + [dependencies] bitvec = "1.0.1" fxhash = "0.2.1" @@ -11,6 +15,7 @@ itertools = "0.14.0" num = "0.4.3" ppvm-runtime = { version = "0.1.0", path = "../ppvm-runtime" } rand = "0.10.0" +rayon = { version = "1.10", optional = true } [dev-dependencies] criterion = "0.7.0" diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index 268ae2df0..fb6ddee20 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -277,6 +277,223 @@ where (alpha & beta).count_ones() } +/// Compute the phase contribution from destabilizer anticommutation and odd-phase masks. +/// This is a pure function of its arguments (no self access needed), extracted to enable +/// use in parallel contexts where borrowing self is not possible. +#[inline] +fn compute_phase_with_mask_static( + destab_anticomm_bits: I, + basis_index: I, + stab_anticomm_bits: I, + odd_phase_mask: I, +) -> u8 { + let mut phase = (2 * symplectic_inner(destab_anticomm_bits, basis_index) as u8) % 4; + let active = basis_index & stab_anticomm_bits; + let parity = (active & odd_phase_mask).count_ones() % 2; + phase = (phase + 2 * parity as u8) % 4; + phase +} + +/// Minimum number of coefficients before engaging rayon parallelism. +/// Below this threshold, the sequential path is always used even with the `rayon` feature, +/// because rayon's thread-pool overhead (~20µs) dominates the per-coefficient work. +#[cfg(feature = "rayon")] +const RAYON_COEFF_THRESHOLD: usize = 4096; + +/// Sequential accumulation of branch coefficients. +fn branch_coefficients_seq( + items: impl IntoIterator, I)>, + stab_anticomm_bits: I, + destab_anticomm_bits: I, + odd_phase_mask: I, + phase_decomp: u8, + coefficient_factor: Complex, + branch_factor: Complex, +) -> HashMap> +where + I: TableauIndex, + CoeffType: One + Zero + Clone + num::Num, + Complex: + std::ops::Mul> + std::ops::AddAssign + From + Copy, +{ + let mut map: HashMap> = HashMap::default(); + for (coeff, idx) in items { + debug_assert!( + !(coeff.re == CoeffType::zero() && coeff.im == CoeffType::zero()), + "Coefficient should not be zero" + ); + let branch_index = idx ^ stab_anticomm_bits; + let branch_phase_contribution = compute_phase_with_mask_static( + destab_anticomm_bits, + idx, + stab_anticomm_bits, + odd_phase_mask, + ); + let branch_phase = (branch_phase_contribution + phase_decomp) % 4; + let phase_factor: Complex = + COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); + let branch_coefficient = phase_factor * coeff * branch_factor; + let nonbranch_coefficient = coeff * coefficient_factor; + *map.entry(branch_index).or_insert(Complex::zero()) += branch_coefficient; + *map.entry(idx).or_insert(Complex::zero()) += nonbranch_coefficient; + } + map +} + +/// Accumulate branch coefficients. Uses rayon parallel fold/reduce when the feature +/// is enabled AND the coefficient count exceeds `RAYON_COEFF_THRESHOLD`. +#[cfg(feature = "rayon")] +fn branch_coefficients_parallel( + items: &[(Complex, I)], + stab_anticomm_bits: I, + destab_anticomm_bits: I, + odd_phase_mask: I, + phase_decomp: u8, + coefficient_factor: Complex, + branch_factor: Complex, +) -> HashMap> +where + I: TableauIndex + Send + Sync, + CoeffType: One + Zero + Clone + Send + Sync + num::Num, + Complex: + std::ops::Mul> + std::ops::AddAssign + From + Copy, +{ + #[cfg(feature = "rayon")] + if items.len() >= RAYON_COEFF_THRESHOLD { + use rayon::prelude::*; + return items + .par_iter() + .fold( + HashMap::>::default, + |mut map, &(coeff, idx)| { + let branch_index = idx ^ stab_anticomm_bits; + let branch_phase_contribution = compute_phase_with_mask_static( + destab_anticomm_bits, + idx, + stab_anticomm_bits, + odd_phase_mask, + ); + let branch_phase = (branch_phase_contribution + phase_decomp) % 4; + let phase_factor: Complex = + COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); + let branch_coefficient = phase_factor * coeff * branch_factor; + let nonbranch_coefficient = coeff * coefficient_factor; + *map.entry(branch_index).or_insert(Complex::zero()) += branch_coefficient; + *map.entry(idx).or_insert(Complex::zero()) += nonbranch_coefficient; + map + }, + ) + .reduce(HashMap::default, |mut a, b| { + for (k, v) in b { + *a.entry(k).or_insert(Complex::zero()) += v; + } + a + }); + } + + branch_coefficients_seq( + items.iter().copied(), + stab_anticomm_bits, + destab_anticomm_bits, + odd_phase_mask, + phase_decomp, + coefficient_factor, + branch_factor, + ) +} + +/// Sequential accumulation of apply coefficients. +fn apply_coefficients_seq( + items: impl IntoIterator, I)>, + stab_anticomm_bits: I, + destab_anticomm_bits: I, + odd_phase_mask: I, + phase_decomp: u8, +) -> HashMap> +where + I: TableauIndex, + CoeffType: One + Zero + Clone + num::Num, + Complex: + std::ops::Mul> + std::ops::AddAssign + From + Copy, +{ + let mut map: HashMap> = HashMap::default(); + for (coeff, idx) in items { + debug_assert!( + !(coeff.re == CoeffType::zero() && coeff.im == CoeffType::zero()), + "Coefficient should not be zero" + ); + let branch_index = idx ^ stab_anticomm_bits; + let branch_phase_contribution = compute_phase_with_mask_static( + destab_anticomm_bits, + idx, + stab_anticomm_bits, + odd_phase_mask, + ); + let branch_phase = (branch_phase_contribution + phase_decomp) % 4; + let phase_factor: Complex = + COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); + let branch_coefficient = phase_factor * coeff; + *map.entry(branch_index).or_insert(Complex::zero()) += branch_coefficient; + } + map +} + +/// Accumulate coefficients for pauli application. Uses rayon when the feature is +/// enabled AND the coefficient count exceeds `RAYON_COEFF_THRESHOLD`. +#[cfg(feature = "rayon")] +fn apply_coefficients_parallel( + items: &[(Complex, I)], + stab_anticomm_bits: I, + destab_anticomm_bits: I, + odd_phase_mask: I, + phase_decomp: u8, +) -> HashMap> +where + I: TableauIndex + Send + Sync, + CoeffType: One + Zero + Clone + Send + Sync + num::Num, + Complex: + std::ops::Mul> + std::ops::AddAssign + From + Copy, +{ + #[cfg(feature = "rayon")] + if items.len() >= RAYON_COEFF_THRESHOLD { + use rayon::prelude::*; + return items + .par_iter() + .fold( + HashMap::>::default, + |mut map, &(coeff, idx)| { + let branch_index = idx ^ stab_anticomm_bits; + let branch_phase_contribution = compute_phase_with_mask_static( + destab_anticomm_bits, + idx, + stab_anticomm_bits, + odd_phase_mask, + ); + let branch_phase = (branch_phase_contribution + phase_decomp) % 4; + let phase_factor: Complex = + COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); + let branch_coefficient = phase_factor * coeff; + *map.entry(branch_index).or_insert(Complex::zero()) += branch_coefficient; + map + }, + ) + .reduce(HashMap::default, |mut a, b| { + for (k, v) in b { + *a.entry(k).or_insert(Complex::zero()) += v; + } + a + }); + } + + apply_coefficients_seq( + items.iter().copied(), + stab_anticomm_bits, + destab_anticomm_bits, + odd_phase_mask, + phase_decomp, + ) +} + // TODO: builder #[derive(Clone)] pub struct GeneralizedTableau< @@ -293,12 +510,13 @@ pub struct GeneralizedTableau< impl, I>> GeneralizedTableau where - T::Coeff: One + Zero + Clone, + T::Coeff: One + Zero + Clone + Send + Sync + num::Num, Complex: std::ops::Mul> + std::ops::AddAssign + From - + ComplexFloat, - I: TableauIndex, + + ComplexFloat + + Copy, + I: TableauIndex + Send + Sync, { pub fn new(n_qubits: usize, coefficient_threshold: T::Coeff) -> Self { let mut coefficients = C::new(); @@ -527,11 +745,12 @@ where stab_anticomm_bits: I, odd_phase_mask: I, ) -> u8 { - let mut phase = (2 * symplectic_inner(destab_anticomm_bits, basis_index) as u8) % 4; - let active = basis_index & stab_anticomm_bits; - let parity = (active & odd_phase_mask).count_ones() % 2; - phase = (phase + 2 * parity as u8) % 4; - phase + compute_phase_with_mask_static( + destab_anticomm_bits, + basis_index, + stab_anticomm_bits, + odd_phase_mask, + ) } /// Keep only coefficients that correspond to the correct eigenvalue of a Z measurement. @@ -585,38 +804,32 @@ where let odd_phase_mask = self.odd_phase_destabilizer_mask(); let old_coefficients = std::mem::replace(&mut self.coefficients, C::new()); - let mut new_coefficients: HashMap> = HashMap::default(); - for (coeff, idx) in old_coefficients.into_iter() { - debug_assert!( - !(coeff.re == T::Coeff::zero() && coeff.im == T::Coeff::zero()), - "Coefficient should not be zero" - ); - - let branch_index = idx ^ stab_anticomm_bits; // stab_anticomm_bits is the index shift - - // get the phase contributions from duplicate destabilizers - // and anti-commuting through destabilizers - let branch_phase_contribution = self.compute_phase_with_mask( - destab_anticomm_bits, - idx, + + // Only collect to Vec when rayon might be used; otherwise iterate directly + #[cfg(feature = "rayon")] + let new_coefficients = { + let items: Vec<_> = old_coefficients.into_iter().collect(); + branch_coefficients_parallel( + &items, stab_anticomm_bits, + destab_anticomm_bits, odd_phase_mask, - ); - - // the total phase is the product of the above with the decomposition phase - let branch_phase = (branch_phase_contribution + phase_decomp) % 4; - - let phase_factor: Complex = - COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); - - let branch_coefficient = phase_factor * coeff * branch_factor; - let nonbranch_coefficient = coeff * coefficient_factor; + phase_decomp, + coefficient_factor, + branch_factor, + ) + }; - *new_coefficients - .entry(branch_index) - .or_insert(Complex::zero()) += branch_coefficient; - *new_coefficients.entry(idx).or_insert(Complex::zero()) += nonbranch_coefficient; - } + #[cfg(not(feature = "rayon"))] + let new_coefficients = branch_coefficients_seq( + old_coefficients, + stab_anticomm_bits, + destab_anticomm_bits, + odd_phase_mask, + phase_decomp, + coefficient_factor, + branch_factor, + ); let cutoff = Complex { re: self.coefficient_threshold.clone(), @@ -646,35 +859,28 @@ where self.compute_decomposition(addr0, pauli); let odd_phase_mask = self.odd_phase_destabilizer_mask(); - let mut new_coefficients: HashMap> = HashMap::default(); let old_coefficients = std::mem::replace(coefficients, C::new()); - for (coeff, idx) in old_coefficients.into_iter() { - debug_assert!( - !(coeff.re == T::Coeff::zero() && coeff.im == T::Coeff::zero()), - "Coefficient should not be zero" - ); - - let branch_index = idx ^ stab_anticomm_bits; // stab_anticomm_bits is the index shift - // get the phase contributions from duplicate destabilizers - // and anti-commuting through destabilizers - let branch_phase_contribution = self.compute_phase_with_mask( - destab_anticomm_bits, - idx, + #[cfg(feature = "rayon")] + let new_coefficients = { + let items: Vec<_> = old_coefficients.into_iter().collect(); + apply_coefficients_parallel( + &items, stab_anticomm_bits, + destab_anticomm_bits, odd_phase_mask, - ); - let branch_phase = (branch_phase_contribution + phase_decomp) % 4; - - let phase_factor: Complex = - COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); - - let branch_coefficient = phase_factor * coeff; + phase_decomp, + ) + }; - *new_coefficients - .entry(branch_index) - .or_insert(Complex::zero()) += branch_coefficient; - } + #[cfg(not(feature = "rayon"))] + let new_coefficients = apply_coefficients_seq( + old_coefficients, + stab_anticomm_bits, + destab_anticomm_bits, + odd_phase_mask, + phase_decomp, + ); let cutoff = Complex { re: self.coefficient_threshold.clone(), diff --git a/crates/ppvm-tableau/src/gates/reset.rs b/crates/ppvm-tableau/src/gates/reset.rs index b1d1abde0..754535604 100644 --- a/crates/ppvm-tableau/src/gates/reset.rs +++ b/crates/ppvm-tableau/src/gates/reset.rs @@ -21,7 +21,7 @@ impl Reset for GeneralizedTableau where T: Config, <::Storage as BitView>::Store: PrimInt, - I: TableauIndex + Debug, + I: TableauIndex + Debug + Send + Sync, C: SparseVector, I> + Debug, T::Coeff: One + Zero @@ -30,13 +30,16 @@ where + ToPrimitive + std::fmt::Debug + std::ops::Mul - + PartialOrd, + + PartialOrd + + Send + + Sync, Complex: std::ops::Mul> + From + std::ops::MulAssign + std::ops::AddAssign + One - + ComplexFloat, + + ComplexFloat + + Copy, { fn reset(&mut self, addr0: usize) { let m = self.measure(addr0); diff --git a/crates/ppvm-tableau/src/gates/rot1.rs b/crates/ppvm-tableau/src/gates/rot1.rs index 33bed76f4..f28a80fba 100644 --- a/crates/ppvm-tableau/src/gates/rot1.rs +++ b/crates/ppvm-tableau/src/gates/rot1.rs @@ -11,12 +11,13 @@ impl, I>> RotationOne for GeneralizedTableau where <::Storage as BitView>::Store: PrimInt, - T::Coeff: Zero + One, - I: TableauIndex, + T::Coeff: Zero + One + Send + Sync + num::Num, + I: TableauIndex + Send + Sync, Complex: std::ops::Mul> + std::ops::AddAssign + From - + ComplexFloat, + + ComplexFloat + + Copy, { fn rotate_1(&mut self, axis: Pauli, addr0: usize, theta: ::Coeff) { if self.is_lost[addr0] { diff --git a/crates/ppvm-tableau/src/gates/rot2.rs b/crates/ppvm-tableau/src/gates/rot2.rs index c108d7d20..d87c463f1 100644 --- a/crates/ppvm-tableau/src/gates/rot2.rs +++ b/crates/ppvm-tableau/src/gates/rot2.rs @@ -15,12 +15,13 @@ impl, I>> RotationTwo for GeneralizedTableau where <::Storage as BitView>::Store: PrimInt, - T::Coeff: Zero + One, - I: TableauIndex, + T::Coeff: Zero + One + Send + Sync + num::Num, + I: TableauIndex + Send + Sync, Complex: std::ops::Mul> + std::ops::AddAssign + From - + ComplexFloat, + + ComplexFloat + + Copy, { fn rotate_2( &mut self, diff --git a/crates/ppvm-tableau/src/gates/tgate.rs b/crates/ppvm-tableau/src/gates/tgate.rs index db089a232..dee72d626 100644 --- a/crates/ppvm-tableau/src/gates/tgate.rs +++ b/crates/ppvm-tableau/src/gates/tgate.rs @@ -19,12 +19,13 @@ where T: Config, <::Storage as BitView>::Store: PrimInt, C: SparseVector, I>, - T::Coeff: One + Zero + Clone, + T::Coeff: One + Zero + Clone + Send + Sync + num::Num, Complex: std::ops::Mul> + std::ops::AddAssign + From - + ComplexFloat, - I: TableauIndex, + + ComplexFloat + + Copy, + I: TableauIndex + Send + Sync, >::Output>>::Output: PartialEq, { fn t(&mut self, index: usize) { diff --git a/crates/ppvm-tableau/src/measure.rs b/crates/ppvm-tableau/src/measure.rs index 1b69319b2..8193564a0 100644 --- a/crates/ppvm-tableau/src/measure.rs +++ b/crates/ppvm-tableau/src/measure.rs @@ -56,14 +56,17 @@ where + ToPrimitive + std::fmt::Debug + std::ops::Mul - + PartialOrd, + + PartialOrd + + Send + + Sync, Complex: std::ops::Mul> + From + std::ops::MulAssign + std::ops::AddAssign + One - + ComplexFloat, - I: TableauIndex + Debug, + + ComplexFloat + + Copy, + I: TableauIndex + Debug + Send + Sync, { fn measure(&mut self, addr0: usize) -> Option { if self.is_lost[addr0] { diff --git a/crates/ppvm-tableau/src/noise.rs b/crates/ppvm-tableau/src/noise.rs index 0b1b39244..5d881f495 100644 --- a/crates/ppvm-tableau/src/noise.rs +++ b/crates/ppvm-tableau/src/noise.rs @@ -210,18 +210,27 @@ where } } -impl, I>> LossChannel +impl, I>> LossChannel for GeneralizedTableau where <::Storage as BitView>::Store: PrimInt, C: std::fmt::Debug, - T::Coeff: PartialOrd + One + Zero + Clone + num::Num + ToPrimitive + std::fmt::Debug, + T::Coeff: PartialOrd + + One + + Zero + + Clone + + num::Num + + ToPrimitive + + std::fmt::Debug + + Send + + Sync, Complex: std::ops::Mul> + From + std::ops::MulAssign + std::ops::AddAssign + One - + ComplexFloat, + + ComplexFloat + + Copy, I: Debug, { fn loss_channel(&mut self, addr0: usize, p: ::Coeff) { @@ -239,18 +248,27 @@ where } } -impl, I>> CorrelatedLossChannel - for GeneralizedTableau +impl, I>> + CorrelatedLossChannel for GeneralizedTableau where <::Storage as BitView>::Store: PrimInt, C: std::fmt::Debug, - T::Coeff: PartialOrd + One + Zero + Clone + num::Num + ToPrimitive + std::fmt::Debug, + T::Coeff: PartialOrd + + One + + Zero + + Clone + + num::Num + + ToPrimitive + + std::fmt::Debug + + Send + + Sync, Complex: std::ops::Mul> + From + std::ops::MulAssign + std::ops::AddAssign + One - + ComplexFloat, + + ComplexFloat + + Copy, I: Debug, { /// Apply a correlated loss channel to qubits at `addr0` and `addr1`. diff --git a/crates/ppvm-tableau/src/stim.rs b/crates/ppvm-tableau/src/stim.rs index 944940da5..d36db888f 100644 --- a/crates/ppvm-tableau/src/stim.rs +++ b/crates/ppvm-tableau/src/stim.rs @@ -64,14 +64,17 @@ where + ToPrimitive + std::fmt::Debug + std::ops::Mul - + PartialOrd, + + PartialOrd + + Send + + Sync, Complex: std::ops::Mul> + From + std::ops::MulAssign + std::ops::AddAssign + One - + ComplexFloat, - I: TableauIndex + Debug, + + ComplexFloat + + Copy, + I: TableauIndex + Debug + Send + Sync, { fn run_stim_string(&mut self, circuit: &str) -> Vec> { let mut results = Vec::>::new(); From 757399945ecafde14e4724a6d0b4b03189759936 Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 14:52:37 -0400 Subject: [PATCH 3/8] chore(autotune): record rayon iteration 1 results Co-Authored-By: Claude Opus 4.6 (1M context) --- .../2026-04-11-tableau-rayon-threading/log.md | 2 ++ .../2026-04-11-tableau-rayon-threading/metric.toml | 11 +++++++++++ 2 files changed, 13 insertions(+) diff --git a/docs/autotune/2026-04-11-tableau-rayon-threading/log.md b/docs/autotune/2026-04-11-tableau-rayon-threading/log.md index 95557ecf0..a60b8af09 100644 --- a/docs/autotune/2026-04-11-tableau-rayon-threading/log.md +++ b/docs/autotune/2026-04-11-tableau-rayon-threading/log.md @@ -29,3 +29,5 @@ 3. Measurement coefficient update (measure.rs) — processing b_keys entries. **Design constraint:** All parallelism behind `features = ["rayon"]`. When off, zero regression — identical codepath. +## 2026-04-12 +- Iteration 1: Added rayon feature with threshold-based parallel coefficient branching. No regression without feature. With feature enabled, rayon overhead dominates at current benchmark sizes (4096 coefficients). The Vec::collect + par_iter + HashMap fold/reduce pattern adds ~20µs per rayon call. Threshold set to 4096 items; below that sequential path is used. At 65536 coefficients (16 T gates), rayon is still ~2x slower than sequential due to HashMap merge overhead. The per-coefficient work is too cheap for rayon. Next steps: explore chunked accumulation or different data structures for the parallel path. diff --git a/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml b/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml index 51072dd39..22b88e30d 100644 --- a/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml +++ b/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml @@ -7,3 +7,14 @@ "t_gate_8_us" = 8.25 "t_gate_12_us" = 50.87 "msd_us" = 124.97 +[[metric]] +"commit" = "b574c45828cc59daa0c00a7084f911536a31013a" +"status" = "keep" +"description" = "rayon feature with threshold=4096 for coefficient branching" +"t_gate_1_ns" = 206.0 +"t_gate_4_us" = 1.08 +"t_gate_8_us" = 8.46 +"t_gate_12_us" = 51.3 +"msd_us" = 127.8 +"t_gate_12_rayon_us" = 52.3 +"msd_rayon_us" = 125.3 From aa17bc063c766bff1ab9c1283f15228be5397a6f Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 14:58:30 -0400 Subject: [PATCH 4/8] bench(tableau): add large-coefficient T-gate benchmarks Adds a `large-tgate` benchmark group that isolates the cost of a single T-gate applied to pre-built states with 2K-131K coefficients. Also adds a profile_rayon example for quick A/B comparisons. Co-Authored-By: Claude Opus 4.6 (1M context) --- crates/ppvm-tableau/benches/tableau.rs | 36 +++++++++ crates/ppvm-tableau/examples/profile_rayon.rs | 77 +++++++++++++++++++ 2 files changed, 113 insertions(+) create mode 100644 crates/ppvm-tableau/examples/profile_rayon.rs diff --git a/crates/ppvm-tableau/benches/tableau.rs b/crates/ppvm-tableau/benches/tableau.rs index d41889423..104e6243e 100644 --- a/crates/ppvm-tableau/benches/tableau.rs +++ b/crates/ppvm-tableau/benches/tableau.rs @@ -91,6 +91,42 @@ pub fn benchmark_suite_tableau(c: &mut Criterion, name: impl AsRef) { } group.finish(); + + // Large-coefficient T-gate benchmarks: measures the cost of a single T gate + // applied to a pre-built state with a known number of coefficients. + // This isolates the coefficient branching hot path for rayon testing. + let mut large_group = c.benchmark_group("large-tgate"); + large_group + .warm_up_time(Duration::from_secs(1)) + .measurement_time(Duration::from_secs(5)) + .sample_size(20); + + // Use 128 qubits (u128 index) so we have enough room for many distinct T targets + type LargeTab = GeneralizedTableau, u128>; + + for n_tgates in [12, 14, 16, 18] { + // Pre-build the state with (n_tgates - 1) T gates already applied, + // then benchmark applying the last one. + let mut setup: LargeTab = GeneralizedTableau::new(128, 1e-10); + for i in 0..n_tgates { + setup.h(i); + } + for i in 0..n_tgates - 1 { + setup.t(i); + } + let last_qubit = n_tgates - 1; + let n_coeffs = setup.coefficients.len(); + + large_group.bench_function(format!("single-t-on-{n_coeffs}-coeffs"), |b| { + b.iter_batched_ref( + || setup.fork(None), + |tab| tab.t(last_qubit), + criterion::BatchSize::LargeInput, + ); + }); + } + + large_group.finish(); } pub fn tableau_scaling_benchmarks(c: &mut Criterion) { diff --git a/crates/ppvm-tableau/examples/profile_rayon.rs b/crates/ppvm-tableau/examples/profile_rayon.rs new file mode 100644 index 000000000..274a50bcd --- /dev/null +++ b/crates/ppvm-tableau/examples/profile_rayon.rs @@ -0,0 +1,77 @@ +//! Profile the rayon vs sequential coefficient branching at various scales. +//! +//! This measures only the T-gate application phase (where coefficients branch), +//! isolating the exact code path that rayon parallelizes. +//! +//! Run without rayon: cargo run --release -p ppvm-tableau --example profile_rayon +//! Run with rayon: cargo run --release -p ppvm-tableau --example profile_rayon --features rayon + +use std::time::Instant; + +use ppvm_runtime::config::fx64hash::Byte8F64; +use ppvm_tableau::prelude::*; + +type Tab = GeneralizedTableau, u128>; + +fn main() { + // Use 128 qubits so we have enough room for many T gates on distinct qubits + let n_qubits = 128; + + println!("=== T-gate branching: sequential vs rayon ({n_qubits} qubits) ==="); + println!( + "{:>4} {:>8} {:>12} {:>12}", + "T#", "coeffs", "total (µs)", "last T (µs)" + ); + + for n_tgates in [1, 4, 8, 12, 14, 16, 18, 20] { + // Prepare: H on each qubit that will get a T gate (ensures branching) + let mut tab: Tab = GeneralizedTableau::new(n_qubits, 1e-10); + for i in 0..n_tgates { + tab.h(i); + } + + // Warm up: one run to prime caches + { + let mut t = tab.fork(Some(42)); + for i in 0..n_tgates { + t.t(i); + } + } + + // Measure: time all T gates together + let n_runs = if n_tgates <= 16 { 5 } else { 3 }; + let mut total_all_ns = 0u128; + let mut total_last_ns = 0u128; + let mut coeff_count = 0; + + for _ in 0..n_runs { + let mut t = tab.fork(Some(42)); + + // Apply T gates 0..n-1 (setup, not timed individually) + let t0 = Instant::now(); + for i in 0..n_tgates.saturating_sub(1) { + t.t(i); + } + let setup_ns = t0.elapsed().as_nanos(); + + // Time only the last T gate (largest coefficient set) + let t1 = Instant::now(); + if n_tgates > 0 { + t.t(n_tgates - 1); + } + let last_ns = t1.elapsed().as_nanos(); + + total_all_ns += setup_ns + last_ns; + total_last_ns += last_ns; + coeff_count = t.coefficients.len(); + } + + println!( + "{:>4} {:>8} {:>12.1} {:>12.1}", + n_tgates, + coeff_count, + total_all_ns as f64 / n_runs as f64 / 1000.0, + total_last_ns as f64 / n_runs as f64 / 1000.0, + ); + } +} From d5db5a83efc4d7838bb0aca607e777e50139dd75 Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 15:02:53 -0400 Subject: [PATCH 5/8] =?UTF-8?q?perf(tableau):=20fix=20rayon=20approach=20?= =?UTF-8?q?=E2=80=94=20parallel=20map=20+=20sequential=20accumulate?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the rayon fold/reduce-with-HashMap pattern (which was 2-3x slower due to thread-local HashMap merge overhead) with parallel map into a Vec of computed pairs followed by sequential accumulation into a pre-sized HashMap. Results at various coefficient counts (single T-gate application): - 32K coefficients: 35% faster than sequential - 131K coefficients: 39% faster - 1M coefficients: 47% faster Also raise threshold from 4096 to 16384 to avoid regressions at smaller coefficient counts where rayon overhead dominates. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../examples/profile_breakdown.rs | 129 ++++++++++++++++++ crates/ppvm-tableau/src/data.rs | 117 ++++++++-------- 2 files changed, 189 insertions(+), 57 deletions(-) create mode 100644 crates/ppvm-tableau/examples/profile_breakdown.rs diff --git a/crates/ppvm-tableau/examples/profile_breakdown.rs b/crates/ppvm-tableau/examples/profile_breakdown.rs new file mode 100644 index 000000000..a12cf1f9c --- /dev/null +++ b/crates/ppvm-tableau/examples/profile_breakdown.rs @@ -0,0 +1,129 @@ +//! Profile the time breakdown inside branch_with_coefficients: +//! computation (phase + multiply) vs HashMap accumulation. + +use std::time::Instant; + +use fxhash::FxHashMap as HashMap; +use num::Zero; +use num::complex::{Complex, Complex64}; + +fn main() { + // Simulate the coefficient branching loop at various sizes + println!( + "{:>10} {:>10} {:>10} {:>10} {:>6}", + "N", "compute", "accum", "total", "%comp" + ); + + for n in [2048, 8192, 32768, 131072] { + // Create fake coefficient data + let items: Vec<(Complex, u128)> = (0..n) + .map(|i| (Complex::new(1.0 / n as f64, 0.0), i as u128)) + .collect(); + + let stab_anticomm_bits: u128 = 0b1010_1010; + let destab_anticomm_bits: u128 = 0b0101_0101; + let odd_phase_mask: u128 = 0b1100_1100; + let phase_decomp: u8 = 1; + let coefficient_factor = Complex::new(0.85, 0.35); + let branch_factor = Complex::new(0.15, -0.35); + + let phase_table: [Complex64; 4] = [ + Complex64::new(1.0, 0.0), + Complex64::new(0.0, 1.0), + Complex64::new(-1.0, 0.0), + Complex64::new(0.0, -1.0), + ]; + + let n_runs = 10; + + // Phase 1: compute only (write to Vec, no HashMap) + let mut compute_ns = 0u128; + for _ in 0..n_runs { + let mut results = Vec::with_capacity(items.len()); + let t0 = Instant::now(); + for &(coeff, idx) in &items { + let branch_index = idx ^ stab_anticomm_bits; + let symplectic = (destab_anticomm_bits & idx).count_ones(); + let mut phase = (2 * symplectic as u8) % 4; + let active = idx & stab_anticomm_bits; + let parity = (active & odd_phase_mask).count_ones() % 2; + phase = (phase + 2 * parity as u8) % 4; + let branch_phase = (phase + phase_decomp) % 4; + let phase_factor: Complex = phase_table[branch_phase as usize].into(); + let branch_coeff = phase_factor * coeff * branch_factor; + let nonbranch_coeff = coeff * coefficient_factor; + results.push((branch_index, branch_coeff, idx, nonbranch_coeff)); + } + compute_ns += t0.elapsed().as_nanos(); + std::hint::black_box(&results); + } + + // Phase 2: accumulate only (from pre-computed data) + let mut precomputed: Vec<(u128, Complex, u128, Complex)> = + Vec::with_capacity(items.len()); + for &(coeff, idx) in &items { + let branch_index = idx ^ stab_anticomm_bits; + let symplectic = (destab_anticomm_bits & idx).count_ones(); + let mut phase = (2 * symplectic as u8) % 4; + let active = idx & stab_anticomm_bits; + let parity = (active & odd_phase_mask).count_ones() % 2; + phase = (phase + 2 * parity as u8) % 4; + let branch_phase = (phase + phase_decomp) % 4; + let phase_factor: Complex = phase_table[branch_phase as usize].into(); + precomputed.push(( + branch_index, + phase_factor * coeff * branch_factor, + idx, + coeff * coefficient_factor, + )); + } + + let mut accum_ns = 0u128; + for _ in 0..n_runs { + let mut map: HashMap> = HashMap::default(); + map.reserve(2 * items.len()); + let t0 = Instant::now(); + for &(branch_idx, branch_coeff, idx, nonbranch_coeff) in &precomputed { + *map.entry(branch_idx).or_insert(Complex::zero()) += branch_coeff; + *map.entry(idx).or_insert(Complex::zero()) += nonbranch_coeff; + } + accum_ns += t0.elapsed().as_nanos(); + std::hint::black_box(&map); + } + + // Phase 3: combined (original pattern) + let mut total_ns = 0u128; + for _ in 0..n_runs { + let mut map: HashMap> = HashMap::default(); + let t0 = Instant::now(); + for &(coeff, idx) in &items { + let branch_index = idx ^ stab_anticomm_bits; + let symplectic = (destab_anticomm_bits & idx).count_ones(); + let mut phase = (2 * symplectic as u8) % 4; + let active = idx & stab_anticomm_bits; + let parity = (active & odd_phase_mask).count_ones() % 2; + phase = (phase + 2 * parity as u8) % 4; + let branch_phase = (phase + phase_decomp) % 4; + let phase_factor: Complex = phase_table[branch_phase as usize].into(); + let branch_coeff = phase_factor * coeff * branch_factor; + let nonbranch_coeff = coeff * coefficient_factor; + *map.entry(branch_index).or_insert(Complex::zero()) += branch_coeff; + *map.entry(idx).or_insert(Complex::zero()) += nonbranch_coeff; + } + total_ns += t0.elapsed().as_nanos(); + std::hint::black_box(&map); + } + + let compute_us = compute_ns as f64 / n_runs as f64 / 1000.0; + let accum_us = accum_ns as f64 / n_runs as f64 / 1000.0; + let total_us = total_ns as f64 / n_runs as f64 / 1000.0; + println!( + "{:>10} {:>9.0}µ {:>9.0}µ {:>9.0}µ {:>5.1}%", + n, + compute_us, + accum_us, + total_us, + compute_us / total_us * 100.0, + ); + } +} diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index fb6ddee20..e94ac973e 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -295,10 +295,11 @@ fn compute_phase_with_mask_static( } /// Minimum number of coefficients before engaging rayon parallelism. -/// Below this threshold, the sequential path is always used even with the `rayon` feature, -/// because rayon's thread-pool overhead (~20µs) dominates the per-coefficient work. +/// Below this threshold, the sequential path is always used even with the `rayon` feature. +/// Benchmarked: at 8K coefficients rayon has ~24% overhead; at 32K it's 35% faster. +/// Set to 16384 to avoid regressions while capturing the large-coefficient wins. #[cfg(feature = "rayon")] -const RAYON_COEFF_THRESHOLD: usize = 4096; +const RAYON_COEFF_THRESHOLD: usize = 16384; /// Sequential accumulation of branch coefficients. fn branch_coefficients_seq( @@ -358,37 +359,42 @@ where Complex: std::ops::Mul> + std::ops::AddAssign + From + Copy, { - #[cfg(feature = "rayon")] if items.len() >= RAYON_COEFF_THRESHOLD { use rayon::prelude::*; - return items + + // Parallel phase: compute all (branch_idx, branch_coeff, idx, nonbranch_coeff) tuples. + // This is pure math with no shared mutable state. + let pairs: Vec<(I, Complex, I, Complex)> = items .par_iter() - .fold( - HashMap::>::default, - |mut map, &(coeff, idx)| { - let branch_index = idx ^ stab_anticomm_bits; - let branch_phase_contribution = compute_phase_with_mask_static( - destab_anticomm_bits, - idx, - stab_anticomm_bits, - odd_phase_mask, - ); - let branch_phase = (branch_phase_contribution + phase_decomp) % 4; - let phase_factor: Complex = - COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); - let branch_coefficient = phase_factor * coeff * branch_factor; - let nonbranch_coefficient = coeff * coefficient_factor; - *map.entry(branch_index).or_insert(Complex::zero()) += branch_coefficient; - *map.entry(idx).or_insert(Complex::zero()) += nonbranch_coefficient; - map - }, - ) - .reduce(HashMap::default, |mut a, b| { - for (k, v) in b { - *a.entry(k).or_insert(Complex::zero()) += v; - } - a - }); + .map(|&(coeff, idx)| { + let branch_index = idx ^ stab_anticomm_bits; + let branch_phase_contribution = compute_phase_with_mask_static( + destab_anticomm_bits, + idx, + stab_anticomm_bits, + odd_phase_mask, + ); + let branch_phase = (branch_phase_contribution + phase_decomp) % 4; + let phase_factor: Complex = + COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); + ( + branch_index, + phase_factor * coeff * branch_factor, + idx, + coeff * coefficient_factor, + ) + }) + .collect(); + + // Sequential phase: accumulate into a pre-sized HashMap. + // HashMap inserts dominate the cost but benefit from cache locality. + let mut map: HashMap> = + HashMap::with_capacity_and_hasher(2 * pairs.len(), Default::default()); + for (branch_idx, branch_coeff, idx, nonbranch_coeff) in pairs { + *map.entry(branch_idx).or_insert(Complex::zero()) += branch_coeff; + *map.entry(idx).or_insert(Complex::zero()) += nonbranch_coeff; + } + return map; } branch_coefficients_seq( @@ -454,35 +460,32 @@ where Complex: std::ops::Mul> + std::ops::AddAssign + From + Copy, { - #[cfg(feature = "rayon")] if items.len() >= RAYON_COEFF_THRESHOLD { use rayon::prelude::*; - return items + + let pairs: Vec<(I, Complex)> = items .par_iter() - .fold( - HashMap::>::default, - |mut map, &(coeff, idx)| { - let branch_index = idx ^ stab_anticomm_bits; - let branch_phase_contribution = compute_phase_with_mask_static( - destab_anticomm_bits, - idx, - stab_anticomm_bits, - odd_phase_mask, - ); - let branch_phase = (branch_phase_contribution + phase_decomp) % 4; - let phase_factor: Complex = - COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); - let branch_coefficient = phase_factor * coeff; - *map.entry(branch_index).or_insert(Complex::zero()) += branch_coefficient; - map - }, - ) - .reduce(HashMap::default, |mut a, b| { - for (k, v) in b { - *a.entry(k).or_insert(Complex::zero()) += v; - } - a - }); + .map(|&(coeff, idx)| { + let branch_index = idx ^ stab_anticomm_bits; + let branch_phase_contribution = compute_phase_with_mask_static( + destab_anticomm_bits, + idx, + stab_anticomm_bits, + odd_phase_mask, + ); + let branch_phase = (branch_phase_contribution + phase_decomp) % 4; + let phase_factor: Complex = + COMPLEX_PHASE_CONVERSION[branch_phase as usize].into(); + (branch_index, phase_factor * coeff) + }) + .collect(); + + let mut map: HashMap> = + HashMap::with_capacity_and_hasher(pairs.len(), Default::default()); + for (branch_idx, branch_coeff) in pairs { + *map.entry(branch_idx).or_insert(Complex::zero()) += branch_coeff; + } + return map; } apply_coefficients_seq( From e1545ca4ee47ae201f5cd651f7ee5b0fe6548678 Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 15:03:23 -0400 Subject: [PATCH 6/8] =?UTF-8?q?chore(autotune):=20record=20iteration=202?= =?UTF-8?q?=20results=20=E2=80=94=20rayon=20now=2047%=20faster=20at=201M?= =?UTF-8?q?=20coeffs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 4.6 (1M context) --- .../2026-04-11-tableau-rayon-threading/log.md | 1 + .../2026-04-11-tableau-rayon-threading/metric.toml | 13 +++++++++++++ 2 files changed, 14 insertions(+) diff --git a/docs/autotune/2026-04-11-tableau-rayon-threading/log.md b/docs/autotune/2026-04-11-tableau-rayon-threading/log.md index a60b8af09..b824e1c1a 100644 --- a/docs/autotune/2026-04-11-tableau-rayon-threading/log.md +++ b/docs/autotune/2026-04-11-tableau-rayon-threading/log.md @@ -31,3 +31,4 @@ **Design constraint:** All parallelism behind `features = ["rayon"]`. When off, zero regression — identical codepath. ## 2026-04-12 - Iteration 1: Added rayon feature with threshold-based parallel coefficient branching. No regression without feature. With feature enabled, rayon overhead dominates at current benchmark sizes (4096 coefficients). The Vec::collect + par_iter + HashMap fold/reduce pattern adds ~20µs per rayon call. Threshold set to 4096 items; below that sequential path is used. At 65536 coefficients (16 T gates), rayon is still ~2x slower than sequential due to HashMap merge overhead. The per-coefficient work is too cheap for rayon. Next steps: explore chunked accumulation or different data structures for the parallel path. +- Iteration 2: Replaced rayon fold/reduce-with-HashMap with parallel map → Vec + sequential accumulate. The HashMap merge in fold/reduce was the bottleneck (not the computation). New approach: rayon par_map computes all (key, value) pairs into a Vec, then sequential loop inserts into pre-sized HashMap. Results: 35% faster at 32K coeffs, 39% at 131K, 47% at 1M. Threshold raised to 16384. Profile breakdown showed accumulation is 70-85% of total time, computation only 15-25%. diff --git a/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml b/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml index 22b88e30d..9b69c27ce 100644 --- a/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml +++ b/docs/autotune/2026-04-11-tableau-rayon-threading/metric.toml @@ -18,3 +18,16 @@ "msd_us" = 127.8 "t_gate_12_rayon_us" = 52.3 "msd_rayon_us" = 125.3 +[[metric]] +"commit" = "d5db5a83efc4d7838bb0aca607e777e50139dd75" +"status" = "keep" +"description" = "parallel map + sequential accumulate, threshold=16384" +"t_gate_12_us" = 53.3 +"single_t_on_2048_us" = 62.0 +"single_t_on_8192_us" = 275.0 +"single_t_on_32768_us" = 748.0 +"single_t_on_131072_us" = 2517.0 +"single_t_on_32768_seq_us" = 1150.0 +"single_t_on_131072_seq_us" = 4140.0 +"total_20t_seq_us" = 45819.0 +"total_20t_rayon_us" = 24126.0 From 790d7ea122236898e67b546fb321d7aeb66f2289 Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 18:33:19 -0400 Subject: [PATCH 7/8] =?UTF-8?q?fix(tableau):=20address=20PR=20review=20?= =?UTF-8?q?=E2=80=94=20avoid=20Vec=20alloc=20below=20threshold,=20fix=20do?= =?UTF-8?q?cs?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Skip Vec::collect when coefficient count is below RAYON_COEFF_THRESHOLD, calling the sequential path directly (fixes unnecessary allocation) - Update doc comments to reflect parallel map/collect strategy (not fold/reduce) - Use deterministic seed in large-tgate benchmark (fork(Some(0)) not fork(None)) Co-Authored-By: Claude Opus 4.6 (1M context) --- crates/ppvm-tableau/benches/tableau.rs | 2 +- crates/ppvm-tableau/src/data.rs | 39 +++++++++++++++++++++----- 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/crates/ppvm-tableau/benches/tableau.rs b/crates/ppvm-tableau/benches/tableau.rs index 104e6243e..f32174796 100644 --- a/crates/ppvm-tableau/benches/tableau.rs +++ b/crates/ppvm-tableau/benches/tableau.rs @@ -119,7 +119,7 @@ pub fn benchmark_suite_tableau(c: &mut Criterion, name: impl AsRef) { large_group.bench_function(format!("single-t-on-{n_coeffs}-coeffs"), |b| { b.iter_batched_ref( - || setup.fork(None), + || setup.fork(Some(0)), |tab| tab.t(last_qubit), criterion::BatchSize::LargeInput, ); diff --git a/crates/ppvm-tableau/src/data.rs b/crates/ppvm-tableau/src/data.rs index e94ac973e..4d277622a 100644 --- a/crates/ppvm-tableau/src/data.rs +++ b/crates/ppvm-tableau/src/data.rs @@ -341,8 +341,9 @@ where map } -/// Accumulate branch coefficients. Uses rayon parallel fold/reduce when the feature -/// is enabled AND the coefficient count exceeds `RAYON_COEFF_THRESHOLD`. +/// Accumulate branch coefficients. When the coefficient count exceeds +/// `RAYON_COEFF_THRESHOLD`, uses parallel map/collect into a Vec followed +/// by sequential accumulation. Below the threshold, falls back to sequential. #[cfg(feature = "rayon")] fn branch_coefficients_parallel( items: &[(Complex, I)], @@ -444,8 +445,9 @@ where map } -/// Accumulate coefficients for pauli application. Uses rayon when the feature is -/// enabled AND the coefficient count exceeds `RAYON_COEFF_THRESHOLD`. +/// Accumulate coefficients for pauli application. When the coefficient count +/// exceeds `RAYON_COEFF_THRESHOLD`, uses parallel map/collect followed by +/// sequential accumulation. Below the threshold, falls back to sequential. #[cfg(feature = "rayon")] fn apply_coefficients_parallel( items: &[(Complex, I)], @@ -807,10 +809,13 @@ where let odd_phase_mask = self.odd_phase_destabilizer_mask(); let old_coefficients = std::mem::replace(&mut self.coefficients, C::new()); + #[cfg(feature = "rayon")] + let n_coefficients = old_coefficients.len(); - // Only collect to Vec when rayon might be used; otherwise iterate directly + // When rayon is enabled and above the threshold, collect to Vec for parallel map; + // otherwise iterate directly with the sequential path. #[cfg(feature = "rayon")] - let new_coefficients = { + let new_coefficients = if n_coefficients >= RAYON_COEFF_THRESHOLD { let items: Vec<_> = old_coefficients.into_iter().collect(); branch_coefficients_parallel( &items, @@ -821,6 +826,16 @@ where coefficient_factor, branch_factor, ) + } else { + branch_coefficients_seq( + old_coefficients, + stab_anticomm_bits, + destab_anticomm_bits, + odd_phase_mask, + phase_decomp, + coefficient_factor, + branch_factor, + ) }; #[cfg(not(feature = "rayon"))] @@ -862,10 +877,12 @@ where self.compute_decomposition(addr0, pauli); let odd_phase_mask = self.odd_phase_destabilizer_mask(); + #[cfg(feature = "rayon")] + let n_coefficients = coefficients.len(); let old_coefficients = std::mem::replace(coefficients, C::new()); #[cfg(feature = "rayon")] - let new_coefficients = { + let new_coefficients = if n_coefficients >= RAYON_COEFF_THRESHOLD { let items: Vec<_> = old_coefficients.into_iter().collect(); apply_coefficients_parallel( &items, @@ -874,6 +891,14 @@ where odd_phase_mask, phase_decomp, ) + } else { + apply_coefficients_seq( + old_coefficients, + stab_anticomm_bits, + destab_anticomm_bits, + odd_phase_mask, + phase_decomp, + ) }; #[cfg(not(feature = "rayon"))] From 23471024491acc16876247bcac487c4eaaef0e3d Mon Sep 17 00:00:00 2001 From: Roger-luo Date: Sun, 12 Apr 2026 19:05:25 -0400 Subject: [PATCH 8/8] bench(tableau): add fused+tgate combined circuit benchmark Adds a `fused-tgate-circuit` benchmark group that combines batch Clifford operations (fusion) with variable T-gates and measurement. Also adds profile_combined example showing time breakdown. Key finding: measurement dominates at 85-90% of total circuit time when coefficient counts are large (>4K). Fusion reduces Clifford overhead to <1%. T-gate branching (rayon target) is ~10%. Co-Authored-By: Claude Opus 4.6 (1M context) --- crates/ppvm-tableau/benches/tableau.rs | 50 ++++++++++++ .../ppvm-tableau/examples/profile_combined.rs | 78 +++++++++++++++++++ 2 files changed, 128 insertions(+) create mode 100644 crates/ppvm-tableau/examples/profile_combined.rs diff --git a/crates/ppvm-tableau/benches/tableau.rs b/crates/ppvm-tableau/benches/tableau.rs index f32174796..c0e30f720 100644 --- a/crates/ppvm-tableau/benches/tableau.rs +++ b/crates/ppvm-tableau/benches/tableau.rs @@ -127,6 +127,56 @@ pub fn benchmark_suite_tableau(c: &mut Criterion, name: impl AsRef) { } large_group.finish(); + + // Combined benchmark: fused Clifford batches + variable T gates + measurement. + // Shows how fusion and rayon compose in a realistic circuit. + let mut combined_group = c.benchmark_group("fused-tgate-circuit"); + combined_group + .warm_up_time(Duration::from_secs(1)) + .measurement_time(Duration::from_secs(5)) + .sample_size(20); + + for n_tgates in [8, 12, 16] { + // Circuit: batched Cliffords → H (ensure branching) → T gates → batched Cliffords → measure. + // H is applied AFTER Cliffords to ensure T gates see non-stabilizer states and branch. + let n_qubits: usize = 85; + + // Pre-build the Clifford portion so the setup cost is excluded + let mut setup = GeneralizedTableau::, u128>::new(n_qubits, 1e-10); + let block1: Vec = (0..17).collect(); + let block2: Vec = (17..34).collect(); + setup.sqrt_y_batch(&block1); + setup.sqrt_x_batch(&block2); + setup.cz_block_pairs(0, 17, 17); + // H on qubits that will get T gates — applied after Cliffords to guarantee branching + for i in 0..n_tgates { + setup.h(i); + } + + combined_group.bench_function(format!("fused-{n_tgates}t-{n_qubits}q"), |b| { + b.iter_batched_ref( + || setup.fork(Some(0)), + |tab| { + // T gates (coefficient branching — rayon target) + for i in 0..n_tgates { + tab.t(i); + } + + // More batched Cliffords after T gates + tab.sqrt_x_adj_batch(&block1); + tab.sqrt_y_adj_batch(&block2); + + // Measure all + for i in 0..n_qubits { + tab.measure(i); + } + }, + criterion::BatchSize::LargeInput, + ); + }); + } + + combined_group.finish(); } pub fn tableau_scaling_benchmarks(c: &mut Criterion) { diff --git a/crates/ppvm-tableau/examples/profile_combined.rs b/crates/ppvm-tableau/examples/profile_combined.rs new file mode 100644 index 000000000..147da6f38 --- /dev/null +++ b/crates/ppvm-tableau/examples/profile_combined.rs @@ -0,0 +1,78 @@ +//! Profile the time breakdown of a fused circuit with variable T gates. + +use std::time::Instant; + +use ppvm_runtime::config::fx64hash::Byte8F64; +use ppvm_tableau::prelude::*; + +type Tab = GeneralizedTableau, u128>; + +fn main() { + let n_qubits = 85; + + println!( + "{:>4} {:>10} {:>10} {:>10} {:>10} {:>6} {:>6} {:>6}", + "T#", "clifford", "t-gates", "measure", "total", "%clif", "%t", "%meas" + ); + + for n_tgates in [4, 8, 12, 14, 16] { + let block1: Vec = (0..17).collect(); + let block2: Vec = (17..34).collect(); + + let n_runs = if n_tgates <= 12 { 5 } else { 2 }; + let mut clif_ns = 0u128; + let mut tgate_ns = 0u128; + let mut meas_ns = 0u128; + + for _ in 0..n_runs { + let mut tab: Tab = GeneralizedTableau::new(n_qubits, 1e-10); + + // Clifford layer 1 (fused) + let t0 = Instant::now(); + tab.sqrt_y_batch(&block1); + tab.sqrt_x_batch(&block2); + tab.cz_block_pairs(0, 17, 17); + for i in 0..n_tgates { + tab.h(i); + } + clif_ns += t0.elapsed().as_nanos(); + + // T gates + let t1 = Instant::now(); + for i in 0..n_tgates { + tab.t(i); + } + tgate_ns += t1.elapsed().as_nanos(); + + // Clifford layer 2 (fused) + let t2 = Instant::now(); + tab.sqrt_x_adj_batch(&block1); + tab.sqrt_y_adj_batch(&block2); + clif_ns += t2.elapsed().as_nanos(); + + // Measure + let t3 = Instant::now(); + for i in 0..n_qubits { + tab.measure(i); + } + meas_ns += t3.elapsed().as_nanos(); + } + + let d = n_runs as f64; + let c = clif_ns as f64 / d / 1000.0; + let t = tgate_ns as f64 / d / 1000.0; + let m = meas_ns as f64 / d / 1000.0; + let total = c + t + m; + println!( + "{:>4} {:>9.0}µ {:>9.0}µ {:>9.0}µ {:>9.0}µ {:>5.1}% {:>5.1}% {:>5.1}%", + n_tgates, + c, + t, + m, + total, + c / total * 100.0, + t / total * 100.0, + m / total * 100.0, + ); + } +}