refactor(ppvm-lindblad): split lib.rs along spec/basis/step lines

crates/ppvm-lindblad/Cargo.toml

@@ -10,7 +10,6 @@ ndarray = "0.17"
 num = "0.4.3"
 ppvm-traits     = { version = "0.1.0", path = "../ppvm-traits" }
 ppvm-pauli-word = { version = "0.1.0", path = "../ppvm-pauli-word" }
-ppvm-pauli-sum  = { version = "0.1.0", path = "../ppvm-pauli-sum" }
 rayon = "1.11"
 # Matrix-exponential action (Al-Mohy & Higham). QuSpin-rust is MIT-licensed;
 # the pinned rev is the commit that added the LICENSE file.
@@ -21,4 +20,4 @@ quspin-expm = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2
 quspin-types = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2e8063208f9ba1c6677150c993bb554" }
 
 [dev-dependencies]
-approx = "0.5.1"
+ppvm-pauli-sum = { version = "0.1.0", path = "../ppvm-pauli-sum" }

crates/ppvm-lindblad/src/algebra.rs

@@ -0,0 +1,78 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Hot-path Pauli product / commutator on raw word chunks.
+//!
+//! Phase encoding: Pauli product `P·Q = ω · R` where `ω = i^phase` and
+//! `phase ∈ {0,1,2,3}` ↔ `ω ∈ {1, i, -1, -i}`. The per-byte XOR/AND
+//! formulas are the same ones used by
+//! [`ppvm_pauli_word::phase::PhasedPauliWord`]'s `MulAssign`. This module
+//! keeps a copy that returns the unpacked `(word, phase)` pair without
+//! constructing a phased wrapper.
+
+use crate::word::{W_CHUNKS, Word};
+use num::Complex;
+use ppvm_traits::PauliWordTrait;
+
+#[inline(always)]
+pub(crate) fn phase_factor(phase: u8) -> Complex<f64> {
+    match phase & 3 {
+        0 => Complex::new(1.0, 0.0),
+        1 => Complex::new(0.0, 1.0),
+        2 => Complex::new(-1.0, 0.0),
+        _ => Complex::new(0.0, -1.0),
+    }
+}
+
+/// `true` if Pauli words `a` and `b` anti-commute.
+///
+/// Two Pauli strings anti-commute iff
+/// `popcount(a.x & b.z) + popcount(a.z & b.x)` is odd.
+#[inline(always)]
+pub(crate) fn anti_commutes(a: &Word, b: &Word) -> bool {
+    let mut bits: u32 = 0;
+    for i in 0..W_CHUNKS {
+        bits += (a.xbits.data[i] & b.zbits.data[i]).count_ones();
+        bits += (a.zbits.data[i] & b.xbits.data[i]).count_ones();
+    }
+    bits & 1 == 1
+}
+
+/// Commutator product `h · p`: returns `(out, eps)` where `out = h ⊕ p` and
+///
+/// - `eps =  0` if `h` and `p` commute (caller should skip — `[h,p] = 0`),
+/// - `eps = -2.0` if `h·p` has phase `+i` (so `i·[h,p] = -2·out`),
+/// - `eps = +2.0` if `h·p` has phase `-i` (so `i·[h,p] = +2·out`).
+#[inline(always)]
+pub(crate) fn comm_product(h: &Word, p: &Word) -> (Word, f64) {
+    let (out, phase) = pauli_mul(h, p);
+    let eps = match phase {
+        1 => -2.0,
+        3 => 2.0,
+        _ => 0.0,
+    };
+    (out, eps)
+}
+
+/// Full Pauli product `p · q`: returns `(out, phase)` where the product
+/// is `ω · out` with `ω = i^phase`.
+#[inline(always)]
+pub(crate) fn pauli_mul(p: &Word, q: &Word) -> (Word, u8) {
+    let mut out = Word::new(p.n_qubits());
+    let mut sign_count: u32 = 0;
+    let mut imag_count: u32 = 0;
+    for i in 0..W_CHUNKS {
+        let a = p.xbits.data[i];
+        let b = p.zbits.data[i];
+        let c = q.xbits.data[i];
+        let d = q.zbits.data[i];
+        let sign = (a & b & c & !d) | (a & !b & !c & d) | (!a & b & c & d);
+        let imag = (a & !b & d) | (a & !c & d) | (!a & b & c) | (b & c & !d);
+        sign_count += sign.count_ones();
+        imag_count += imag.count_ones();
+        out.xbits.data[i] = a ^ c;
+        out.zbits.data[i] = b ^ d;
+    }
+    out.rehash();
+    (out, ((2 * sign_count + imag_count) & 3) as u8)
+}

crates/ppvm-lindblad/src/basis.rs

@@ -0,0 +1,256 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Basis-level `L*` operators: in-basis generator and off-basis leakage.
+
+use crate::Error;
+use crate::spec::LindbladSpec;
+use crate::word::{Word, word_hash};
+use fxhash::{FxBuildHasher, FxHashMap};
+use num::Complex;
+use rayon::prelude::*;
+
+/// Build a `word → row` map for a basis assumed to contain unique Pauli
+/// words; debug-asserts the uniqueness invariant.
+pub fn build_basis_index(basis: &[Word]) -> FxHashMap<Word, u32> {
+    let mut index: FxHashMap<Word, u32> = FxHashMap::default();
+    for (i, w) in basis.iter().enumerate() {
+        let prev = index.insert(*w, i as u32);
+        debug_assert!(
+            prev.is_none(),
+            "basis contains duplicate Pauli word at positions {} and {}",
+            prev.unwrap(),
+            i,
+        );
+    }
+    index
+}
+
+impl LindbladSpec {
+    /// Off-basis component of `L*( Σ_j coeffs[j] · basis[j] )`. Output
+    /// strings that lie in `basis` or in `protected` are dropped.
+    pub fn leakage(
+        &self,
+        basis: &[Word],
+        coeffs: &[f64],
+        protected: &[Word],
+    ) -> Result<Vec<(Word, f64)>, Error> {
+        self.leakage_with_prune(basis, coeffs, protected, usize::MAX, 0.0)
+    }
+
+    /// Like [`Self::leakage`], but caps the live off-basis leakage map to
+    /// the *available room* `room = max_basis − basis.len()` — only the
+    /// strings we could actually add to the basis are worth keeping. The
+    /// cap is applied during accumulation (after each chunk), keeping the
+    /// `room` largest-magnitude entries.
+    ///
+    /// Basis indices are processed in descending-`|c|` order so the
+    /// running cap keeps the entries that are most likely to be the true
+    /// largest contributors. When `max_basis` is large enough that
+    /// `room ≥ all candidates`, nothing is dropped — the near-exact case.
+    pub fn leakage_with_prune(
+        &self,
+        basis: &[Word],
+        coeffs: &[f64],
+        protected: &[Word],
+        max_basis: usize,
+        tau_add: f64,
+    ) -> Result<Vec<(Word, f64)>, Error> {
+        if basis.len() != coeffs.len() {
+            return Err(Error::LengthMismatch {
+                what: "basis and coeffs",
+                a: basis.len(),
+                b: coeffs.len(),
+            });
+        }
+        // Hash-only membership tables: storing 8-byte `u64` keys instead
+        // of 48-byte Words shrinks the in-basis structure ~6×, keeping it
+        // in L3 (and often L2) at basis sizes where the full-Word version
+        // would spill to DRAM.
+        let in_basis: FxHashMap<u64, ()> = basis.iter().map(|w| (word_hash(w), ())).collect();
+        let protected_set: FxHashMap<u64, ()> =
+            protected.iter().map(|w| (word_hash(w), ())).collect();
+
+        // Descending sort by |c|: process largest-magnitude contributors
+        // first so the running room-cap keeps the right entries.
+        let mut order: Vec<usize> = (0..basis.len()).collect();
+        order.sort_by(|&a, &b| {
+            coeffs[b]
+                .abs()
+                .partial_cmp(&coeffs[a].abs())
+                .unwrap_or(std::cmp::Ordering::Equal)
+        });
+
+        const CHUNK_SIZE: usize = 4096;
+        let room = max_basis.saturating_sub(basis.len());
+        let n_qubits = self.n_qubits();
+        let mut merged: FxHashMap<Word, f64> = FxHashMap::default();
+        for chunk_indices in order.chunks(CHUNK_SIZE) {
+            let local: Vec<Vec<(Word, f64)>> = chunk_indices
+                .par_iter()
+                .map_init(
+                    || {
+                        (
+                            Vec::<u32>::with_capacity(n_qubits),
+                            Vec::<u32>::with_capacity(128),
+                            FxHashMap::<Word, Complex<f64>>::with_capacity_and_hasher(
+                                128,
+                                FxBuildHasher::default(),
+                            ),
+                        )
+                    },
+                    |(s1, s2, lm), &i| {
+                        let p = &basis[i];
+                        let c = coeffs[i];
+                        let terms = self.compute_action_terms(p, s1, s2, lm);
+                        let mut out = Vec::with_capacity(terms.len());
+                        for (w, v) in terms.iter() {
+                            let h = word_hash(w);
+                            if !in_basis.contains_key(&h) && !protected_set.contains_key(&h) {
+                                out.push((*w, c * *v));
+                            }
+                        }
+                        out
+                    },
+                )
+                .collect();
+            for v in local {
+                for (k, val) in v {
+                    *merged.entry(k).or_insert(0.0) += val;
+                }
+            }
+
+            // Room-cap: keep only the `room` largest-magnitude entries.
+            if merged.len() > room {
+                if room == 0 {
+                    merged.clear();
+                } else {
+                    let mut mags: Vec<f64> = merged.values().map(|v| v.abs()).collect();
+                    let k = room.min(mags.len() - 1);
+                    mags.select_nth_unstable_by(k, |a, b| {
+                        b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal)
+                    });
+                    let cutoff = mags[k];
+                    merged.retain(|_, &mut v| v.abs() >= cutoff);
+                }
+            }
+        }
+        // Rate-based admission: keep only candidates whose leakage rate
+        // exceeds `tau_add`. `tau_add = 0` admits everything except exact
+        // zeros.
+        Ok(merged
+            .into_iter()
+            .filter(|(_, c)| c.abs() > tau_add)
+            .collect())
+    }
+
+    /// Sparse generator matrix in COO form: returns `(row, col, val)`
+    /// triplets. Row = output Pauli's position in `basis`; col = input
+    /// Pauli's position. Output Paulis not in `basis` are silently dropped.
+    ///
+    /// Precondition: `basis` must not contain duplicate Pauli words
+    /// (asserted in debug builds).
+    pub fn generator(&self, basis: &[Word]) -> Vec<(usize, usize, f64)> {
+        let index = build_basis_index(basis);
+        let n_qubits = self.n_qubits();
+
+        // `compute_action_terms` returns a deduplicated `Vec<(Word, f64)>`,
+        // so it can be scattered directly into COO triplets.
+        let local: Vec<Vec<(usize, usize, f64)>> = basis
+            .par_iter()
+            .enumerate()
+            .map_init(
+                || {
+                    (
+                        Vec::<u32>::with_capacity(n_qubits),
+                        Vec::<u32>::with_capacity(128),
+                        FxHashMap::<Word, Complex<f64>>::with_capacity_and_hasher(
+                            128,
+                            FxBuildHasher::default(),
+                        ),
+                    )
+                },
+                |(s1, s2, lm), (col, p)| {
+                    let terms = self.compute_action_terms(p, s1, s2, lm);
+                    let mut out = Vec::with_capacity(terms.len());
+                    for (w, v) in terms.iter() {
+                        if let Some(&row) = index.get(w) {
+                            out.push((row as usize, col, *v));
+                        }
+                    }
+                    out
+                },
+            )
+            .collect();
+
+        // Pre-allocate the flat output to avoid sequential push reallocation.
+        let total: usize = local.iter().map(|v| v.len()).sum();
+        let mut flat = Vec::with_capacity(total);
+        for v in local {
+            flat.extend(v);
+        }
+        flat
+    }
+
+    /// Complex-coefficient variant of [`Self::leakage`]: off-basis
+    /// component of `L*( Σ_j coeffs[j] · basis[j] )` with complex `coeffs`.
+    pub fn leakage_complex(
+        &self,
+        basis: &[Word],
+        coeffs: &[Complex<f64>],
+        protected: &[Word],
+    ) -> Result<Vec<(Word, Complex<f64>)>, Error> {
+        if basis.len() != coeffs.len() {
+            return Err(Error::LengthMismatch {
+                what: "basis and coeffs",
+                a: basis.len(),
+                b: coeffs.len(),
+            });
+        }
+        let in_basis: FxHashMap<u64, ()> = basis.iter().map(|w| (word_hash(w), ())).collect();
+        let protected_set: FxHashMap<u64, ()> =
+            protected.iter().map(|w| (word_hash(w), ())).collect();
+
+        const CHUNK_SIZE: usize = 4096;
+        let n_qubits = self.n_qubits();
+        let mut merged: FxHashMap<Word, Complex<f64>> = FxHashMap::default();
+        for chunk_start in (0..basis.len()).step_by(CHUNK_SIZE) {
+            let chunk_end = (chunk_start + CHUNK_SIZE).min(basis.len());
+            let chunk_basis = &basis[chunk_start..chunk_end];
+            let chunk_coeffs = &coeffs[chunk_start..chunk_end];
+            let local: Vec<Vec<(Word, Complex<f64>)>> = chunk_basis
+                .par_iter()
+                .zip(chunk_coeffs.par_iter())
+                .map_init(
+                    || {
+                        (
+                            Vec::<u32>::with_capacity(n_qubits),
+                            Vec::<u32>::with_capacity(128),
+                            FxHashMap::<Word, Complex<f64>>::with_capacity_and_hasher(
+                                128,
+                                FxBuildHasher::default(),
+                            ),
+                        )
+                    },
+                    |(s1, s2, lm), (p, &c)| {
+                        let terms = self.compute_action_terms(p, s1, s2, lm);
+                        let mut out = Vec::with_capacity(terms.len());
+                        for (w, v) in terms.iter() {
+                            let h = word_hash(w);
+                            if !in_basis.contains_key(&h) && !protected_set.contains_key(&h) {
+                                out.push((*w, c * *v));
+                            }
+                        }
+                        out
+                    },
+                )
+                .collect();
+            for v in local {
+                for (k, val) in v {
+                    *merged.entry(k).or_insert(Complex::new(0.0, 0.0)) += val;
+                }
+            }
+        }
+        Ok(merged.into_iter().filter(|(_, c)| c.norm() > 0.0).collect())
+    }
+}

crates/ppvm-lindblad/src/config.rs

@@ -4,12 +4,10 @@
 //! Configuration objects for the predictor-corrector stepper.
 
 /// Truncation and execution policy for a single predictor-corrector step
-/// ([`crate::LindbladSpec::pc_step`], [`crate::LindbladSpec::pc_step_timed`],
-/// [`crate::orbit_rep::pc_step_orbit_rep`]).
+/// ([`crate::LindbladSpec::pc_step`], [`crate::LindbladSpec::pc_step_timed`]).
 ///
 /// These are the per-run *tuning knobs*, kept separate from the per-call data
-/// (`basis`, `coeffs`, `dt`, `protected`, and — on the orbit path — the
-/// translation group and momentum).
+/// (`basis`, `coeffs`, `dt`, `protected`).
 ///
 /// `max_basis` is the primary accuracy/cost dial; `admit_basis` selects the
 /// displacement scheme; `drop_tol` is the churn valve of the admission-bound