Symmetry-merged evolution: Trotter merging + momentum-sector CTPP

crates/ppvm-lindblad/src/lib.rs

@@ -47,6 +47,10 @@ pub(crate) mod expm;
 /// Matrix-free / quspin-expm-backed `exp(dt·L*)·b` engine. See module docs.
 pub(crate) mod mf_expm;
 
+/// Per-step orbit-rep evolution under translation symmetry, with a
+/// phase-aware complex action. See module docs.
+pub mod orbit_rep;
+
 /// Pauli-word storage chunk: `u64` on 64-bit targets, `u32` elsewhere
 /// (`bitvec` implements `BitStore` for `u64` only on 64-bit targets, so
 /// e.g. wasm32 builds use four 32-bit chunks instead of two 64-bit ones).
@@ -1135,6 +1139,109 @@ mod tests {
         assert_eq!(out.as_slice(), &codes);
     }
 
+    /// Per-step orbit-rep evolution gives the SAME final orbit-rep
+    /// state as full-basis complex evolution followed by a single
+    /// projection at the end. Validates that the phase-aware complex
+    /// action machinery is consistent with the full-basis reference.
+    #[test]
+    fn pc_step_orbit_rep_matches_full_basis_projection() {
+        use ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex;
+        use std::f64::consts::PI;
+        let n = 4usize;
+        let dt = 0.01f64;
+        let n_steps = 3usize;
+        let mut h_terms: Vec<(String, f64)> = Vec::new();
+        for j in 0..n {
+            let nxt = (j + 1) % n;
+            for op in ["X", "Y"] {
+                let mut s = vec!['I'; n];
+                s[j] = op.chars().next().unwrap();
+                s[nxt] = op.chars().next().unwrap();
+                h_terms.push((s.into_iter().collect(), 1.0));
+            }
+        }
+        let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap();
+        let group = ppvm_pauli_sum::symmetry::TranslationGroup::chain_1d(n);
+        let k_mode: i32 = 1;
+        let k = vec![k_mode];
+
+        // Build the k=1 eigenstate in FULL basis form.
+        let basis_full: Vec<Word> = (0..n)
+            .map(|j| {
+                let mut s = vec!['I'; n];
+                s[j] = 'Z';
+                let (w, _) = parse_pauli_string(&s.into_iter().collect::<String>(), n).unwrap();
+                w
+            })
+            .collect();
+        let coeffs_full: Vec<Complex<f64>> = (0..n as i32)
+            .map(|a| {
+                Complex::from_polar(1.0, -2.0 * PI * (k_mode as f64) * (a as f64) / (n as f64))
+            })
+            .collect();
+
+        // ----- Full-basis path -----
+        let mut bf = basis_full.clone();
+        let mut cf = coeffs_full.clone();
+        let protected: Vec<Word> = Vec::new();
+        for _ in 0..n_steps {
+            // Full enrichment (tau_add = 0.0 adds every leakage string):
+            // for a momentum eigenstate the leakage is pure-sector, so the
+            // full-basis and orbit-rep paths build corresponding bases and
+            // the projection theorem gives an exact match. The orbit-rep
+            // side uses a large max_basis so its rank cap never binds.
+            pc_step_complex_full(&spec, &mut bf, &mut cf, dt);
+        }
+        // Project at the end.
+        canonicalize_pauli_sum_complex(&mut bf, &mut cf, &group, &k);
+
+        // ----- Orbit-rep path -----
+        // Initial orbit-rep form: project the full-basis input.
+        let mut br = basis_full.clone();
+        let mut cr = coeffs_full.clone();
+        canonicalize_pauli_sum_complex(&mut br, &mut cr, &group, &k);
+        // Evolve in orbit-rep form (max_basis large ⇒ full enrichment).
+        for _ in 0..n_steps {
+            orbit_rep::pc_step_orbit_rep(
+                &spec,
+                &mut br,
+                &mut cr,
+                dt,
+                &protected,
+                &group,
+                &k,
+                &PcStepConfig {
+                    max_basis: 10_000_000,
+                    ..Default::default()
+                },
+            )
+            .unwrap();
+        }
+
+        // Compare.
+        let mf: FxHashMap<Word, Complex<f64>> = bf.into_iter().zip(cf).collect();
+        let mr: FxHashMap<Word, Complex<f64>> = br.into_iter().zip(cr).collect();
+        assert_eq!(
+            mf.len(),
+            mr.len(),
+            "orbit-rep ({}) and full-basis-projected ({}) basis sizes differ",
+            mr.len(),
+            mf.len()
+        );
+        let mut max_diff = 0.0_f64;
+        for (w, cm) in &mr {
+            let cf_val = mf
+                .get(w)
+                .copied()
+                .unwrap_or_else(|| panic!("rep {:?} in orbit-rep but not in full-basis", w));
+            max_diff = max_diff.max((cm - cf_val).norm());
+        }
+        assert!(
+            max_diff < 1e-9,
+            "orbit-rep diverged from full-basis: max |Δc| = {max_diff:e}"
+        );
+    }
+
     /// The full-space complex step at momentum k=0 must reproduce the real
     /// pc_step on the same trajectory exactly.
     #[test]