wai/CHANGELOG.md

@@ -6,6 +6,33 @@ wire format may change and such changes are recorded here.
 
 ## v0.4.0 — current
 
+### wai-quantum 0.6.7: tensor networks that split controlled gates and fix projections
+
+Additive. `Network::amplitude` and `Network::simplify` are unchanged, so every pinned digest is too.
+
+- **`Network::amplitude_split`** (`quantum_tn`) splits each gate on two or more qubits at
+  every qubit it is block-diagonal on: a control, or a qubit it is diagonal on. There
+  `G = Σ_a |a⟩⟨a| ⊗ G_a`, so that qubit gets the delta `δ(out, in, k)` on its wire, and the
+  blocks `G_k` keep the rest, joined by a bond of dimension 2. Only zeros of the matrix are
+  read. A cut through a controlled-Z or controlled-X then costs 2, not 4.
+- **`Network::reduce`** fixes every index that some tensor allows only one value of: each
+  tensor holding it keeps its entries at that value and drops the index, until nothing
+  changes; then it simplifies. Only zeros are read and entries selected. Input and output
+  projections run through diagonal gates and split deltas, so a controlled-Z beside a
+  projected output becomes scalars. Open indices are never fixed.
+- On 61-qubit controlled-Z circuits of depth 8, 12, 16 and 20, the search (64 trials) finds
+  single-amplitude trees of log₁₀ 5.35, 8.25, 11.22 and 14.26 multiply-adds on split, reduced
+  networks, against 8.03, 12.46, 16.09 and 19.64 on whole-gate, simplified ones.
+- `examples/tn_cost`, the CLI's `wai-quantum tn` (`--whole-gates` for the old networks) and the
+  service's `tensor_network` method build split, reduced networks.
+- Tests: split networks, reduced or not, and reduced whole-gate networks give the state
+  vector's amplitudes on circuits of controlled-Z, controlled-X both ways, a controlled
+  random unitary, random diagonal gates, Toffoli, CCZ and random two-qubit unitaries, with
+  outputs projected and open. CZ and CX split and a random unitary does not; CX splits at its
+  control. On a 4 × 5 grid of controlled-Z, the split tree costs under an eighth of the
+  whole-gate one. Three mutants are caught: the delta's entry, the block-diagonal test, and
+  the value kept by `reduce`.
+
 ### wai-quantum 0.6.6: Clifford circuits measured in a product basis, exactly, by cut-rank
 
 Additive. Every pinned digest is unchanged.

wai/tools/quantum-meter/Cargo.lock

@@ -1017,7 +1017,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/tools/quantum-wgsl-check/Cargo.lock

@@ -599,7 +599,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-pack/Cargo.lock

@@ -713,7 +713,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum-qpu/Cargo.lock

@@ -422,7 +422,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum-serve/Cargo.lock

@@ -599,7 +599,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum-service/Cargo.lock

@@ -857,7 +857,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum-service/src/run.rs

@@ -370,8 +370,9 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
             for b in &r.amplitudes {
                 // qubit q's bit, in qubit order.
                 let bits: Vec<u8> = (0..n).map(|q| b.as_bytes()[(n - 1 - q) as usize] - b'0').collect();
-                let mut net = Network::amplitude(n, &gates, &bits, &[]).map_err(|e| format!("{e:?}"))?;
-                net.simplify();
+                // Controlled and diagonal gates split at their block-diagonal qubits.
+                let mut net = Network::amplitude_split(n, &gates, &bits, &[]).map_err(|e| format!("{e:?}"))?;
+                net.reduce();
                 let plan = find_path(&net.shape(), &net.dims, &net.open, &PathOptions::default());
                 flops += plan.total_flops();
                 let t = contract_plan(&net, &plan).map_err(|e| format!("{e:?}"))?;

wai/wai-quantum/Cargo.lock

@@ -1079,7 +1079,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "block2",

wai/wai-quantum/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 edition = "2024"
 rust-version = "1.97"
 license = "Apache-2.0"

wai/wai-quantum/README.md

@@ -424,6 +424,17 @@ wai-quantum tn --file circuit.jsonl --plan-only --cap-log2 30
 
 Circuits are JSON lines, one gate per line: `{"q":[a,b],"m":[re,im, …]}`.
 
+Each gate is split at the qubits it is block-diagonal on: a control, or a
+qubit it is diagonal on. That qubit's part is an exact three-index delta,
+joined to the rest of the gate by a bond of dimension 2, so a cut through a
+controlled-Z costs 2 instead of the 4 of its two wires. Then every index
+that some tensor allows only one value of is fixed there, so the input and
+output projections run through diagonal gates: a controlled-Z beside a
+projected output becomes scalars. On 61-qubit controlled-Z circuits of depth
+8 to 20, the trees found on split, reduced networks cost 2.7 to 5.4 decades
+less than on whole-gate ones. `--whole-gates` keeps every gate one tensor and
+the network as built.
+
 On single amplitudes of random circuits on the 54-qubit advantage-experiment
 grid, the search finds trees of `10¹⁰·⁵` operations at 10 cycles and `10¹⁸·⁸` at
 20, in seconds. An independently developed hyper-optimiser, given the same

wai/wai-quantum/examples/tn_cost.rs

@@ -8,7 +8,9 @@
 //! `CIRCUIT.qasm`: OpenQASM 3 with `rz`, `rx`, `ry`, `h`, `x`, `sx`, `cz`,
 //! `cx` (barriers and measurements skipped). Prints log10 of the complex
 //! multiply-adds for one amplitude (every slice), log2 of the largest
-//! intermediate, and the number of slices.
+//! intermediate, and the number of slices. The network splits each `cz` and
+//! `cx` at the qubits it is block-diagonal on ([`Network::amplitude_split`]);
+//! `TN_DENSE=1` keeps every gate whole.
 
 use std::time::Instant;
 use wai_quantum::gates::Gate;
@@ -50,8 +52,16 @@ fn main() {
         });
     }
     let bits = vec![0u8; n as usize];
-    let mut net = Network::amplitude(n, &gates, &bits, &[]).expect("the network builds");
-    net.simplify();
+    // TN_DENSE=1: every gate one tensor, as before splitting.
+    let dense = std::env::var("TN_DENSE").as_deref() == Ok("1");
+    let mut net = if dense { Network::amplitude(n, &gates, &bits, &[]) } else { Network::amplitude_split(n, &gates, &bits, &[]) }.expect("the network builds");
+    // Split networks are reduced (projections fixed through diagonal gates);
+    // TN_DENSE=1 is the pipeline before both.
+    if dense {
+        net.simplify();
+    } else {
+        net.reduce();
+    }
     let t = Instant::now();
     let plan = find_path(&net.shape(), &net.dims, &net.open, &PathOptions { trials, seed, max_size, ..PathOptions::default() });
     println!(

wai/wai-quantum/src/quantum_tn.rs

@@ -14,9 +14,17 @@
 //! - **Build.** [`Network::amplitude`] turns gates (any unitary on any qubits,
 //!   given as a matrix) and an output bitstring into a network. Any subset of
 //!   the outputs can be left open, giving a batch of amplitudes at once.
+//!   [`Network::amplitude_split`] also splits each gate at the qubits it is
+//!   block-diagonal on, a control or a qubit it is diagonal on. Each such
+//!   qubit gets an exact three-index delta joined to the rest by a bond of
+//!   dimension 2, so a cut through a controlled-Z or controlled-X costs 2
+//!   instead of the 4 of its two wires.
 //! - **Simplify.** Vectors and one-qubit gates are absorbed into their
 //!   neighbours before any search. They never cost more than the tensor they
-//!   join.
+//!   join. [`Network::reduce`] first fixes every index some tensor allows only
+//!   one value of, so input and output projections run through diagonal gates
+//!   and split deltas; a controlled-Z beside a projected output becomes
+//!   scalars.
 //! - **Order.** [`find_path`] runs a seeded hyper-optimiser over two
 //!   families of contraction trees:
 //!   - randomised greedy, scoring a pair by
@@ -81,7 +89,12 @@
 //!   second on large products. Plans worth `10¹⁴` operations or more are for
 //!   planning and verification arithmetic here, not for contraction.
 //! - **Networks are simple graphs.** Every index joins at most two tensors.
-//!   Diagonal-gate hyperedges are not exploited.
+//!   Splitting gives a controlled or diagonal gate's qubit a three-index
+//!   delta, not a hyperedge: two deltas in a row on one wire contract to a
+//!   four-index delta, held densely, where a hyperedge would stay one index.
+//!   On 61-qubit controlled-Z circuits of depth 8 to 20, split and reduced
+//!   networks cost 2.7 to 5.4 decades less than whole-gate, simplified ones,
+//!   with the same search.
 
 pub use crate::linalg::C;
 use crate::repro::ln;
@@ -122,6 +135,23 @@ impl Network {
     /// `2^|open|` amplitudes over those qubits (first listed most
     /// significant), with the other outputs fixed to `bits`.
     pub fn amplitude(n: u32, gates: &[Gate], bits: &[u8], open_qubits: &[u32]) -> Result<Network, TnError> {
+        Self::build(n, gates, bits, open_qubits, false)
+    }
+
+    /// [`Network::amplitude`], with each gate on two or more qubits split at
+    /// every qubit it is block-diagonal on: a qubit it controls, or on which it
+    /// is diagonal. There `G = Σ_a |a⟩⟨a| ⊗ G_a`, so the qubit's part is the
+    /// three-index delta `δ(out, in, k)` and the rest is `G_k`, joined by a bond
+    /// `k` of dimension 2 (at least one qubit stays in the rest). The same
+    /// amplitudes, by the same entries, with no arithmetic: only zeros of the
+    /// matrix are read. A cut through such a gate then costs the bond's 2,
+    /// not the 4 of its two wires, which is most of what a contraction of
+    /// controlled-Z or controlled-X circuits costs.
+    pub fn amplitude_split(n: u32, gates: &[Gate], bits: &[u8], open_qubits: &[u32]) -> Result<Network, TnError> {
+        Self::build(n, gates, bits, open_qubits, true)
+    }
+
+    fn build(n: u32, gates: &[Gate], bits: &[u8], open_qubits: &[u32], split: bool) -> Result<Network, TnError> {
         if bits.len() != n as usize || bits.iter().any(|&b| b > 1) || open_qubits.iter().any(|&q| q >= n) {
             return Err(TnError::BadOutput);
         }
@@ -143,6 +173,18 @@ impl Network {
             if k == 0 || seen.len() != k || g.qubits.iter().any(|&q| q >= n) || g.matrix.len() != 1 << (2 * k) {
                 return Err(TnError::BadGate(gi));
             }
+            if split && k >= 2 {
+                let (bonds, rest, data) = peel(&g.qubits, &g.matrix, &mut fresh, &mut wire, &mut tensors);
+                let outs: Vec<u32> = rest.iter().map(|_| fresh()).collect();
+                let mut inds = bonds;
+                inds.extend(&outs);
+                inds.extend(rest.iter().map(|&q| wire[q as usize]));
+                tensors.push(Tensor { inds, data });
+                for (i, &q) in rest.iter().enumerate() {
+                    wire[q as usize] = outs[i];
+                }
+                continue;
+            }
             let outs: Vec<u32> = (0..k).map(|_| fresh()).collect();
             let mut inds = outs.clone();
             inds.extend(g.qubits.iter().map(|&q| wire[q as usize]));
@@ -208,6 +250,68 @@ impl Network {
         }
     }
 
+    /// [`Network::simplify`], after fixing every index that some tensor
+    /// allows one value of: where a tensor is zero unless index `i` takes `v`,
+    /// every term with `i ≠ v` vanishes, so each tensor holding `i` keeps only
+    /// its entries at `v` and drops `i`. Repeated until nothing changes. Only
+    /// zeros are read and entries selected, so the network's value is the same
+    /// sum. Input and output projections then run through diagonal gates and
+    /// the deltas of split ones: a controlled-Z next to a projected output
+    /// collapses to scalars. Open indices are never fixed. Tensors left with
+    /// no index are multiplied into one of the others.
+    pub fn reduce(&mut self) {
+        loop {
+            let mut changed = false;
+            'scan: for t in 0..self.tensors.len() {
+                for p in 0..self.tensors[t].inds.len() {
+                    let i = self.tensors[t].inds[p];
+                    if self.open.contains(&i) {
+                        continue;
+                    }
+                    let dims: Vec<usize> = self.tensors[t].inds.iter().map(|&j| self.dims[j as usize]).collect();
+                    let inner: usize = dims[p + 1..].iter().product();
+                    let d = dims[p];
+                    let mut values = (0..d).filter(|&v| self.tensors[t].data.iter().enumerate().any(|(x, c)| (x / inner) % d == v && (c.re != 0.0 || c.im != 0.0)));
+                    let (Some(v), None) = (values.next(), values.next()) else { continue };
+                    for u in &mut self.tensors {
+                        if let Some(q) = u.inds.iter().position(|&j| j == i) {
+                            let udims: Vec<usize> = u.inds.iter().map(|&j| self.dims[j as usize]).collect();
+                            let uin: usize = udims[q + 1..].iter().product();
+                            let ud = udims[q];
+                            u.data = u.data.iter().enumerate().filter(|(x, _)| (x / uin) % ud == v).map(|(_, c)| *c).collect();
+                            u.inds.remove(q);
+                        }
+                    }
+                    changed = true;
+                    break 'scan;
+                }
+            }
+            if !changed {
+                break;
+            }
+        }
+        // Scalars: each into the first tensor that has an index (or, if none
+        // has, into the first).
+        let scalars: Vec<usize> = (0..self.tensors.len()).filter(|&t| self.tensors[t].inds.is_empty()).collect();
+        if !scalars.is_empty() && self.tensors.len() > 1 {
+            let host = (0..self.tensors.len()).find(|&t| !self.tensors[t].inds.is_empty()).unwrap_or(scalars[0]);
+            let mut factor = C::ONE;
+            for &t in scalars.iter().filter(|&&t| t != host) {
+                factor = factor.mul(self.tensors[t].data[0]);
+            }
+            for c in &mut self.tensors[host].data {
+                *c = c.mul(factor);
+            }
+            let mut k = 0;
+            self.tensors.retain(|_| {
+                let keep = !scalars.contains(&k) || k == host;
+                k += 1;
+                keep
+            });
+        }
+        self.simplify();
+    }
+
     fn owners(&self) -> BTreeMap<u32, Vec<usize>> {
         let mut o: BTreeMap<u32, Vec<usize>> = BTreeMap::new();
         for (t, x) in self.tensors.iter().enumerate() {
@@ -224,6 +328,58 @@ impl Network {
     }
 }
 
+/// Split a gate at each of its qubits it is block-diagonal on, while two or
+/// more remain: each such qubit gets the delta `δ(out, in, k)` on its wire (its
+/// wire moving to `out`), and the rest keeps the blocks, indexed by the new
+/// bonds. Returns the bonds (most significant first), the qubits left, and the
+/// rest's entries over `[bonds…, outs…, ins…]`.
+fn peel(qubits: &[u32], matrix: &[C], fresh: &mut impl FnMut() -> u32, wire: &mut [u32], tensors: &mut Vec<Tensor>) -> (Vec<u32>, Vec<u32>, Vec<C>) {
+    let mut bonds: Vec<u32> = Vec::new();
+    let mut rest: Vec<u32> = qubits.to_vec();
+    let mut data: Vec<C> = matrix.to_vec();
+    while rest.len() >= 2 {
+        let m = rest.len();
+        let side = 1usize << m;
+        // Block-diagonal on rest[j]: no entry, for any bonds, has that qubit's
+        // output bit differ from its input bit.
+        let found = (0..m).find(|&j| {
+            let bit = 1usize << (m - 1 - j);
+            data.iter().enumerate().all(|(x, c)| {
+                let (row, col) = ((x / side) % side, x % side);
+                (row & bit) == (col & bit) || (c.re == 0.0 && c.im == 0.0)
+            })
+        });
+        let Some(j) = found else { break };
+        let q = rest[j];
+        let (k, out) = (fresh(), fresh());
+        let mut delta = vec![C::ZERO; 8];
+        delta[0] = C::ONE; // out 0, in 0, k 0
+        delta[7] = C::ONE; // out 1, in 1, k 1
+        tensors.push(Tensor { inds: vec![out, wire[q as usize], k], data: delta });
+        wire[q as usize] = out;
+        // The blocks: the entries with that qubit's bits both equal to k, the
+        // new bond after the old ones.
+        let half = side / 2;
+        let bit = m - 1 - j;
+        let insert = |v: usize, b: usize| ((v >> bit) << (bit + 1)) | (b << bit) | (v & ((1 << bit) - 1));
+        let groups = data.len() / (side * side);
+        let mut next = vec![C::ZERO; groups * 2 * half * half];
+        for g in 0..groups {
+            for b in 0..2 {
+                for r in 0..half {
+                    for c in 0..half {
+                        next[((g * 2 + b) * half + r) * half + c] = data[(g * side + insert(r, b)) * side + insert(c, b)];
+                    }
+                }
+            }
+        }
+        data = next;
+        bonds.push(k);
+        rest.remove(j);
+    }
+    (bonds, rest, data)
+}
+
 /// Indices in exactly one of `a`, `b`, plus open indices in either.
 fn sym_diff(a: &[u32], b: &[u32], open: &[u32]) -> Vec<u32> {
     let mut out: Vec<u32> = a.iter().copied().filter(|i| !b.contains(i) || open.contains(i)).collect();
@@ -1537,6 +1693,137 @@ pub(crate) mod tests {
         }
     }
 
+    /// Gates that split and gates that do not: controlled-Z and controlled-X
+    /// both ways, a controlled random unitary, a diagonal gate at random
+    /// phases, Toffoli and controlled-controlled-Z, and random two-qubit
+    /// unitaries, between random one-qubit layers.
+    fn mixed_circuit(n: u32, depth: usize, seed: u64) -> Vec<Gate> {
+        let mut rng = Rng::new(seed);
+        let (o, z) = (C::ONE, C::ZERO);
+        let mut gates = Vec::new();
+        for _ in 0..depth {
+            for q in 0..n {
+                gates.push(Gate { qubits: vec![q], matrix: random_unitary(1, &mut rng) });
+            }
+            for _ in 0..n / 2 {
+                let a = rng.below(n as usize) as u32;
+                let b = (a + 1 + rng.below(n as usize - 1) as u32) % n;
+                let c = (0..n).find(|&c| c != a && c != b).unwrap();
+                let u = random_unitary(1, &mut rng);
+                let ph: Vec<C> = (0..4).map(|_| { let t = 6.28 * rng.unit(); C::new(t.cos(), t.sin()) }).collect();
+                let matrix = match rng.below(7) {
+                    0 => vec![o, z, z, z, z, o, z, z, z, z, o, z, z, z, z, C::new(-1.0, 0.0)],
+                    1 => vec![o, z, z, z, z, o, z, z, z, z, z, o, z, z, o, z],
+                    2 => vec![o, z, z, z, z, o, z, z, z, z, u[0], u[1], z, z, u[2], u[3]],
+                    3 => vec![ph[0], z, z, z, z, ph[1], z, z, z, z, ph[2], z, z, z, z, ph[3]],
+                    4 => {
+                        // Toffoli on (a, b, c), c the target.
+                        let mut m = vec![z; 64];
+                        for x in 0..8 {
+                            let y = if x >= 6 { x ^ 1 } else { x };
+                            m[y * 8 + x] = o;
+                        }
+                        gates.push(Gate { qubits: vec![a, b, c], matrix: m });
+                        continue;
+                    }
+                    5 => {
+                        let mut m = vec![z; 64];
+                        for x in 0..8 {
+                            m[x * 8 + x] = if x == 7 { C::new(-1.0, 0.0) } else { o };
+                        }
+                        gates.push(Gate { qubits: vec![c, a, b], matrix: m });
+                        continue;
+                    }
+                    _ => random_unitary(2, &mut rng),
+                };
+                gates.push(Gate { qubits: vec![a, b], matrix });
+            }
+        }
+        gates
+    }
+
+    #[test]
+    fn split_networks_give_the_same_amplitudes() {
+        for seed in 0..8 {
+            let n = 7;
+            let gates = mixed_circuit(n, 5, seed);
+            let psi = dense(n, &gates);
+            for x in [0usize, 9, 64, 127] {
+                let bits: Vec<u8> = (0..n).map(|q| ((x >> (n - 1 - q)) & 1) as u8).collect();
+                for (split, reduce) in [(true, false), (true, true), (false, true)] {
+                    let mut net = if split { Network::amplitude_split(n, &gates, &bits, &[]) } else { Network::amplitude(n, &gates, &bits, &[]) }.unwrap();
+                    if reduce {
+                        net.reduce();
+                    } else {
+                        net.simplify();
+                    }
+                    let path = greedy_path(&net.shape(), &net.dims, &net.open, 1.0, 0.2, seed);
+                    let t = contract_path(&net, &path).unwrap();
+                    assert_eq!(t.data.len(), 1);
+                    assert!(close(t.data[0], psi[x]), "seed {seed} x {x} split {split} reduce {reduce}: {:?} vs {:?}", t.data[0], psi[x]);
+                }
+            }
+            // Open outputs: qubits 6 and 2 open (6 most significant), the rest 0.
+            let mut net = Network::amplitude_split(n, &gates, &[0; 7], &[6, 2]).unwrap();
+            net.reduce();
+            let t = contract_path(&net, &greedy_path(&net.shape(), &net.dims, &net.open, 1.0, 0.0, 0)).unwrap();
+            for (j, got) in t.data.iter().enumerate() {
+                // j = b6·2 + b2; qubit q sits at bit n − 1 − q of x.
+                let x = ((j >> 1) & 1) | ((j & 1) << 4);
+                assert!(close(*got, psi[x]), "seed {seed} open j {j}");
+            }
+        }
+    }
+
+    #[test]
+    fn controlled_and_diagonal_gates_split_and_others_do_not() {
+        let rank_counts = |g: Gate| {
+            let net = Network::amplitude_split(2, &[g], &[0, 0], &[]).unwrap();
+            // Two input vectors and two output vectors around the gate's parts.
+            let mut ranks: Vec<usize> = net.tensors[2..net.tensors.len() - 2].iter().map(|t| t.inds.len()).collect();
+            ranks.sort_unstable();
+            ranks
+        };
+        assert_eq!(rank_counts(Gate::cz(0, 1)), vec![3, 3]);
+        assert_eq!(rank_counts(Gate::cx(0, 1)), vec![3, 3]);
+        let mut rng = Rng::new(1);
+        assert_eq!(rank_counts(Gate { qubits: vec![0, 1], matrix: random_unitary(2, &mut rng) }), vec![4]);
+        // CX splits at its control: the delta sits on qubit 0's wire.
+        let net = Network::amplitude_split(2, &[Gate::cx(0, 1)], &[0, 0], &[]).unwrap();
+        let delta = &net.tensors[2];
+        assert_eq!((delta.inds[1], delta.data.iter().filter(|c| c.re == 1.0).count()), (0, 2));
+    }
+
+    #[test]
+    fn splitting_lowers_the_cost_of_a_controlled_z_circuit() {
+        // A 4 × 5 grid: random one-qubit layers between layers of CZ.
+        let (w, h) = (5u32, 4u32);
+        let n = w * h;
+        let mut rng = Rng::new(7);
+        let mut gates = Vec::new();
+        for layer in 0..12u32 {
+            for q in 0..n {
+                gates.push(Gate { qubits: vec![q], matrix: random_unitary(1, &mut rng) });
+            }
+            for y in 0..h {
+                for x in 0..w {
+                    let (a, b, ok) = if layer % 2 == 0 { (y * w + x, y * w + x + 1, x + 1 < w && (x + layer / 2) % 2 == 0) } else { (y * w + x, (y + 1) * w + x, y + 1 < h && (y + layer / 2) % 2 == 0) };
+                    if ok {
+                        gates.push(Gate::cz(a, b));
+                    }
+                }
+            }
+        }
+        let bits = vec![0u8; n as usize];
+        let cost = |mut net: Network| {
+            net.simplify();
+            find_path(&net.shape(), &net.dims, &net.open, &PathOptions { trials: 16, seed: 3, max_size: None, reconfigure: 4, families: Families::Both }).total_flops()
+        };
+        let dense_cost = cost(Network::amplitude(n, &gates, &bits, &[]).unwrap());
+        let split_cost = cost(Network::amplitude_split(n, &gates, &bits, &[]).unwrap());
+        assert!(split_cost * 8.0 < dense_cost, "split {split_cost:.3e} against dense {dense_cost:.3e}");
+    }
+
     #[test]
     fn open_outputs_give_a_batch_of_amplitudes() {
         let n = 6;

wai/wai-rs/Cargo.lock

@@ -4310,7 +4310,7 @@ version = "0.1.0"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-rs/src/bin/wai_quantum.rs

@@ -207,7 +207,7 @@ fn usage() -> ! {
          STATE VECTOR (any gate, any angle):\n  \
          wai-quantum sv --file CIRCUIT.jsonl [-n N] [--fuse W] [--top K | --sample S [--seed X] | --expect \"Z0 Z1; X3\"] [--json]\n\n\
          TENSOR NETWORKS:\n  \
-         wai-quantum tn --file CIRCUIT.jsonl [-n N] [--bits B] [--open q,q] [--trials T] [--seed S] [--cap-log2 L] [--plan-only] [--json]\n\n\
+         wai-quantum tn --file CIRCUIT.jsonl [-n N] [--bits B] [--open q,q] [--trials T] [--seed S] [--cap-log2 L] [--whole-gates] [--plan-only] [--json]\n\n\
          Channels: identity, bitflip, dephasing, depolarizing, damping.\n\
          Run `wai-quantum list` for the algorithm catalog."
     );
@@ -884,9 +884,18 @@ fn main() {
             if open.len() > 20 {
                 die("at most 20 open qubits (2^20 amplitudes)");
             }
-            let mut net = Network::amplitude(n, &gates, &bits, &open).unwrap_or_else(|e| die(&format!("{e:?}")));
+            // Gates split where they are block-diagonal (a control, or diagonal
+            // on a qubit), unless --whole-gates keeps each gate one tensor.
+            let whole = args.iter().any(|a| a == "--whole-gates");
+            let built = if whole { Network::amplitude(n, &gates, &bits, &open) } else { Network::amplitude_split(n, &gates, &bits, &open) };
+            let mut net = built.unwrap_or_else(|e| die(&format!("{e:?}")));
             let raw = net.tensors.len();
-            net.simplify();
+            // Projections fixed through diagonal gates, then absorption.
+            if whole {
+                net.simplify();
+            } else {
+                net.reduce();
+            }
             let shape = net.shape();
             let cap: Option<f64> = opt(args, "--cap-log2", "").map(|v| 2f64.powf(v.parse().unwrap_or_else(|_| die("bad --cap-log2"))));
             let options = PathOptions {

wai/wai-web/Cargo.lock

@@ -1241,7 +1241,7 @@ version = "0.1.0"
 
 [[package]]
 name = "wai-quantum"
-version = "0.6.6"
+version = "0.6.7"
 dependencies = [
  "blake3",
  "ed25519-dalek",