//! `wai-quantum lab`: the research tools of the command line, the deterministic //! quantum stack's every capability. //! //! Every quantum capability in this crate is pure Rust with no host bindings, so //! this binary targets **native and cloud WASM from one source**: //! `cargo build --release --target wasm32-wasip2` produces a component that //! `wasmtime` runs with **byte-identical output** — the same circuit reconstructs //! to the same statevector hash on a laptop or at the edge. That equality is the //! point: a quantum computation you can put in a receipt and re-run anywhere. //! The one exception is `fabric`, which describes the host it runs on: its //! certificates name the member they ran on, and agree everywhere in parity //! and corpus digest. //! //! ```text //! wai-quantum lab list //! wai-quantum lab run ghz --qubits 5 [--top 8] [--json] //! wai-quantum lab run --qasm # read OpenQASM on stdin //! wai-quantum lab qasm bell # emit OpenQASM //! wai-quantum lab receipt ghz -n 3 [--shots 4096 --seed 7] //! wai-quantum lab verify # read a receipt JSON on stdin //! ``` //! //! A circuit is built in `wai.quantum.circuit2`, whose phase table is correctly //! rounded. `--gateset v1` builds it in the frozen `wai.quantum.circuit` instead, //! to reproduce or re-verify content made with it. Every command that builds a //! circuit from an algorithm, OpenQASM or QIR (`run`, `receipt`, `qasm`, `qir`, //! `entropy`) takes it, and so does `fabric`. //! //! Beyond running circuits, the same binary answers what a state *carries* and //! what transport does to it: //! //! ```text //! wai-quantum lab entropy bell --keep 0 # half a Bell pair: 1 bit, and no state vector //! wai-quantum lab channel depolarizing --prob 0.5 //! wai-quantum lab bb84 --eavesdrop # interception shows up in the error rate //! wai-quantum lab chsh # 2√2, exactly, not sampled //! wai-quantum lab qutrit # d = 3, not just qubits //! ``` //! //! At utility scale, sparse Pauli dynamics computes observables of circuits far //! past any state vector, with arbitrary rotation angles: //! //! ```text //! wai-quantum lab utility --theta 0.2 --steps 20 # the 127-qubit heavy-hex experiment //! wai-quantum lab utility --method bptn --chi 32 # the same, by tensor networks //! wai-quantum lab utility --receipt | wai-quantum verify //! ``` //! //! And for fault-tolerant machines that do not exist yet, what a program would //! need — physical qubits and time — on six reference qubits, as a re-checkable //! signed estimate: //! //! ```text //! wai-quantum lab estimate factoring # 2048-bit factoring, all six qubits //! wai-quantum lab estimate chemistry --qubit ns-e4 --frontier //! wai-quantum lab estimate --qubits 100 --ccz 1e6 --budget 0.01 //! wai-quantum lab estimate factoring --qubit ns-e4 --receipt | wai-quantum verify //! ``` //! //! Error correction itself, at circuit level: sample a noisy Clifford circuit, //! derive its detector error model, decode it exactly, and measure the code //! model that estimates otherwise take on trust: //! //! ```text //! wai-quantum lab qec -d 5 --p 0.003 --shots 100000 //! wai-quantum lab qec --file circuit.txt //! wai-quantum lab qec --fit # then: estimate --prefactor A --threshold P //! ``` //! //! The `channel` command prints the Holevo bound before and after, which can only //! fall: noise never creates information. Two results there are easy to miss — //! fully damping `|1⟩` yields a *pure* output that still carries nothing (both //! inputs land on `|0⟩`), and dephasing leaves a Z-basis ensemble's Holevo bound //! untouched, because loss is relative to the encoding rather than absolute. use wai_quantum::ed25519_dalek::SigningKey; use std::io::Read; use wai_quantum::quantum::{BaseGate, Circuit, Gateset, Measure, ONE}; use wai_quantum::quantum_channel::{tomography, Channel}; use wai_quantum::quantum_comm::{bb84, clone_fidelity}; use wai_quantum::quantum_info::{chsh, DensityMatrix}; use wai_quantum::quantum_qudit::QuditState; use wai_quantum::quantum_energy::{EnergyClass, MaybeLabelled}; use wai_quantum::quantum_receipt::QuantumReceipt; use wai_quantum::quantum_receipt::did_key; use wai_quantum::quantum_resource::{self as qr, EstimateConfig, LogicalCounts, PhysicalQubit, QecModel}; use wai_quantum::quantum_resource_receipt::{Program, ResourceJob, ResourceReceipt}; use wai_quantum::quantum_spd_receipt::{SpdJob, SpdReceipt}; use wai_quantum::quantum_source::{holevo_bound, schumacher_limit, von_neumann_entropy}; use wai_quantum::quantum_toolchain as qt; /// The three single-qubit states the channel command sends: the ground state, the /// excited one, and the superposition whose coherence noise actually destroys. fn one_qubit_state(which: &str) -> DensityMatrix { let mut c = Circuit::new(1); match which { "zero" => {} "one" => { c.x(0); } "plus" => { c.h(0); } other => die(&format!("unknown state '{other}' (zero, one, plus)")), } DensityMatrix::from_pure(&c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}")))) } /// The joule figure as a line shows it: `0` is the unmetered marker, shown as /// unmetered and never as `0 uJ` (energy-measurement §2). fn joules_line(joules_micro: u64) -> String { if joules_micro == 0 { "unmetered (no figure; 0 is the unmetered marker, not zero joules)".into() } else { format!("{joules_micro} uJ") } } /// One line saying what the joule figure is: its class and evidence, or that /// it has none. fn energy_class_line(class: Option<&EnergyClass>, joules_micro: u64) -> String { let Some(c) = class else { return if joules_micro == 0 { "none (unmetered)".into() } else { "unlabelled (legacy form: read as no class, not comparable with a labelled figure)".into() }; }; let mut line = c.label().to_string(); if let Some(u) = c.uncertainty() { line += &format!( " (declared u {:.2} %, k=1, over {:.3} s", u.relative_ppm as f64 / 1e4, u.window_us as f64 / 1e6 ); match c { EnergyClass::OnChipCounter { filtering, .. } => line += &format!(", filtering {})", filtering.label()), EnergyClass::CalibratedInstrument { calibration_ref, .. } => line += &format!(", calibration {calibration_ref})"), _ => line += ")", } } line } fn hx(b: &[u8]) -> String { b.iter().map(|x| format!("{x:02x}")).collect() } fn die(msg: &str) -> ! { eprintln!("wai-quantum: {msg}"); std::process::exit(2); } fn usage() -> ! { eprintln!( "wai-quantum lab — the research tools (run `wai-quantum --help` for run, plan, lint, draw, verify and the rest)\n\n\ USAGE:\n \ wai-quantum lab list\n \ wai-quantum lab run [-n QUBITS] [-p PARAM] [--top K] [--json] [--gateset v1|v2]\n \ wai-quantum lab run --qasm [--top K] [--json] [--gateset v1|v2] (OpenQASM on stdin)\n \ wai-quantum lab run --qir [--top K] [--json] [--gateset v1|v2] (QIR on stdin)\n \ wai-quantum lab qasm [-n QUBITS] [-p PARAM] [--gateset v1|v2]\n \ wai-quantum lab qir [-n QUBITS] [-p PARAM] [--no-measure] [--gateset v1|v2]\n \ wai-quantum lab receipt [-n QUBITS] [-p PARAM] [--shots N] [--seed S] [--gateset v1|v2] [--key FILE | --demo-key]\n \ wai-quantum lab verify (receipt JSON on stdin)\n \ --gateset v2 (the default) is wai.quantum.circuit2; v1 is the frozen wai.quantum.circuit.\n\n\ INFORMATION AND TRANSPORT:\n \ wai-quantum lab entropy [-n QUBITS] [--keep 0,1] [--json] [--gateset v1|v2]\n \ wai-quantum lab channel [--prob P] [--state zero|one|plus] [--json]\n \ wai-quantum lab bb84 [--rounds N] [--seed S] [--eavesdrop] [--json]\n \ wai-quantum lab chsh [--json]\n \ wai-quantum lab qutrit [-n QUTRITS] [--json]\n\n\ UTILITY SCALE:\n \ wai-quantum lab utility [--theta T] [--steps N] [--threshold D] [-q QUBIT] [--json] [--receipt]\n \ wai-quantum lab utility --method bptn [--chi X] [--theta T] [--steps N] [-q QUBIT] [--json]\n\n\ FAULT-TOLERANT RESOURCES:\n \ wai-quantum lab estimate [factoring|chemistry] [--qubit NAME|all] [--code surface|surface-meas|floquet]\n \ [--budget E] [--slowdown N] [--frontier] [--json] [--receipt]\n \ wai-quantum lab estimate --qubits N [--ccz N] [--t N] [--rotations N] [--depth N] [--measurements N] ...\n \ Qubits: us-e3, us-e4, ns-e3, ns-e4, meas-e4, meas-e6. A measured code: --prefactor A --threshold P.\n \ From a program's own counts, with the code model fitted to the decoders: wai-quantum estimate prog.qasm\n\n\ ERROR CORRECTION AT CIRCUIT LEVEL:\n \ wai-quantum lab qec [-d D] [--rounds R] [--p P] [--shots N] [--seed S] [--decoder exact|uf] [--json]\n \ wai-quantum lab qec --file CIRCUIT [--shots N] (the ecosystem's circuit text)\n \ wai-quantum lab qec [--circuit | --dem] [-d D | --file CIRCUIT]\n \ wai-quantum lab qec --fit [--shots N] (measure a and p* of the code model)\n\n\ GAUSSIAN BOSON SAMPLING:\n \ wai-quantum lab gbs --squeezing r1,r2,… [--eta E] [--seed S | --unitary U.json] [--json]\n \ wai-quantum lab gbs … --sample K (exact click patterns)\n \ wai-quantum lab gbs … --check SAMPLES.txt (observed patterns vs exact marginals)\n\n\ SPIN CHAINS (DMRG ground states, exact referees):\n \ wai-quantum lab dmrg [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--bond D] [--sweeps S]\n \ [--cutoff C] [--exact] [--json]\n \ wai-quantum lab tdvp [--chain ising|heisenberg] [-n N] [--j J] [--jz JZ] [--h H] [--start plus|neel|up]\n \ [--bond D] [--dt DT] [--t T] [--every K] [--cutoff C] [--json] (time evolution by TDVP)\n\n\ MOLECULES (DMRG from FCIDUMP integrals, full CI as referee):\n \ wai-quantum lab chem --fcidump FILE [--bond D] [--sweeps S] [--penalty P] [--no-fci] [--json]\n\n\ NEURAL QUANTUM STATES (RBM, variational Monte Carlo, stochastic reconfiguration):\n \ wai-quantum lab nqs [--model heisenberg|j1j2|ising] [--lx L] [--ly L] [--j2 J] [--h H] [--periodic]\n \ [--complex] [--project [--kpi] [--point-group a1|a2|b1|b2 | --reflections a1|a2|b1|b2]\n \ [--odd]] [--alpha A] [--steps N] [--samples S] [--lr R]\n \ [--init W] [--guide ε] [--exact] [--seed S] [--json]\n\n\ THE FABRIC (members, parity against the CPU reference, placement):\n \ wai-quantum lab fabric [--circuits N] [--seed S] [--gateset v1|v2] [--json]\n\n\ OPTIMISATION (simulated bifurcation, annealing, exact):\n \ wai-quantum lab qio (--gset FILE | --sk N [--instance S]) [--method dsb|bsb|hdsb|hbsb|sa|exact]\n \ [--steps N] [--trials T] [--dt D] [--c1 C] [--gamma G] [--beta0 B --beta1 B] [--descend]\n \ [--target V] [--seed S] [--threads T] [--spins] [--json]\n\n\ QLDPC MEMORY (circuit level, BP-OSD):\n \ wai-quantum lab ldpc [--code 72|90|108|144|288] [--p P] [--cycles N] [--shots S] [--basis z|x|both]\n \ [--max-iter I] [--osd-order K] [--threads T] [--seed S] [--json]\n \ wai-quantum lab ldpc --circuit [--code …] [--p P] [--cycles N] [--basis z|x] (the circuit text)\n\n\ ARCHITECTURE:\n \ wai-quantum lab arch factoring [--bits 1024|1536|2048|3072|4096|6144|8192] [--json]\n \ wai-quantum lab arch storage [--target T] [--fits refit|uniform|si1000] [--json]\n \ wai-quantum lab arch cultivate [--p 5e-4|1e-3|2e-3] [--target T] [--json]\n\n\ NOISY CIRCUITS (trajectories, exact density matrix, cross-entropy):\n \ wai-quantum lab noisy --file CIRCUIT.jsonl [-n N] [--depol1 P] [--depol2 P] [--damping G] [--dephasing P]\n \ [--readout P] [--trajectories T] [--seed S] [--expect \"Z0 Z1; X3\"] [--xeb] [--exact] [--json]\n\n\ STATE VECTOR (any gate, any angle):\n \ wai-quantum lab 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 lab 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 lab list` for the algorithm catalog." ); std::process::exit(2); } fn read_stdin() -> String { let mut s = String::new(); if std::io::stdin().read_to_string(&mut s).is_err() { die("could not read stdin"); } s } /// The key a receipt is signed with ([`crate::keys`]), from `--key FILE` or /// `--demo-key`. fn signing_key(args: &[String]) -> SigningKey { let given = opt(args, "--key", "").map(std::path::PathBuf::from); let k = wai_quantum_cli::keys::load(given.as_deref(), flag(args, "--demo-key")).unwrap_or_else(|e| die(&e)); wai_quantum_cli::keys::notes(&k); k.key } /// Pull `--flag VALUE` out of the argument list. fn opt(args: &[String], long: &str, short: &str) -> Option { args.iter() .position(|a| a == long || (!short.is_empty() && a == short)) .and_then(|i| args.get(i + 1).cloned()) } fn flag(args: &[String], long: &str) -> bool { args.iter().any(|a| a == long) } fn num(args: &[String], long: &str, short: &str, dflt: T) -> T { match opt(args, long, short) { Some(v) => v.parse().unwrap_or_else(|_| die("bad numeric argument")), None => dflt, } } /// What `--gateset v1` prints on standard error, never on standard output or in /// JSON, for a circuit whose phases circuit/1 does not reproduce. const V1_NOTE: &str = "circuit/1 table: P(k≥8) deviates by up to 12868 units of 2^-30, P(k≥19) applies no phase; use --gateset v2"; /// The gateset `--gateset` names: `v2` (`wai.quantum.circuit2`, the default) or /// `v1` (the frozen `wai.quantum.circuit`). fn gateset_from(args: &[String]) -> Gateset { match opt(args, "--gateset", "").as_deref() { None | Some("v2") => Gateset::V2, Some("v1") => Gateset::V1, Some(other) => die(&format!("unknown gateset '{other}' (v1, v2)")), } } /// Build the circuit a command operates on: either a named algorithm or OpenQASM /// read from stdin (`--qasm`), in the gateset `--gateset` names. A circuit/1 /// circuit with a `P(k ≥ 8)` gets the note on standard error. fn circuit_from(args: &[String], algo: Option<&String>) -> Circuit { let gs = gateset_from(args); let mut c = if flag(args, "--qasm") { wai_quantum::qasm::dyadic_circuit(&read_stdin()).unwrap_or_else(|e| die(&format!("QASM parse: {e}"))) } else if flag(args, "--qir") { wai_quantum::quantum_qir::from_qir(&read_stdin()).unwrap_or_else(|e| die(&format!("QIR parse: {e}"))) } else { let id = algo.unwrap_or_else(|| usage()); let n: u8 = num(args, "--qubits", "-n", 3); let p: u32 = num(args, "--param", "-p", 0); qt::build_algorithm(id, n, p) .unwrap_or_else(|| die(&format!("unknown algorithm '{id}' (try `wai-quantum list`)"))) }; // Every circuit2 circuit is a valid circuit/1 circuit, so the switch is safe. c.gateset = gs; if gs == Gateset::V1 && c.ops.iter().any(|g| g.base == BaseGate::P && g.param >= 8) { eprintln!("{V1_NOTE}"); } c } fn describe(c: &Circuit) -> (String, String, u64, String) { let (backend, why) = qt::recommend_backend(c); ( format!("{backend:?}"), why.to_string(), qt::statevector_work(c), hx(&qt::circuit_id(c).unwrap_or_else(|e| die(&format!("circuit id: {e}")))), ) } /// The six reference qubits by command-line name. fn reference_qubit(name: &str) -> PhysicalQubit { match name { "us-e3" => PhysicalQubit::us_e3(), "us-e4" => PhysicalQubit::us_e4(), "ns-e3" => PhysicalQubit::ns_e3(), "ns-e4" => PhysicalQubit::ns_e4(), "meas-e4" => PhysicalQubit::meas_ns_e4(), "meas-e6" => PhysicalQubit::meas_ns_e6(), other => die(&format!("unknown qubit '{other}' (us-e3, us-e4, ns-e3, ns-e4, meas-e4, meas-e6, all)")), } } /// A duration in the unit a person would say it in. pub(crate) fn human_time(ns: f64) -> String { let s = ns / 1e9; let (m, h, d, y) = (60.0, 3600.0, 86_400.0, 365.25 * 86_400.0); if s < 1e-3 { format!("{:.1} us", ns / 1e3) } else if s < 1.0 { format!("{:.1} ms", s * 1e3) } else if s < m { format!("{s:.2} s") } else if s < h { format!("{:.1} min", s / m) } else if s < d { format!("{} h {} min", (s / h) as u64, ((s % h) / m) as u64) } else if s < y { format!("{:.1} days", s / d) } else { format!("{:.1} years", s / y) } } /// A circuit file: JSON lines, one gate per line, `{"q":[a,b],"m":[re,im, …]}` /// (rows are outputs, the first qubit most significant). Returns the gates /// and the register size they need. fn read_gate_file(path: &str) -> (Vec, u32) { use wai_quantum::gates::{complex, Gate}; let text = std::fs::read_to_string(path).unwrap_or_else(|e| die(&format!("{path}: {e}"))); let mut gates = Vec::new(); let mut n: u32 = 0; for (i, line) in text.lines().enumerate().filter(|(_, l)| !l.trim().is_empty()) { let v: serde_json::Value = serde_json::from_str(line).unwrap_or_else(|_| die(&format!("line {}: not JSON", i + 1))); let qubits: Vec = v["q"] .as_array() .unwrap_or_else(|| die(&format!("line {}: missing \"q\"", i + 1))) .iter() .map(|x| x.as_u64().unwrap_or_else(|| die("qubits are integers")) as u32) .collect(); let m: Vec = v["m"] .as_array() .unwrap_or_else(|| die(&format!("line {}: missing \"m\"", i + 1))) .iter() .map(|x| x.as_f64().unwrap_or_else(|| die("matrix entries are numbers"))) .collect(); n = n.max(qubits.iter().max().map_or(0, |q| q + 1)); gates.push(Gate { qubits, matrix: m.chunks(2).map(|c| complex(c[0], *c.get(1).unwrap_or(&0.0))).collect() }); } (gates, n) } /// A program for `estimate`: a named example, or counts from flags. /// Returns the counts and the error budget the example is stated with. fn estimate_program(args: &[String], name: Option<&String>) -> (LogicalCounts, f64) { let counted = ["--qubits", "--ccz", "--t", "--rotations", "--depth", "--measurements"] .iter() .any(|f| opt(args, f, "").is_some()); if counted { // Counts may be written as 1e9; they must still be whole numbers. let count = |f: &str| -> u64 { let v: f64 = num(args, f, "", 0.0); if !(v >= 0.0 && v.fract() == 0.0 && v < 1.8e19) { die(&format!("{f} must be a whole number")); } v as u64 }; let c = LogicalCounts { qubits: count("--qubits"), measurements: count("--measurements"), rotations: count("--rotations"), rotation_depth: count("--depth"), t_gates: count("--t"), ccz: count("--ccz"), }; if c.rotations > 0 && c.rotation_depth == 0 { die("--rotations needs --depth (non-Clifford layers containing a rotation)"); } return (c, 0.001); } match name.map(String::as_str).unwrap_or("factoring") { // 2048-bit factoring (windowed, coset representation), ε = 1/3: the // answer is checkable, so a failed run is simply repeated. "factoring" => ( LogicalCounts { qubits: 12_581, measurements: 1_080_000_000, rotations: 12, rotation_depth: 12, t_gates: 12, ccz: 3_730_000_000 }, 1.0 / 3.0, ), // Ground-state energy of a 76-orbital catalyst model by double- // factorized qubitization, ε = 0.01. "chemistry" => ( LogicalCounts { qubits: 1318, measurements: 1_370_000_000, rotations: 206_000_000, rotation_depth: 205_000_000, t_gates: 55_300_000, ccz: 135_000_000_000, }, 0.01, ), other => die(&format!("unknown program '{other}' (factoring, chemistry, or counts via --qubits ...)")), } } fn estimate_json(e: &qr::Estimate) -> String { let rounds: Vec = e .factory .iter() .flat_map(|f| f.rounds.iter()) .map(|r| { format!( "{{\"copies\":{},\"distance\":{},\"unit\":\"{}\"}}", r.copies, r.distance, match r.kind { qr::UnitKind::SpaceEfficient => "space-efficient", qr::UnitKind::ReedMuller => "reed-muller", } ) }) .collect(); format!( "{{\"code\":{},\"distance\":{},\"factories\":{},\"factory_fraction\":{},\"factory_rounds\":[{}],\ \"logical_steps\":{},\"physical_qubits\":{},\"qubit\":{},\"runtime_s\":{},\"t_states\":{},\"tiles\":{}}}", serde_json::to_string(&e.code).unwrap(), e.distance, e.factories, e.factory_fraction(), rounds.join(","), e.steps, e.physical_qubits, serde_json::to_string(&e.qubit).unwrap(), e.runtime_ns / 1e9, e.isa.t_states, e.isa.tiles, ) } /// Check a receipt of any kind this tool or the library signs: `src` is its /// JSON. Prints what it binds, and exits 1 when it does not verify. pub fn verify_receipt(src: &str) { // Either form: a receipt as sealed, or one whose joule figure // carries its acquisition class inside the signature // (`quantum_energy::Labelled`, as `wai_quantum_meter` emits). if let Some(any) = MaybeLabelled::::from_json(src) { let ok = any.verify(); let r = any.receipt(); let e = &r.estimate; println!("signature {}", if ok { "VALID" } else { "INVALID" }); println!("kind fault-tolerant resource estimate"); println!("job-hash {}", hx(&r.job_hash)); println!("qubit {} ({}, d = {})", e.qubit, e.code, e.distance); println!("physical {} qubits", e.physical_qubits); println!("runtime {}", human_time(e.runtime_ns)); println!("factories {}", e.factories); println!("joules {}", joules_line(r.joules_micro)); println!("energy class {}", energy_class_line(any.energy_class(), r.joules_micro)); println!("signer {}", r.signer_id); if !ok { std::process::exit(1); } return; } if let Some(any) = MaybeLabelled::::from_json(src) { let ok = any.verify(); let r = any.receipt(); println!("signature {}", if ok { "VALID" } else { "INVALID" }); println!("kind sparse Pauli dynamics"); println!("circuit-hash {}", hx(&r.circuit_hash)); println!("query-hash {}", hx(&r.query_hash)); println!("value {}", r.result.value); println!("bound {}", r.result.truncation_bound); println!("work {} term-updates (peak {} terms)", r.result.term_updates, r.result.peak_terms); println!("joules {}", joules_line(r.joules_micro)); println!("energy class {}", energy_class_line(any.energy_class(), r.joules_micro)); println!("signer {}", r.signer_id); if !ok { std::process::exit(1); } return; } match MaybeLabelled::::from_json(src) { None => { println!("INVALID (could not parse receipt JSON, or its energy class is malformed)"); std::process::exit(1); } Some(any) => { let ok = any.verify(); let r = any.receipt(); println!("signature {}", if ok { "VALID" } else { "INVALID" }); println!("circuit-hash {}", hx(&r.circuit_hash)); println!("state-hash {}", hx(&r.statevector_hash)); println!("qubits {}", r.n_qubits); println!("gates {}", r.n_ops); println!("work {} amplitude-updates", r.work_amp_updates); println!("joules {}", joules_line(r.joules_micro)); println!("energy class {}", energy_class_line(any.energy_class(), r.joules_micro)); println!("signer {}", r.signer_id); if !ok { std::process::exit(1); } } } } /// The research tools, with `args` the words after `wai-quantum lab`. pub fn main(args: Vec) { let argv: Vec = std::iter::once("wai-quantum lab".to_string()).chain(args).collect(); if argv.len() < 2 { usage(); } let cmd = argv[1].as_str(); let args = &argv[2..]; let algo = args.iter().find(|a| !a.starts_with('-')); match cmd { "list" => { println!("{:<14} {:<26} DESCRIPTION", "ID", "NAME"); for (id, name, desc) in qt::catalog() { println!("{id:<14} {name:<26} {desc}"); } } "qasm" => { let c = circuit_from(args, algo); print!("{}", qt::to_qasm(&c)); } "qir" => { let c = circuit_from(args, algo); print!("{}", wai_quantum::quantum_qir::to_qir(&c, !flag(args, "--no-measure"))); } "run" => { let c = circuit_from(args, algo); if let Err(e) = c.validate() { die(&format!("invalid circuit: {e:?}")); } let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}"))); let top: usize = num(args, "--top", "", 8); let (backend, why, work, cid) = describe(&c); let hash = hx(&sv.statevector_hash()); let outs = qt::top_outcomes(&sv, top); let n = c.n_qubits as usize; if flag(args, "--json") { let items: Vec = outs .iter() .map(|(k, p)| format!("{{\"state\":\"{k:0n$b}\",\"p\":{p}}}", n = n)) .collect(); println!( "{{\"circuit_id\":\"{cid}\",\"capability\":\"{}\",\"n_qubits\":{},\"n_gates\":{},\"backend\":\"{backend}\",\ \"work_amp_updates\":{work},\"statevector_hash\":\"{hash}\",\"top\":[{}]}}", c.gateset.capability(), c.n_qubits, c.ops.len(), items.join(",") ); } else { println!("circuit-id {cid}"); println!("capability {}", c.gateset.capability()); println!("qubits {}", c.n_qubits); println!("gates {}", c.ops.len()); println!("backend {backend} ({why})"); println!("work {work} amplitude-updates"); println!("state-hash {hash}"); println!("top outcomes"); for (k, p) in outs { if p > 0.0 { println!(" |{k:0n$b}> {p:.6}", n = n); } } } } "receipt" => { let c = circuit_from(args, algo); let shots: u64 = num(args, "--shots", "", 0); let seed: u64 = num(args, "--seed", "", 0); let measure = (shots > 0).then_some(Measure { seed, shots }); let key = signing_key(args); let r = QuantumReceipt::seal(&key, did_key(&key.verifying_key()), &c, measure, 0, None) .unwrap_or_else(|e| die(&format!("seal: {e:?}"))); println!("{}", r.to_json()); } "verify" => verify_receipt(&read_stdin()), // How much quantum information a state actually carries. A state vector // cannot answer this for a subsystem — half a Bell pair has no state // vector — so `--keep` traces the rest out and reports the entropy that // only a density matrix can see. // The 127-qubit heavy-hex kicked-Ising experiment, by sparse Pauli // dynamics: ⟨Z_q⟩ after `--steps` Trotter steps of RX(θ_h) on every // qubit and RZZ(−π/2) on every edge, and the bound on how far // truncation can have moved it. `--receipt` seals the run instead. "utility" => { use wai_quantum::quantum_spd::{heavy_hex_127, kicked_ising, Pauli, SpdConfig}; let theta: f64 = num(args, "--theta", "", std::f64::consts::FRAC_PI_4); let steps: u32 = num(args, "--steps", "", 20); let threshold: f64 = num(args, "--threshold", "", 1e-3); let qubit: u32 = num(args, "--qubit", "-q", 62); // The same circuit by belief-propagation tensor networks: an // independent method, whose errors are unrelated to truncating a // Pauli sum. if opt(args, "--method", "").as_deref() == Some("bptn") { if flag(args, "--receipt") { die("--receipt seals sparse Pauli dynamics runs; tensor-network runs have no receipt yet"); } let chi: usize = num(args, "--chi", "", 16); let edges = heavy_hex_127(); let mut tns = wai_quantum::quantum_bptn::Tns::new(127, &edges, chi).unwrap_or_else(|e| die(&e.to_string())); tns.run(&kicked_ising(127, &edges, theta, steps)).unwrap_or_else(|e| die(&e.to_string())); let value = tns.expect1(qubit, Pauli::Z).unwrap_or_else(|e| die(&e.to_string())); if flag(args, "--json") { println!( "{{\"chi\":{chi},\"discarded\":{},\"max_bond\":{},\"method\":\"bptn\",\"qubit\":{qubit},\ \"steps\":{steps},\"theta_h\":{theta},\"value\":{value}}}", tns.discarded(), tns.max_bond() ); } else { println!(" {value:.6}"); println!("discarded {:.3e} (squared singular-value weight truncated)", tns.discarded()); println!("max bond {}", tns.max_bond()); } return; } if opt(args, "--method", "").as_deref().is_some_and(|m| m != "spd") { die("--method is spd (the default) or bptn"); } let job = SpdJob { circuit: kicked_ising(127, &heavy_hex_127(), theta, steps), observable: vec![(1.0, vec![(qubit, Pauli::Z)])], ones: vec![], config: SpdConfig { threshold, max_weight: None }, }; if flag(args, "--receipt") { let key = signing_key(args); let r = SpdReceipt::seal(&key, did_key(&key.verifying_key()), &job, None, |run| { run(); 0 }) .unwrap_or_else(|e| die(&format!("seal: {e}"))); println!("{}", r.to_json()); return; } let r = job.run().unwrap_or_else(|e| die(&e.to_string())); if flag(args, "--json") { println!( "{{\"bound\":{},\"final_terms\":{},\"peak_terms\":{},\"qubit\":{qubit},\"steps\":{steps},\ \"term_updates\":{},\"theta_h\":{theta},\"threshold\":{threshold},\"value\":{}}}", r.truncation_bound, r.final_terms, r.peak_terms, r.term_updates, r.value ); } else { println!(" {:.6}", r.value); println!("bound {:.3e} (|exact - value| is at most this)", r.truncation_bound); println!("peak terms {}", r.peak_terms); println!("work {} term-updates", r.term_updates); } } // What a fault-tolerant machine would need to run a program: code // distance, T factories, physical qubits and wall-clock time, on the // six reference qubits. `--receipt` seals one estimate, re-checkable // bit for bit by anyone who re-runs it. "estimate" => { let (counts, default_budget) = estimate_program(args, algo); let budget: f64 = num(args, "--budget", "", default_budget); let slowdown: u64 = num(args, "--slowdown", "", 1); let cfg = EstimateConfig { error_budget: budget, slowdown, ..EstimateConfig::new(budget) }; let code = opt(args, "--code", "").map(|c| match c.as_str() { "surface" => QecModel::surface_gate(), "surface-meas" => QecModel::surface_measurement(), "floquet" => QecModel::floquet(), other => die(&format!("unknown code '{other}' (surface, surface-meas, floquet)")), }); // A code model with given constants, e.g. as `qec --fit` measured them. let a = opt(args, "--prefactor", "").map(|v| v.parse::().unwrap_or_else(|_| die("bad --prefactor"))); let t = opt(args, "--threshold", "").map(|v| v.parse::().unwrap_or_else(|_| die("bad --threshold"))); let code = if a.is_some() || t.is_some() { let base = code.unwrap_or_else(QecModel::surface_gate); Some(QecModel { name: format!("{} (given)", base.name), prefactor: a.unwrap_or(base.prefactor), threshold: t.unwrap_or(base.threshold), ..base }) } else { code }; let which = opt(args, "--qubit", "").unwrap_or_else(|| { if flag(args, "--receipt") || flag(args, "--frontier") { "ns-e4".into() } else { "all".into() } }); let qubits: Vec = if which == "all" { PhysicalQubit::presets() } else { vec![reference_qubit(&which)] }; if flag(args, "--receipt") { if qubits.len() != 1 { die("--receipt seals one estimate: pick a --qubit"); } let job = ResourceJob { program: Program::Logical(counts), qubit: qubits[0].clone(), code, config: cfg }; let key = signing_key(args); let r = ResourceReceipt::seal(&key, did_key(&key.verifying_key()), &job, None, |run| { run(); 0 }) .unwrap_or_else(|e| die(&format!("estimate: {e:?}"))); println!("{}", r.to_json()); return; } let run = |q: &PhysicalQubit, cfg: &EstimateConfig| match &code { Some(c) => qr::estimate(&counts, q, c, cfg), None => qr::estimate_best(&counts, q, cfg), }; if flag(args, "--frontier") { // Slow the program down to need fewer factories: each point // uses fewer qubits than every faster one. let q = &qubits[0]; let c = code.clone().unwrap_or_else(|| { run(q, &cfg).map(|e| { QecModel::for_qubit(q).into_iter().find(|m| m.name == e.code).expect("code of the estimate") }) .unwrap_or_else(|e| die(&format!("estimate: {e:?}"))) }); let isa = qr::layout(&counts, budget); let json = flag(args, "--json"); let points = qr::frontier(&isa, q, &c, &cfg, 1024); if json { let rows: Vec = points.iter().map(estimate_json).collect(); println!("[{}]", rows.join(",")); return; } println!("{:>8} {:>3} {:>9} {:>14} runtime", "slowdown", "d", "factories", "physical"); for e in &points { println!( "{:>8} {:>3} {:>9} {:>14} {}", e.steps / isa.min_steps, e.distance, e.factories, e.physical_qubits, human_time(e.runtime_ns) ); } return; } let estimates: Vec<_> = qubits.iter().map(|q| (q, run(q, &cfg))).collect(); if flag(args, "--json") { let rows: Vec = estimates .iter() .map(|(q, e)| match e { Ok(e) => estimate_json(e), Err(err) => format!("{{\"error\":\"{err:?}\",\"qubit\":{}}}", serde_json::to_string(&q.name).unwrap()), }) .collect(); println!("[{}]", rows.join(",")); return; } let isa = qr::layout(&counts, budget); println!("program {} tiles, {} logical steps, {} T states (budget {budget})", isa.tiles, isa.min_steps, isa.t_states); println!( "{:<12} {:<22} {:>3} {:>9} {:>8} {:>14} runtime", "qubit", "code", "d", "factories", "share", "physical" ); for (q, e) in &estimates { match e { Ok(e) => println!( "{:<12} {:<22} {:>3} {:>9} {:>7.1}% {:>14} {}", e.qubit, e.code, e.distance, e.factories, 100.0 * e.factory_fraction(), e.physical_qubits, human_time(e.runtime_ns) ), Err(err) => println!("{:<12} {err:?}", q.name), } } } // Error correction at circuit level: sample a noisy Clifford circuit, // derive its detector error model, decode it exactly. "qec" => { use wai_quantum::quantum_frame::{error_model, parse, sample, surface_code_memory}; use wai_quantum::quantum_match::{fit_code_model, memory_experiment, MatchingGraph}; let shots: usize = num(args, "--shots", "", 100_000); let seed: u64 = num(args, "--seed", "", 1); let json = flag(args, "--json"); if flag(args, "--fit") { // The surface code's logical error per round on a grid, and // the code model a resource estimate would otherwise assume. let mut pts = Vec::new(); if !json { println!("{:>3} {:>7} {:>9} {:>9} {:>11}", "d", "p", "shots", "failures", "per round"); } for d in [3u32, 5, 7] { for p in [0.002, 0.003, 0.004, 0.005] { let n = if d == 7 { shots / 4 } else { shots }; let t = memory_experiment(d, d, p, n, seed).unwrap_or_else(|e| die(&format!("{e:?}"))); if !json { println!("{:>3} {:>7} {:>9} {:>9} {:>11.3e}", d, p, t.shots, t.failures, t.per_round()); } pts.push(t); } } let f = fit_code_model(&pts).unwrap_or_else(|| die("too few failures to fit; raise --shots")); if json { println!("{{\"prefactor\":{},\"threshold\":{}}}", f.prefactor, f.threshold); } else { println!("fit P(d) = a·(p/p*)^((d+1)/2) per round, a = {:.4}, p* = {:.5}", f.prefactor, f.threshold); println!("estimate wai-quantum estimate factoring --prefactor {:.4} --threshold {:.5}", f.prefactor, f.threshold); } return; } let d: u32 = num(args, "--distance", "-d", 3); let rounds: u32 = num(args, "--rounds", "", d); let p: f64 = num(args, "--p", "", 0.001); let text = match opt(args, "--file", "") { Some(path) => std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))), None => { if d < 2 || rounds < 1 || !(0.0..=0.5).contains(&p) { die("need -d ≥ 2, --rounds ≥ 1 and 0 ≤ --p ≤ 0.5"); } surface_code_memory(d, rounds, p) } }; if flag(args, "--circuit") { print!("{text}"); return; } let c = parse(&text).unwrap_or_else(|e| die(&format!("circuit: {e}"))); let model = error_model(&c).unwrap_or_else(|e| die(&format!("model: {e}"))); if flag(args, "--dem") { print!("{}", model.to_text()); return; } let g = MatchingGraph::from_model(&model).unwrap_or_else(|e| die(&format!("matching graph: {e:?}"))); let s = sample(&c, shots, seed); let decoder = opt(args, "--decoder", "").unwrap_or_else(|| "exact".into()); let failures = match decoder.as_str() { "exact" => g.failures(&s), "uf" => wai_quantum::quantum_uf::UnionFind::new(&g).failures(&s), _ => die("--decoder: exact or uf"), } .unwrap_or_else(|e| die(&format!("decode: {e:?}"))); let rate = failures as f64 / shots as f64; if json { println!( "{{\"decoder\":\"{decoder}\",\"detectors\":{},\"failures\":{failures},\"mechanisms\":{},\"observables\":{},\"qubits\":{},\"seed\":{seed},\"shots\":{shots}}}", model.detectors, model.mechanisms.len(), model.observables, c.qubits ); } else { println!("circuit {} qubits, {} detectors, {} observables", c.qubits, model.detectors, model.observables); println!("model {} error mechanisms", model.mechanisms.len()); println!("failures {failures} of {shots} shots ({rate:.3e} per shot, {decoder} decoder)"); if opt(args, "--file", "").is_none() { let per_round = wai_quantum::quantum_match::MemoryPoint { distance: d, rounds, p, shots: shots as u64, failures }.per_round(); println!("per round {per_round:.3e} (d = {d}, {rounds} rounds, p = {p})"); } } } // A circuit's amplitudes by tensor-network contraction: a seeded // hyper-optimised contraction tree, sliced to a memory cap, contracted // in a fixed order. Circuits are JSON lines, one gate per line: // {"q":[0,1],"m":[re,im, …]} (rows are outputs, first qubit most // significant). "tn" => { use wai_quantum::quantum_tn::{contract_plan, find_path, Families, Network, PathOptions}; let path = opt(args, "--file", "").unwrap_or_else(|| die("--file CIRCUIT.jsonl is required")); let (gates, n) = read_gate_file(&path); let n: u32 = num(args, "--qubits", "-n", n); let bits: Vec = match opt(args, "--bits", "") { Some(b) => b.chars().map(|c| if c == '1' { 1 } else if c == '0' { 0 } else { die("--bits is a 0/1 string") }).collect(), None => vec![0; n as usize], }; let open: Vec = opt(args, "--open", "") .map(|l| l.split(',').map(|q| q.trim().parse().unwrap_or_else(|_| die("bad --open qubit"))).collect()) .unwrap_or_default(); if open.len() > 20 { die("at most 20 open qubits (2^20 amplitudes)"); } // 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(); // Projections fixed through diagonal gates, then absorption. if whole { net.simplify(); } else { net.reduce(); } let shape = net.shape(); let cap: Option = opt(args, "--cap-log2", "").map(|v| 2f64.powf(v.parse().unwrap_or_else(|_| die("bad --cap-log2")))); let options = PathOptions { trials: num(args, "--trials", "", 64), seed: num(args, "--seed", "", 0), max_size: cap, reconfigure: 8, families: Families::Both, }; let plan = find_path(&shape, &net.dims, &net.open, &options); let flops = plan.total_flops(); let plan_line = if flops > 0.0 { format!("10^{:.3} operations, largest intermediate 2^{:.1}, {} slices", flops.log10(), plan.per_slice.max_size.log2(), plan.slices) } else { "no contractions needed".to_string() }; if flag(args, "--plan-only") || flops > 1e13 { if !flag(args, "--plan-only") { eprintln!("wai-quantum: the plan costs 10^{:.1} operations; printing the plan only (pass a smaller circuit, or --plan-only)", flops.log10()); } if flag(args, "--json") { println!( "{{\"flops\":{flops},\"log2_max_size\":{},\"slices\":{},\"tensors\":{}}}", plan.per_slice.max_size.log2(), plan.slices, shape.len() ); } else { println!("network {raw} tensors, {} after simplification", shape.len()); println!("plan {plan_line}"); } return; } let t = contract_plan(&net, &plan).unwrap_or_else(|e| die(&format!("{e:?}"))); if flag(args, "--json") { let amps: Vec = t.data.iter().map(|c| format!("[{},{}]", c.re, c.im)).collect(); println!( "{{\"amplitudes\":[{}],\"flops\":{flops},\"log2_max_size\":{},\"slices\":{}}}", amps.join(","), plan.per_slice.max_size.log2(), plan.slices ); } else { println!("network {raw} tensors, {} after simplification", shape.len()); println!("plan {plan_line}"); for (i, c) in t.data.iter().enumerate() { println!("amp[{i}] {:+.12e} {:+.12e}i", c.re, c.im); } } } // Gaussian boson sampling, exactly: squeezed light through an // interferometer with loss, and the click statistics it must show. "gbs" => { use wai_quantum::quantum_gbs::{random_unitary, GaussianState, Mat, C}; let r: Vec = opt(args, "--squeezing", "") .unwrap_or_else(|| die("--squeezing r1,r2,… is required (one value per mode)")) .split(',') .map(|x| x.trim().parse().unwrap_or_else(|_| die("bad --squeezing value"))) .collect(); let n = r.len(); let eta: f64 = num(args, "--eta", "", 1.0); let u = match opt(args, "--unitary", "") { Some(path) => { // JSON: an array of N rows, each N [re, im] pairs. let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))); let v: serde_json::Value = serde_json::from_str(&text).unwrap_or_else(|_| die("--unitary is not JSON")); let rows = v.as_array().unwrap_or_else(|| die("--unitary: an array of rows")); if rows.len() != n { die("--unitary needs one row per mode"); } let mut m = Mat::zeros(n); for (i, row) in rows.iter().enumerate() { let row = row.as_array().filter(|r| r.len() == n).unwrap_or_else(|| die("--unitary rows need one entry per mode")); for (j, z) in row.iter().enumerate() { let re = z[0].as_f64().unwrap_or_else(|| die("--unitary entries are [re, im]")); let im = z[1].as_f64().unwrap_or(0.0); m.a[i * n + j] = C { re, im }; } } m } None => random_unitary(n, num(args, "--seed", "", 0)), }; let st = GaussianState::squeezed(&r, &u, eta).unwrap_or_else(|e| die(&format!("{e:?}"))); let json = flag(args, "--json"); if let Some(path) = opt(args, "--check", "") { // Lines of 0/1, one click pattern per line. let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))); let samples: Vec> = text .lines() .filter(|l| !l.trim().is_empty()) .map(|l| l.trim().chars().map(|c| c == '1').collect()) .collect(); let c = st.check_samples(&samples).unwrap_or_else(|e| die(&format!("{e:?} (every line needs {n} characters)"))); if json { println!( "{{\"mean_clicks_exact\":{},\"mean_clicks_observed\":{},\"shots\":{},\"worst_pair_z\":{},\"worst_single_z\":{}}}", c.mean_clicks.1, c.mean_clicks.0, c.shots, c.worst_pair, c.worst_single ); } else { println!("shots {}", c.shots); println!("mean clicks {:.4} observed, {:.4} exact", c.mean_clicks.0, c.mean_clicks.1); println!("worst |z| {:.2} over single modes, {:.2} over pairs", c.worst_single, c.worst_pair); let verdict = if c.worst_single < 5.0 && c.worst_pair < 5.0 { "consistent" } else { "INCONSISTENT" }; println!("verdict {verdict} with the state's exact one- and two-mode click statistics"); } return; } if let Some(k) = opt(args, "--sample", "") { let k: usize = k.parse().unwrap_or_else(|_| die("bad --sample count")); for s in st.sample_clicks(k, num(args, "--seed", "", 0) ^ 0x9e37) { println!("{}", s.iter().map(|&b| if b { '1' } else { '0' }).collect::()); } return; } let (p1, _) = st.click_marginals(); let none = st.click_probability(&vec![false; n]).unwrap_or(f64::NAN); let all = if n <= 22 { st.click_probability(&vec![true; n]).ok() } else { None }; if json { let p1s: Vec = p1.iter().map(|p| p.to_string()).collect(); println!( "{{\"click_marginals\":[{}],\"modes\":{n},\"p_all_click\":{},\"p_no_click\":{none}}}", p1s.join(","), all.map_or("null".to_string(), |p| p.to_string()) ); } else { println!("modes {n}, transmission {eta}"); println!("mean clicks {:.4}", p1.iter().sum::()); println!("P(no click) {none:.6e}"); if let Some(p) = all { println!("P(all click) {p:.6e}"); } for (i, p) in p1.iter().enumerate() { println!("mode {i:>3} P(click) = {p:.6}"); } } } // Ground states of spin chains by DMRG, against the exact energy // where one is known (the Ising chain at any length, any chain by // exact diagonalisation up to 22 sites). "dmrg" => { use wai_quantum::quantum_dmrg::*; let n: usize = num(args, "--sites", "-n", 40); if n < 2 { die("-n: at least 2 sites"); } let j: f64 = num(args, "--j", "", 1.0); let h: f64 = num(args, "--h", "", if opt(args, "--chain", "").as_deref() == Some("heisenberg") { 0.0 } else { 1.0 }); let kind = opt(args, "--chain", "").unwrap_or_else(|| "ising".into()); let chain = match kind.as_str() { "ising" => Chain::ising(n, j, h), "heisenberg" => Chain::heisenberg(n, j, num(args, "--jz", "", j), h), _ => die("--chain: ising or heisenberg"), }; let cfg = DmrgConfig { max_bond: num(args, "--bond", "", 64), sweeps: num(args, "--sweeps", "", 5), cutoff: num(args, "--cutoff", "", 1e-12), ..DmrgConfig::default() }; let r = dmrg(&chain, &cfg); let exact = if kind == "ising" { Some(ising_exact_energy(n, j, h)) } else if flag(args, "--exact") || n <= 16 { exact_ground_energy(&chain) } else { None }; let max_bond = r.mps.dims.iter().copied().max().unwrap_or(1); if flag(args, "--json") { let ent: Vec = r.entropies.iter().map(|e| e.to_string()).collect(); let sw: Vec = r.sweep_energies.iter().map(|e| e.to_string()).collect(); println!( "{{\"bond\":{max_bond},\"chain\":\"{kind}\",\"discarded\":{},\"energy\":{},\"entropies\":[{}],\"exact\":{},\"sites\":{n},\"sweeps\":[{}]}}", r.discarded, r.energy, ent.join(","), exact.map_or("null".into(), |e| e.to_string()), sw.join(",") ); } else { println!("{kind} chain, {n} sites, bond ≤ {}, {} sweeps", cfg.max_bond, cfg.sweeps); println!("energy {:.12} ({:.12} per site)", r.energy, r.energy / n as f64); if let Some(e) = exact { println!("exact {e:.12} (error per site {:.1e})", (r.energy - e) / n as f64); } println!("discarded {:.1e} largest bond {max_bond}", r.discarded); println!("entropy, mid {:.6}", r.entropies[n / 2 - 1]); } } // Time evolution of a spin chain by TDVP from a product state. The // Ising chain started with every σˣ = +1 is held at each report to // its free-fermion solution. "tdvp" => { use wai_quantum::quantum_dmrg::{Chain, Mps, Op, Scalar, C}; use wai_quantum::quantum_tdvp::*; let n: usize = num(args, "--sites", "-n", 40); if n < 2 { die("-n: at least 2 sites"); } let kind = opt(args, "--chain", "").unwrap_or_else(|| "ising".into()); let j: f64 = num(args, "--j", "", 1.0); let h: f64 = num(args, "--h", "", if kind == "heisenberg" { 0.0 } else { 1.0 }); let chain = match kind.as_str() { "ising" => Chain::ising(n, j, h), "heisenberg" => Chain::heisenberg(n, j, num(args, "--jz", "", j), h), _ => die("--chain: ising or heisenberg"), }; let start_kind = opt(args, "--start", "").unwrap_or_else(|| if kind == "ising" { "plus".into() } else { "neel".into() }); let start = match start_kind.as_str() { "plus" => Mps::product(&vec![[C::ONE, C::ONE]; n]), "neel" => Mps::basis(&(0..n).map(|i| i % 2 == 1).collect::>()), "up" => Mps::basis(&vec![false; n]), _ => die("--start: plus, neel or up"), }; let dt: f64 = num(args, "--dt", "", 0.05); let t_end: f64 = num(args, "--t", "", 2.0); if dt.is_nan() || dt <= 0.0 || t_end.is_nan() || t_end < 0.0 { die("--dt must be positive and --t not negative"); } let cfg = TdvpConfig { dt, max_bond: num(args, "--bond", "", 32), cutoff: num(args, "--cutoff", "", 1e-20), ..TdvpConfig::default() }; let steps = (t_end / dt).round() as usize; let every: usize = num(args, "--every", "", ((0.5 / dt).round() as usize).max(1)); // The observable at the middle site: σˣ for the Ising chain, σᶻ // for the Heisenberg chain. let (obs, obs_name) = if kind == "ising" { (Op::x(), "x") } else { (Op::z(), "z") }; let referee = kind == "ising" && start_kind == "plus"; let mut ev = Evolution::new(&chain, &start, cfg); let e0 = ev.energy(); let mut discarded = 0.0; let mut rows: Vec = Vec::new(); let json = flag(args, "--json"); if !json { println!("{kind} chain, {n} sites, start {start_kind}, bond ≤ {}, dt {dt}", cfg.max_bond); print!(" t {obs_name}(mid) "); if referee { print!("exact max error rate exact rate "); } println!("S(mid) bond |ΔE| discarded"); } for s in 1..=steps { discarded += ev.step().discarded; if s % every != 0 && s != steps { continue; } let t = ev.time(); let o = ev.expect(obs); let s_mid = ev.entropies()[n / 2 - 1]; let de = (ev.energy() - e0).abs(); let rate = -ev.overlap(&start).norm2().ln() / n as f64; let ff = referee.then(|| ising_quench(n, j, h, t)); let max_err = ff.as_ref().map(|f| o.iter().zip(&f.x).map(|(p, q)| (p - q).abs()).fold(0.0, f64::max)); if json { rows.push(format!( "{{\"bond\":{},\"discarded\":{discarded},\"energy_drift\":{de},\"entropy\":{s_mid},\"exact\":{},\"exact_rate\":{},\"max_error\":{},\"mid\":{},\"rate\":{rate},\"t\":{t}}}", ev.max_bond(), ff.as_ref().map_or("null".into(), |f| f.x[n / 2].to_string()), ff.as_ref().map_or("null".into(), |f| (-f.echo.ln() / n as f64).to_string()), max_err.map_or("null".into(), |e| e.to_string()), o[n / 2] )); } else { print!("{t:6.2} {:+.10} ", o[n / 2]); if let (Some(f), Some(e)) = (&ff, max_err) { print!("{:+.10} {e:.2e} {rate:.8} {:.8} ", f.x[n / 2], -f.echo.ln() / n as f64); } println!("{s_mid:.4} {:4} {de:.1e} {discarded:.1e}", ev.max_bond()); } } if json { println!("{{\"chain\":\"{kind}\",\"dt\":{dt},\"rows\":[{}],\"sites\":{n},\"start\":\"{start_kind}\"}}", rows.join(",")); } } // A molecule's ground state by DMRG from its integrals, against // Hartree–Fock and full configuration interaction where it fits. "chem" => { use wai_quantum::quantum_chem::*; let path = opt(args, "--fcidump", "").unwrap_or_else(|| die("--fcidump FILE is required")); let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))); let fd = read_fcidump(&text).unwrap_or_else(|e| die(&format!("{path}: {e}"))); let cfg = ChemConfig { max_bond: num(args, "--bond", "", 64), sweeps: num(args, "--sweeps", "", 4), penalty: num(args, "--penalty", "", 1.0), ..ChemConfig::default() }; let mpo_bond = molecular_mpo(&fd, &cfg).max_bond(); let hf = hartree_fock_energy(&fd); let r = ground_state(&fd, &cfg); let full = if flag(args, "--no-fci") { None } else { fci(&fd) }; if flag(args, "--json") { let sw: Vec = r.sweep_energies.iter().map(|e| e.to_string()).collect(); println!( "{{\"bond\":{},\"discarded\":{},\"electrons\":{},\"energy\":{},\"fci\":{},\"fci_determinants\":{},\"hartree_fock\":{hf},\"mpo_bond\":{mpo_bond},\"orbitals\":{},\"sweeps\":[{}]}}", cfg.max_bond, r.discarded, fd.nelec, r.energy, full.as_ref().map_or("null".into(), |f| f.energy.to_string()), full.as_ref().map_or("null".into(), |f| f.dim.to_string()), fd.norb, sw.join(",") ); } else { println!("{path}: {} orbitals, {} electrons, MPO bond {mpo_bond}, DMRG bond ≤ {}", fd.norb, fd.nelec, cfg.max_bond); println!("Hartree–Fock {hf:.10}"); for (i, e) in r.sweep_energies.iter().enumerate() { println!("sweep {:>2} {e:.10}", i + 1); } println!("discarded {:.1e}", r.discarded); match &full { Some(f) => println!("full CI {:.10} ({} determinants; DMRG − full CI {:.2e})", f.energy, f.dim, r.energy - f.energy), None => println!("full CI not run"), } } } // A neural quantum state for a spin model on a chain or an open // square lattice, against exact diagonalisation where it fits. "nqs" => { use wai_quantum::quantum_nqs::*; let lx: usize = num(args, "--lx", "", 10); let ly: usize = num(args, "--ly", "", 1); let n = lx * ly; if !(2..=400).contains(&n) { die("--lx × --ly: 2 to 400 spins"); } let periodic = flag(args, "--periodic"); let kind = opt(args, "--model", "").unwrap_or_else(|| "heisenberg".into()); let nn = if ly == 1 { chain_edges(n, periodic) } else { square_edges(lx, ly, periodic) }; let nnn = if ly == 1 { chain_next_nearest(n, periodic) } else { square_diagonals(lx, ly, periodic) }; let model = match kind.as_str() { "heisenberg" => SpinModel::heisenberg(n, &nn, 1.0), "j1j2" => SpinModel::j1j2(n, &nn, &nnn, 1.0, num(args, "--j2", "", 0.5)), "ising" => SpinModel::ising(n, &nn, 1.0, num(args, "--h", "", 1.0)), _ => die("--model: heisenberg, j1j2 or ising"), }; let sector = if model.conserves_magnetisation() { Some(0) } else { None }; let exact = flag(args, "--exact"); if exact && n > 20 { die("--exact enumerates every configuration: at most 20 spins"); } let complex = flag(args, "--complex"); // Marshall's sign for the real network on the unfrustrated model; // a complex network finds its own. let marshall = kind == "heisenberg" && !complex; let mut cfg = VmcConfig { samples: num(args, "--samples", "", 1000), exact, magnetisation: sector, lr: num(args, "--lr", "", 0.02), seed: num(args, "--seed", "", 7), ..VmcConfig::default() }; let steps: usize = num(args, "--steps", "", 300); let alpha: usize = num(args, "--alpha", "", 2); let seed: u64 = num(args, "--seed", "", 7); let reference = if n <= 20 { model.exact_ground_energy(sector) } else { None }; // Projection onto translations at momentum 0 or π, on a square // lattice optionally its rotations and reflections too, and a // parity under the global flip; the lattice must be periodic. let point = opt(args, "--point-group", "").map(|g| match g.as_str() { "a1" => PointIrrep::A1, "a2" => PointIrrep::A2, "b1" => PointIrrep::B1, "b2" => PointIrrep::B2, _ => die("--point-group: a1, a2, b1 or b2"), }); let reflections = opt(args, "--reflections", "").map(|g| match g.as_str() { "a1" => RectIrrep::A1, "a2" => RectIrrep::A2, "b1" => RectIrrep::B1, "b2" => RectIrrep::B2, _ => die("--reflections: a1, a2, b1 or b2"), }); let symmetry = flag(args, "--project").then(|| { if !periodic { die("--project needs --periodic"); } let kpi = flag(args, "--kpi"); if kpi && (lx % 2 == 1 || (ly > 1 && ly % 2 == 1)) { die("--kpi needs even lengths"); } let base = match (point, reflections) { (Some(_), Some(_)) => die("--point-group or --reflections, not both"), (Some(_), None) if lx != ly || lx < 3 => die("--point-group needs a square lattice, --lx = --ly ≥ 3"), (Some(irrep), None) => Symmetry::square(lx, kpi, irrep), (None, Some(_)) if lx < 3 || ly < 3 => die("--reflections needs both lengths ≥ 3"), (None, Some(irrep)) => Symmetry::rectangle(lx, ly, kpi, kpi, irrep), (None, None) if ly == 1 => Symmetry::ring(n, kpi), (None, None) => Symmetry::torus(lx, ly, kpi, kpi), }; base.with_spin_flip(flag(args, "--odd")) }); if (point.is_some() || reflections.is_some()) && symmetry.is_none() { die("--point-group and --reflections need --project"); } // Parameters start uniform in ±init. A projected network that // starts near zero is nearly the projection of a constant, the // ferromagnet: an eigenstate, which reconfiguration barely moves. let init: f64 = num(args, "--init", "", if symmetry.is_some() { 0.1 } else { 0.01 }); // A projected network's amplitude is a sum that cancels near its // nodes, where the local energy is large and |ψ|² rarely // samples; a guide floor visits there. cfg.guide = num(args, "--guide", "", if symmetry.is_some() { 0.1 } else { 0.0 }); // The exact ground energy in the projection's own sector. let in_sector = match &symmetry { Some(sym) if n <= 24 => model.exact_ground_energy_in(sector, sym), _ => None, }; fn run(model: SpinModel, rbm: Rbm, cfg: VmcConfig, steps: usize) -> (usize, Estimate, Estimate) { let params = rbm.params(); let mut vmc = Vmc::new(model, rbm, cfg); let mut last = vmc.estimate(); for _ in 0..steps { last = vmc.step(); } (params, last, vmc.estimate()) } let (params, last, fin) = if complex { let mut rbm = Rbm::complex(n, alpha, init, seed); if marshall { rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); } if let Some(sym) = symmetry { rbm = rbm.with_symmetry(sym); } run(model, rbm, cfg, steps) } else { let mut rbm = Rbm::new(n, alpha, init, seed); if marshall { rbm = rbm.with_marshall_sign(checkerboard(lx, ly)); } if let Some(sym) = symmetry { rbm = rbm.with_symmetry(sym); } run(model, rbm, cfg, steps) }; if flag(args, "--json") { println!( "{{\"complex\":{complex},\"energy\":{},\"error\":{},\"exact\":{},\"exact_in_sector\":{},\"last_step\":{},\"model\":\"{kind}\",\"parameters\":{params},\"sites\":{n},\"steps\":{steps},\"variance\":{}}}", fin.energy, fin.error, reference.map_or("null".into(), |e| e.to_string()), in_sector.map_or("null".into(), |e| e.to_string()), last.energy, fin.variance ); } else { println!("{kind}, {lx}×{ly}{}, {} RBM with {params} parameters, {steps} steps{}", if periodic { " periodic" } else { "" }, if complex { "complex" } else { "real" }, if exact { ", exact expectations" } else { "" }); println!("energy {:.8} ± {:.8} variance {:.2e}", fin.energy, fin.error, fin.variance); if let Some(e) = reference { println!("exact {e:.8} relative error {:.2e}", (fin.energy - e) / e.abs()); } if let Some(e) = in_sector { println!("in its sector {e:.8} relative error {:.2e}", (fin.energy - e) / e.abs()); } } } // The fabric this process stands on: its detected members, and the // fixed-point core's 32-bit-lane port certified against the native // simulator on random circuits of the gateset `--gateset` names, then // placed. The member lines and the certificate's member describe the // host, so they differ between platforms; the parity and the corpus // digest do not. "fabric" => { use wai_quantum::quantum_fabric::*; let circuits: usize = num(args, "--circuits", "", 200); let seed: u64 = num(args, "--seed", "", 7); let gs = gateset_from(args); let kernel = apply_kernel(gs); let mut member = detect_local(); let cert = certify_lanes(gs, &member.id, circuits, seed); member.kernels.push(KernelEntry { kernel: kernel.into(), implementation: "lanes32".into(), certificate: cert, energy: None }); let placed = place(std::slice::from_ref(&member), kernel, &Parity::BitExact); let parity = |p: &Parity| match p { Parity::BitExact => "bit-exact".to_string(), Parity::Within { max_ulp, max_abs } => format!("within {max_ulp} ulp ({max_abs:.3e})"), Parity::Statistical { fidelity, shots } => format!("statistical, fidelity {fidelity:.6} over {shots} shots"), Parity::Uncertified => "uncertified".to_string(), }; let hex = |b: &[u8]| b.iter().map(|x| format!("{x:02x}")).collect::(); if flag(args, "--json") { let kernels: Vec = member .kernels .iter() .map(|k| format!("{{\"implementation\":{:?},\"kernel\":{:?},\"parity\":{:?},\"generator\":{:?},\"cases\":{},\"corpus\":\"{}\"}}", k.implementation, k.kernel, parity(&k.certificate.parity), k.certificate.generator, k.certificate.cases, hex(&k.certificate.corpus))) .collect(); println!( "{{\"capability\":\"{}\",\"member\":{:?},\"class\":\"{}\",\"presence\":\"{}\",\"detail\":{:?},\"kernels\":[{}],\"placed\":{}}}", gs.capability(), member.id, member.class.label(), member.presence.label(), member.detail, kernels.join(","), placed.map_or("null".to_string(), |(m, k)| format!("{:?}", format!("{} {}", m.id, k.implementation))) ); } else { println!("member {} ({}, {})", member.id, member.class.label(), member.presence.label()); println!("detail {}", member.detail); println!("kernel {kernel}"); for k in member.kernels.iter().filter(|k| k.kernel == kernel) { if k.implementation == REFERENCE { println!(" {:<11}{} (the reference)", k.implementation, parity(&k.certificate.parity)); } else { println!(" {:<11}{} over {} cases", k.implementation, parity(&k.certificate.parity), k.certificate.cases); println!(" {:<11}corpus {} from {}", "", &hex(&k.certificate.corpus)[..16], k.certificate.generator); } } match placed { Some((m, k)) => println!("placed {} · {} (bit-exact required)", m.id, k.implementation), None => println!("placed nothing meets bit-exact"), } } } // Quantum-inspired optimisation: Max-Cut (G-set text) or a spin // glass, by simulated bifurcation, annealing, or exact enumeration. // Max-Cut reports cuts; a spin glass reports energies. "qio" => { use wai_quantum::quantum_qio::*; let (problem, graph) = if let Some(path) = opt(args, "--gset", "") { let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}"))); let g = read_gset(&text).unwrap_or_else(|e| die(&e.to_string())); (g.ising(), Some(g)) } else if let Some(n) = opt(args, "--sk", "") { let n: usize = n.parse().unwrap_or_else(|_| die("bad --sk size")); (Ising::sk(n, num(args, "--instance", "", 0)), None) } else { die("qio needs --gset FILE or --sk N") }; let score = |s: &[i8]| graph.as_ref().map_or_else(|| problem.energy(s), |g| g.cut(s)); let maximise = graph.is_some(); let method = opt(args, "--method", "").unwrap_or_else(|| "dsb".into()); let threads: usize = num(args, "--threads", "", std::thread::available_parallelism().map_or(1, |n| n.get())); let seed: u64 = num(args, "--seed", "", 1); let trials: usize = num(args, "--trials", "", 64); let run = match method.as_str() { "exact" => Run { trials: vec![ground_state(&problem, threads).unwrap_or_else(|e| die(&e.to_string()))] }, "sa" => { let cfg = SaConfig { sweeps: num(args, "--steps", "", 10_000), beta0: num(args, "--beta0", "", 0.1), beta1: num(args, "--beta1", "", 3.0), trials, seed }; simulated_annealing(&problem, &cfg, threads) } m => { let v = match m { "bsb" => Variant::Ballistic, "dsb" => Variant::Discrete, "hbsb" => Variant::HeatedBallistic(0.0), "hdsb" => Variant::HeatedDiscrete(0.0), _ => die("--method: dsb, bsb, hdsb, hbsb, sa or exact"), }; let mut cfg = SbConfig::published(v, num(args, "--steps", "", 10_000), trials, seed); cfg.dt = num(args, "--dt", "", cfg.dt); cfg.c1 = num(args, "--c1", "", cfg.c1); if let Some(g) = opt(args, "--gamma", "") { let g: f64 = g.parse().unwrap_or_else(|_| die("bad --gamma")); cfg.variant = match cfg.variant { Variant::HeatedBallistic(_) => Variant::HeatedBallistic(g), Variant::HeatedDiscrete(_) => Variant::HeatedDiscrete(g), other => other, }; } simulated_bifurcation(&problem, &cfg, threads) } }; let finals: Vec> = run.trials.iter().map(|t| if flag(args, "--descend") { descend(&problem, &t.spins).spins } else { t.spins.clone() }).collect(); let scores: Vec = finals.iter().map(|s| score(s)).collect(); let best_k = (0..scores.len()).fold(0, |b, k| if (maximise && scores[k] > scores[b]) || (!maximise && scores[k] < scores[b]) { k } else { b }); let best = scores[best_k]; let mean = scores.iter().sum::() / scores.len() as f64; let target: Option = opt(args, "--target", "").map(|t| t.parse().unwrap_or_else(|_| die("bad --target"))); let hit = target.map(|t| scores.iter().filter(|&&v| if maximise { v >= t - 1e-9 } else { v <= t + 1e-9 }).count() as f64 / scores.len() as f64); let bits: String = finals[best_k].iter().map(|&s| if s > 0 { '0' } else { '1' }).collect(); let what = if maximise { "cut" } else { "energy" }; if flag(args, "--json") { let hit_s = hit.map_or("null".into(), |h| h.to_string()); let spins_s = if flag(args, "--spins") { format!(",\"spins\":\"{bits}\"") } else { String::new() }; println!("{{\"best\":{best},\"mean\":{mean},\"method\":\"{method}\",\"n\":{},\"objective\":\"{what}\",\"success\":{hit_s},\"trials\":{}{spins_s}}}", problem.n(), scores.len()); } else { println!("{} spins, {method}, {} trial(s)", problem.n(), scores.len()); println!("best {what} {best}"); println!("mean {what} {mean:.3}"); if let (Some(t), Some(h)) = (target, hit) { println!("reached {t} {:.2}% of trials", 100.0 * h); } if flag(args, "--spins") { println!("spins {bits} (1 = −1; re-check: one pass over the edges)"); } } } // qLDPC memory: a bivariate-bicycle code's depth-8 syndrome cycle, // sampled with circuit-level noise and decoded by BP-OSD. "ldpc" => { use wai_quantum::quantum_frame::Basis; use wai_quantum::quantum_ldpc::{memory_circuit, memory_experiment, BbCode, BpOsdConfig}; let code = match opt(args, "--code", "").as_deref().unwrap_or("72") { "72" => BbCode::bb72(), "90" => BbCode::bb90(), "108" => BbCode::bb108(), "144" => BbCode::bb144(), "288" => BbCode::bb288(), _ => die("--code: 72, 90, 108, 144 or 288"), }; let p: f64 = num(args, "--p", "", 0.003); let cycles: u32 = num(args, "--cycles", "", 6); let bases: Vec<(&str, Basis)> = match opt(args, "--basis", "").as_deref().unwrap_or("both") { "z" | "Z" => vec![("z", Basis::Z)], "x" | "X" => vec![("x", Basis::X)], "both" => vec![("z", Basis::Z), ("x", Basis::X)], _ => die("--basis: z, x or both"), }; if flag(args, "--circuit") { print!("{}", memory_circuit(&code, cycles, p, bases[0].1)); return; } let shots: usize = num(args, "--shots", "", 1000); let seed: u64 = num(args, "--seed", "", 0); let threads: usize = num(args, "--threads", "", std::thread::available_parallelism().map_or(1, |n| n.get())); let cfg = BpOsdConfig { max_iter: num(args, "--max-iter", "", 10_000), osd_order: num(args, "--osd-order", "", 7), ..BpOsdConfig::default() }; let mut rows = Vec::new(); for (i, (name, basis)) in bases.iter().enumerate() { let r = memory_experiment(&code, cycles, p, *basis, shots, seed ^ ((i as u64 + 1) << 32), cfg, threads).unwrap_or_else(|e| die(&e.to_string())); rows.push((*name, r)); } let either = 1.0 - rows.iter().map(|(_, r)| 1.0 - r.per_shot()).product::(); let either_cycle = 1.0 - (1.0 - either).powf(1.0 / f64::from(cycles)); if flag(args, "--json") { let items: Vec = rows .iter() .map(|(b, r)| format!("{{\"basis\":\"{b}\",\"bp_converged\":{},\"failures\":{},\"per_cycle\":{},\"per_shot\":{},\"shots\":{}}}", r.converged, r.failures, r.per_cycle(), r.per_shot(), r.shots)) .collect(); println!("{{\"code\":[{},{}],\"cycles\":{cycles},\"either_per_cycle\":{either_cycle},\"either_per_shot\":{either},\"p\":{p},\"runs\":[{}]}}", code.n(), code.k(), items.join(",")); } else { println!("[[{},{}]] bivariate-bicycle memory, p = {p}, {cycles} noisy cycles", code.n(), code.k()); for (b, r) in &rows { println!( "{b}-basis {} failures in {} shots {:.3e} per shot {:.3e} per cycle BP alone {:.1}%", r.failures, r.shots, r.per_shot(), r.per_cycle(), 100.0 * r.converged as f64 / r.shots as f64 ); } if rows.len() == 2 { println!("either {either:.3e} per shot {either_cycle:.3e} per cycle (the two halves as independent)"); } } } // The architecture around the code: yoked storage, cultivated magic // states, CCZ factories and the clock, priced; and the published // factoring machine re-run region by region. "arch" => { use wai_quantum::quantum_arch::*; let json = flag(args, "--json"); match algo.map(String::as_str) { Some("factoring") | None => { let bits: u32 = num(args, "--bits", "", 2048); let params = FactoringParams::published(bits).unwrap_or_else(|| die("--bits: one of 1024, 1536, 2048, 3072, 4096, 6144, 8192")); let e = factoring_estimate(¶ms, &Machine::published()).unwrap_or_else(|e| die(&e.to_string())); let blocks: Vec<(u32, u32, u32, u64)> = e.cold.blocks.iter().map(|(b, k)| (b.rows - 1, b.distance, b.logical, *k)).collect(); if json { let bs: Vec = blocks.iter().map(|(w, d, l, k)| format!("{{\"count\":{k},\"distance\":{d},\"logical\":{l},\"width\":{w}}}")).collect(); println!( "{{\"bits\":{bits},\"ccz_error\":{},\"ccz_rounds\":{},\"cold_blocks\":[{}],\"cold_logical\":{},\"cold_qubits\":{},\"compute_qubits\":{},\"days\":{},\"expected_shots\":{},\"hot_distance\":{},\"hot_logical\":{},\"hot_qubits\":{},\"hours_per_shot\":{},\"hours_per_shot_rounded\":{},\"physical_qubits\":{},\"shot_survival\":{},\"step_us\":{},\"t_state_d1\":{},\"t_state_volume\":{},\"toffolis_per_shot\":{}}}", e.ccz_error, e.ccz_rounds, bs.join(","), e.cold_logical, e.cold_qubits, e.compute_qubits, e.days, e.expected_shots, e.hot_distance, e.hot_logical, e.hot_qubits, e.hours_per_shot, e.hours_per_shot_rounded, e.physical_qubits(), e.shot_survival, e.timing.step_us(), e.t_state.d1, e.t_state.volume, e.toffolis_per_shot ); } else { println!("{bits}-bit factoring (arXiv:2505.15917 parameters), re-run"); println!("logical qubits {} cold + {} hot", e.cold_logical, e.hot_logical); for (w, d, l, k) in &blocks { println!("cold storage {k} × {w}×{w} grid-yoked block, d = {d}, {l} logical each"); } println!("hot storage d = {}", e.hot_distance); println!("physical qubits {} = {} cold + {} hot + {} compute", e.physical_qubits(), e.cold_qubits, e.hot_qubits, e.compute_qubits); println!("T states cultivated at d1 = {}, {:.0} qubit·rounds each", e.t_state.d1, e.t_state.volume); println!("CCZ states {:.1} rounds per factory, error {:.1e}", e.ccz_rounds, e.ccz_error); println!("clock {:.1} µs per step ({:?}-bound)", e.timing.step_us(), e.timing.bound()); println!("Toffolis {:.4e} per shot", e.toffolis_per_shot); println!("time {:.2} h per shot ({:.2} h at the clock), {:.2} shots, {:.2} days", e.hours_per_shot_rounded, e.hours_per_shot, e.expected_shots, e.days); } } Some("storage") => { let target: f64 = num(args, "--target", "", 1e-15); let model = match opt(args, "--fits", "").as_deref() { None | Some("refit") => StorageModel::uniform_refit(), Some("uniform") => StorageModel::uniform(), Some("si1000") => StorageModel::si1000(), Some(_) => die("--fits: refit, uniform or si1000"), }; let rows = [ ("cold plain", best_block(&model, Storage::Cold, Yoke::None, target, 250)), ("cold rows", best_block(&model, Storage::Cold, Yoke::Rows, target, 250)), ("cold grid", best_block(&model, Storage::Cold, Yoke::Grid, target, 250)), ("hot plain", best_block(&model, Storage::Hot, Yoke::None, target, 250)), ("hot rows", best_block(&model, Storage::Hot, Yoke::Rows, target, 250)), ]; if json { let items: Vec = rows .iter() .filter_map(|(name, b)| b.map(|b| format!("{{\"distance\":{},\"logical\":{},\"name\":\"{name}\",\"qubits_per_logical\":{}}}", b.distance, b.logical, b.qubits_per_logical()))) .collect(); println!("{{\"layouts\":[{}],\"target\":{target}}}", items.join(",")); } else { println!("storage at {target:e} per logical qubit per round"); for (name, b) in rows { match b { Some(b) => println!("{name:<11} d = {:>3} {:>4} × {:<4} patches {:>4} logical {:>8.2} qubits each", b.distance, b.rows, b.cols, b.logical, b.qubits_per_logical()), None => println!("{name:<11} unreachable"), } } } } Some("cultivate") => { let p: f64 = num(args, "--p", "", 1e-3); let target: f64 = num(args, "--target", "", 1e-7); match cultivate(p, target) { Some(t) if json => println!("{{\"d1\":{},\"error\":{},\"volume\":{}}}", t.d1, t.error, t.volume), Some(t) => println!("T state at p = {p:e}: error {:.2e}, {:.0} qubit·rounds (injection distance {})", t.error, t.volume, t.d1), None => die("no cultivation curve reaches that target at that noise (simulated: 5e-4, 1e-3, 2e-3)"), } } Some(other) => die(&format!("arch: unknown '{other}' (factoring, storage, cultivate)")), } } // Noisy circuits: trajectories on the state vector, optionally the // exact density matrix (up to 12 qubits) and cross-entropy fidelity // against the ideal circuit and the digital error model. "noisy" => { use wai_quantum::quantum_noisy::*; let path = opt(args, "--file", "").unwrap_or_else(|| die("--file CIRCUIT.jsonl is required")); let (gates, n) = read_gate_file(&path); let n: u32 = num(args, "--qubits", "-n", n); let noise = NoiseModel { depol1: num(args, "--depol1", "", 0.0), depol2: num(args, "--depol2", "", 0.0), damping: num(args, "--damping", "", 0.0), dephasing: num(args, "--dephasing", "", 0.0), readout: num(args, "--readout", "", 0.0), }; let count: usize = num(args, "--trajectories", "", 1000); let seed: u64 = num(args, "--seed", "", 1); let threads: usize = num(args, "--threads", "", std::thread::available_parallelism().map_or(1, |n| n.get())); let mut names = Vec::new(); let mut observables = Vec::new(); if let Some(spec) = opt(args, "--expect", "") { for term in spec.split(';').map(str::trim).filter(|t| !t.is_empty()) { let string: Vec<(u32, char)> = term .split_whitespace() .map(|f| { let (p, q) = f.split_at(1); let p = p.chars().next().unwrap().to_ascii_uppercase(); if !"XYZ".contains(p) { die("Pauli factors are X, Y or Z followed by a qubit, as in Z0"); } (q.parse().unwrap_or_else(|_| die("bad qubit in --expect")), p) }) .collect(); if string.iter().any(|&(q, _)| q >= n) { die("--expect names a qubit outside the register"); } names.push(term.to_string()); observables.push(vec![(1.0, string)]); } } let run = trajectories(&gates, n, &noise, &observables, count, seed, threads).unwrap_or_else(|e| die(&format!("{e:?}"))); let exact = if flag(args, "--exact") || (flag(args, "--xeb") && n <= MAX_DENSITY_QUBITS) { Some(DensityMatrix::run(&gates, n, &noise).unwrap_or_else(|e| die(&format!("{e:?} (the exact density matrix stops at {MAX_DENSITY_QUBITS} qubits)")))) } else { None }; let xeb = if flag(args, "--xeb") { let mut ideal = wai_quantum::quantum_sv::StateVector::zero(n).unwrap_or_else(|e| die(&format!("{e:?}"))); ideal.run(&gates, 2).unwrap_or_else(|e| die(&format!("{e:?}"))); let sampled = linear_xeb(&ideal, &run.samples); let exact_x = exact.as_ref().map(|rho| exact_xeb(&ideal, rho, noise.readout)); Some((sampled, exact_x, digital_fidelity(&gates, n, &noise), ideal_xeb(&ideal))) } else { None }; if flag(args, "--json") { let mut parts = vec![format!("\"qubits\":{n}"), format!("\"trajectories\":{count}")]; let ex: Vec = names .iter() .enumerate() .map(|(k, t)| { let e = run.expectations[k]; let exact_v = exact.as_ref().map_or("null".into(), |rho| rho.expectation(&observables[k]).to_string()); format!("{{\"exact\":{exact_v},\"mean\":{},\"stderr\":{},\"term\":{}}}", e.mean, e.stderr, serde_json::to_string(t).unwrap()) }) .collect(); parts.push(format!("\"expectations\":[{}]", ex.join(","))); if let Some((s, e, d, i)) = xeb { parts.push(format!("\"xeb\":{{\"digital_model\":{d},\"exact\":{},\"ideal\":{i},\"sampled\":{},\"stderr\":{}}}", e.map_or("null".into(), |v| v.to_string()), s.mean, s.stderr)); } parts.sort(); println!("{{{}}}", parts.join(",")); } else { println!("{n} qubits, {} gates, {count} trajectories", gates.len()); for (k, t) in names.iter().enumerate() { let e = run.expectations[k]; match &exact { Some(rho) => println!("<{t}> = {:+.6} ± {:.6} exact {:+.6}", e.mean, e.stderr, rho.expectation(&observables[k])), None => println!("<{t}> = {:+.6} ± {:.6}", e.mean, e.stderr), } } if let Some((s, e, d, i)) = xeb { println!("XEB noiseless {i:.4} (the circuit's own; 1 when fully scrambling)"); println!("XEB sampled {:.4} ± {:.4} normalised {:.4}", s.mean, s.stderr, s.mean / i); if let Some(v) = e { println!("XEB exact {v:.4} normalised {:.4}", v / i); } println!("digital model {d:.4}"); } } } // The dense state vector, any gate at any angle: run a circuit file // (fused to --fuse qubits), then show the likeliest outcomes, draw // samples, or evaluate Pauli observables. "sv" => { use wai_quantum::quantum_sv::StateVector; let path = opt(args, "--file", "").unwrap_or_else(|| die("--file CIRCUIT.jsonl is required")); let (gates, n) = read_gate_file(&path); let n: u32 = num(args, "--qubits", "-n", n); let width: usize = num(args, "--fuse", "", 2); let mut sv = StateVector::zero(n).unwrap_or_else(|e| die(&format!("{e:?} (at most 32 qubits)"))); sv.run(&gates, width).unwrap_or_else(|e| die(&format!("{e:?}"))); let json = flag(args, "--json"); if let Some(k) = opt(args, "--sample", "") { let k: usize = k.parse().unwrap_or_else(|_| die("bad --sample count")); for x in sv.sample(k, num(args, "--seed", "", 0)) { println!("{:0width$b}", x, width = n as usize); } return; } if let Some(spec) = opt(args, "--expect", "") { // "Z0 X3; Y1 Y2": terms separated by ';', each a product of Pauli factors. let mut values = Vec::new(); for term in spec.split(';').map(str::trim).filter(|t| !t.is_empty()) { let string: Vec<(u32, char)> = term .split_whitespace() .map(|f| { let (p, q) = f.split_at(1); let p = p.chars().next().unwrap().to_ascii_uppercase(); if !"XYZ".contains(p) { die("Pauli factors are X, Y or Z followed by a qubit, as in Z0"); } (q.parse().unwrap_or_else(|_| die("bad qubit in --expect")), p) }) .collect(); if string.iter().any(|&(q, _)| q >= n) { die("--expect names a qubit outside the register"); } values.push((term.to_string(), sv.expectation(&[(1.0, string)]))); } if json { let items: Vec = values.iter().map(|(t, v)| format!("{}:{v}", serde_json::to_string(t).unwrap())).collect(); println!("{{{}}}", items.join(",")); } else { for (t, v) in values { println!("<{t}> = {v:+.12}"); } } return; } let top: usize = num(args, "--top", "", 8); let mut idx: Vec = (0..sv.re.len()).collect(); let p = |x: usize| sv.re[x] * sv.re[x] + sv.im[x] * sv.im[x]; idx.sort_by(|&a, &b| p(b).partial_cmp(&p(a)).unwrap().then(a.cmp(&b))); if json { let rows: Vec = idx .iter() .take(top) .map(|&x| format!("{{\"amplitude\":[{},{}],\"bits\":\"{:0w$b}\",\"p\":{}}}", sv.re[x], sv.im[x], x, p(x), w = n as usize)) .collect(); println!("{{\"gates\":{},\"norm\":{},\"qubits\":{n},\"top\":[{}]}}", gates.len(), sv.norm2(), rows.join(",")); } else { println!("{n} qubits, {} gates (fused to {width}), norm {:.15}", gates.len(), sv.norm2()); for &x in idx.iter().take(top) { println!("|{:0w$b}> p = {:.10} amp = {:+.10} {:+.10}i", x, p(x), sv.re[x], sv.im[x], w = n as usize); } } } "entropy" => { let c = circuit_from(args, algo); let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}"))); let full = DensityMatrix::from_pure(&sv); let rho = match opt(args, "--keep", "") { None => full, Some(list) => { let keep: Vec = list .split(',') .map(|q| q.trim().parse().unwrap_or_else(|_| die("bad --keep qubit"))) .collect(); full.partial_trace(&keep) } }; let s = von_neumann_entropy(&rho); let purity = rho.purity() as f64 / ONE as f64; if flag(args, "--json") { println!( "{{\"qubits\":{},\"entropy_bits\":{s},\"purity\":{purity},\ \"schumacher_qubits_per_symbol\":{}}}", rho.n_qubits, schumacher_limit(&rho) ); } else { println!("qubits {}", rho.n_qubits); println!("purity {purity:.6} (1 = pure, 1/d = maximally mixed)"); println!("entropy {s:.6} bits"); println!("schumacher {:.6} qubits per symbol", schumacher_limit(&rho)); } } // Transport. Push a state through a noisy channel and report what arrives, // including the Holevo bound before and after — which can only fall. "channel" => { let name = algo.map(|s| s.as_str()).unwrap_or("depolarizing"); let p: f64 = num(args, "--prob", "", 0.5); let p_fx = (p * ONE as f64) as i64; let ch = match name { "identity" => Channel::identity(), "bitflip" => Channel::bit_flip(p_fx), "dephasing" => Channel::dephasing(p_fx), "depolarizing" => Channel::depolarizing(p_fx), "damping" => Channel::amplitude_damping(p_fx), other => die(&format!( "unknown channel '{other}' \ (identity, bitflip, dephasing, depolarizing, damping)" )), }; let which = opt(args, "--state", "").unwrap_or_else(|| "plus".into()); let input = one_qubit_state(&which); let out = ch.apply(&input).unwrap_or_else(|| die("channel/state dimension mismatch")); // Tomography is how you learn what arrived; report how far the // reconstruction strays, which is fixed-point rounding and nothing more. let seen = tomography(&out); let mut worst = 0i64; for r in 0..out.d { for c2 in 0..out.d { let (a, b) = (out.get(r, c2), seen.get(r, c2)); worst = worst.max((a.re - b.re).abs()).max((a.im - b.im).abs()); } } // The law, on one line: noise cannot increase what a receiver can read. let zero = one_qubit_state("zero"); let one = one_qubit_state("one"); let before = holevo_bound(&[(ONE / 2, zero.clone()), (ONE / 2, one.clone())]).unwrap(); let after = holevo_bound(&[ (ONE / 2, ch.apply(&zero).unwrap()), (ONE / 2, ch.apply(&one).unwrap()), ]) .unwrap(); let f = |x: i64| x as f64 / ONE as f64; let entropy = von_neumann_entropy(&out); if flag(args, "--json") { println!( "{{\"channel\":\"{name}\",\"p\":{p},\"input\":\"{which}\",\ \"p0\":{},\"p1\":{},\"coherence_re\":{},\"coherence_im\":{},\ \"purity\":{},\"entropy_bits\":{entropy},\ \"holevo_before_bits\":{before},\"holevo_after_bits\":{after},\ \"tomography_max_deviation\":{}}}", f(out.get(0, 0).re), f(out.get(1, 1).re), f(out.get(0, 1).re), f(out.get(0, 1).im), out.purity() as f64 / ONE as f64, f(worst) ); } else { // The identity channel takes no parameter; printing one would imply // a knob that does nothing. if name == "identity" { println!("channel {name}"); } else { println!("channel {name} (p = {p:.3})"); } println!("input |{which}>"); println!("defect {:.9} (0 = probability exactly conserved)", f(ch.trace_preservation_defect())); println!("populations |0> {:.6} |1> {:.6}", f(out.get(0, 0).re), f(out.get(1, 1).re)); println!( "coherence {:.6} {:+.6}i (0 = decohered)", f(out.get(0, 1).re), f(out.get(0, 1).im) ); println!("purity {:.6}", out.purity() as f64 / ONE as f64); println!("entropy {entropy:.6} bits"); println!("holevo {before:.6} -> {after:.6} bits (noise cannot raise it)"); println!("tomography reconstructs to {:.9}", f(worst)); } } // Key distribution. A quiet channel yields a perfect key; an eavesdropper // announces themselves in the error rate, because measuring disturbs. "bb84" => { let rounds: u32 = num(args, "--rounds", "", 2048); let seed: u64 = num(args, "--seed", "", 7); let r = bb84(rounds, seed, flag(args, "--eavesdrop")); if flag(args, "--json") { println!( "{{\"rounds\":{},\"sifted\":{},\"errors\":{},\"qber_ppm\":{},\ \"eavesdropper\":{}}}", r.rounds, r.sifted, r.errors, r.qber_ppm, r.eavesdropper ); } else { println!("rounds {}", r.rounds); println!("sifted {} (bases agreed)", r.sifted); println!("errors {}", r.errors); println!("qber {:.6}", r.qber_ppm as f64 / 1e6); println!("eavesdropper {}", if r.eavesdropper { "yes" } else { "no" }); println!( "no-cloning fidelity {:.6} on a conjugate-basis state \ (1.0 would break the protocol)", clone_fidelity(true) ); } } // The experiment that rules out local hidden variables, computed exactly // rather than sampled — so the violation is arithmetic, not statistics. "chsh" => { let mut c = Circuit::new(2); c.h(0).cx(0, 1); let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}"))); let r = chsh(&sv).unwrap_or_else(|| die("CHSH needs a two-qubit state")); if flag(args, "--json") { println!( "{{\"s\":{},\"classical_bound\":{},\"tsirelson_bound\":{},\ \"violates\":{}}}", r.s, r.classical_bound, r.tsirelson_bound, r.s > r.classical_bound ); } else { println!("state (|00> + |11>)/sqrt(2)"); println!("S {:.6}", r.s); println!("classical {:.6} (any local hidden-variable theory)", r.classical_bound); println!("tsirelson {:.6} (the quantum maximum)", r.tsirelson_bound); println!( "verdict {}", if r.s > r.classical_bound { "violates the classical bound" } else { "no violation" } ); } } // Qutrits. The Fourier gate over Z_3 spreads one basis state evenly across // three, which is the d = 3 analogue of what H does for a qubit. "qutrit" => { let n: u8 = num(args, "--qubits", "-n", 2); let mut q = QuditState::new(3, n).unwrap_or_else(|| die("bad qutrit register")); q.fourier(0); if n > 1 { q.csum(0, 1); } let probs = q.probabilities(); if flag(args, "--json") { let items: Vec = probs.iter().map(|p| format!("{p}")).collect(); println!( "{{\"d\":3,\"qutrits\":{n},\"hash\":\"{}\",\"probabilities\":[{}]}}", hx(&q.hash()), items.join(",") ); } else { println!("d 3 (qutrits)"); println!("register {n}"); println!("circuit F on qutrit 0{}", if n > 1 { ", then CSUM 0->1" } else { "" }); println!("state-hash {}", hx(&q.hash())); println!("outcomes"); for (i, p) in probs.iter().enumerate() { if *p > 1e-9 { let digits: String = (0..n).rev().map(|k| char::from(b'0' + q.digit(i, k))).collect(); println!(" |{digits}> {p:.6}"); } } } } "-h" | "--help" | "help" => usage(), other => die(&format!("unknown command '{other}' (try --help)")), } } #[cfg(test)] mod tests { use super::joules_line; /// The unmetered marker prints as unmetered, never as `0 uJ`; a figure prints as itself. #[test] fn the_unmetered_marker_is_never_printed_as_zero_joules() { assert!(joules_line(0).starts_with("unmetered")); assert!(!joules_line(0).contains("0 uJ")); assert_eq!(joules_line(1_250), "1250 uJ"); } }