From 047e06b88dbf1a8665e90310c71f3c07138def1b Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Fri, 2 Oct 2026 00:26:13 -0400 Subject: [PATCH] wai-quantum 0.3.24: fault-tolerant resource estimation, reproducible and signed Every claim about a fault-tolerant machine (so many physical qubits, so much time) is a resource estimate, and the arithmetic behind one is usually a script run once. quantum_resource makes it a deterministic function of its inputs. - Logical layout and step count by parallel-synthesis sequential Pauli computation; rotation synthesis cost; error budget shared over the parts the program uses. - Planar codes (surface on gate-based and measurement-based qubits, Floquet honeycomb) with distance chosen against the logical budget. - T-state factories by exhaustive search over 15-to-1 unit type, physical or logical rounds, non-decreasing distances and binomially sized copy counts, choosing the smallest space-time volume; the code itself chosen the same way. - The qubits-versus-time frontier. It is held to the planar-ISA model's published tables (arXiv:2211.07629): every distance and factory of the per-qubit factoring table, the distance, factory count, qubits and runtime of all twelve factoring and chemistry summary rows, and all twelve dynamics factory counts. Factoring on the (ns, 1e-4) qubit: d = 13, 18 factories, 8,716,258 qubits, 17 h 43 min. Mutating the model's unit costs, copy-count rule or code selection breaks those tests. Holding it to the tables also showed where the source's prose and numbers disagree (a Toffoli count, two unit durations, a runtime unit), and one printed factory that is not the minimum of its own model; each is documented and settled by the tables. quantum_resource_receipt seals an estimate over its whole job, to the bit, and verify_computation re-runs it and demands every figure back. Receipts and estimate output are byte-identical natively and under wasm32-wasip2, and tests/bit_identity.rs pins the figures of twelve estimates on both. CLI: wai-quantum estimate [factoring|chemistry|--qubits ...] [--qubit ...] [--frontier] [--json] [--receipt]; verify accepts the receipt. repro's sin/cos moved to a gated submodule so the estimator, which needs only ln, builds alone without warnings. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 9 +- wai/wai-quantum/README.md | 35 + wai/wai-quantum/src/lib.rs | 9 +- wai/wai-quantum/src/quantum_resource.rs | 1024 +++++++++++++++++ .../src/quantum_resource_receipt.rs | 584 ++++++++++ wai/wai-quantum/src/repro.rs | 153 +-- wai/wai-quantum/src/repro/trig.rs | 116 ++ wai/wai-quantum/tests/bit_identity.rs | 28 +- wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/Cargo.toml | 5 + wai/wai-rs/src/bin/wai_quantum.rs | 263 +++++ wai/wai-rs/src/lib.rs | 6 + wai/wai-web/Cargo.lock | 2 +- 15 files changed, 2108 insertions(+), 132 deletions(-) create mode 100644 wai/wai-quantum/src/quantum_resource.rs create mode 100644 wai/wai-quantum/src/quantum_resource_receipt.rs create mode 100644 wai/wai-quantum/src/repro/trig.rs diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 5bdd91421..5e9fe6ed9 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -589,7 +589,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.23" +version = "0.3.24" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 8016907e4..e382e5e30 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.3.23" +version = "0.3.24" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 17792bacf..12b401b8e 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.3.23" +version = "0.3.24" edition = "2024" rust-version = "1.97" license = "Apache-2.0" @@ -48,6 +48,11 @@ quantum_spd_receipt = ["quantum_spd", "_signing"] # Tensor-network states on the hardware graph: simple-update evolution, belief- # propagation environments; exact on trees, an independent check on loopy graphs. quantum_bptn = ["quantum_spd"] +# Fault-tolerant resource estimation: logical layout, code distance, exhaustive +# T-factory search and the space-time frontier, reproducible bit for bit. +quantum_resource = [] +# Signed, energy-accounted receipts for it, re-checkable bit for bit anywhere. +quantum_resource_receipt = ["quantum_resource", "_signing"] # Mixed states and the observables that distinguish quantum theory from a classical # account: density matrices, the partial trace, exact Pauli expectations, CHSH. quantum_info = ["quantum"] @@ -87,7 +92,7 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", + "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec", "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom", "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata", diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 9359ff085..9ec24cdb5 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -62,6 +62,8 @@ never creates information. | `quantum_spd` | **sparse Pauli dynamics at utility scale**: arbitrary-angle Pauli rotations on up to 1024 qubits, exact Clifford and quarter-turn rotations, a rigorous bound on what truncation can move a result, and the 127-qubit heavy-hex kicked-Ising circuits — whose exactly computed 5-step magnetization it reproduces on all 127 qubits | | `quantum_spd_receipt` | a signed, energy-accounted receipt over a sparse-Pauli-dynamics run, re-checkable bit for bit on any machine | | `quantum_bptn` | **tensor-network states on the hardware graph**: simple-update evolution with belief-propagation environments and re-gauging, the device's couplings applied a layer at a time in parallel (bit-identical to one at a time). Exact on trees; on the 127-qubit lattice it gives the exactly computed 5-step magnetization to 4e-7 at bond dimension 32 — an independent check on sparse Pauli dynamics, whose errors it does not share | +| `quantum_resource` | **fault-tolerant resource estimation**: from a program's logical counts (qubits, T, Toffoli and rotation gates, measurements) to code distance, T-state factories found by exhaustive search, physical qubits and runtime, on six reference qubits and three codes, with the qubits-versus-time frontier. Reproduces the published worked examples to their printed precision | +| `quantum_resource_receipt` | a signed, energy-accounted receipt over an estimate, re-checkable bit for bit on any machine | | `quantum_receipt` | signed, energy-accounted receipt over a reconstruction | | `quantum_ops` | the operations/attestation layer | | `quantum_cal`, `quantum_control`, `quantum_noise` | calibration, filter-function robust control, DD noise spectroscopy, Cycle-Benchmarking noise learning | @@ -113,6 +115,39 @@ published references. `quantum_bptn` computes the same circuits as tensor networks (`--example utility_bptn`), so a result can be checked by a second method that does not share the first one's approximation. +## Fault-tolerant resources + +What a program would need on a machine that does not exist yet is a resource +estimate, and claims about that machine stand on one. `quantum_resource` makes +the estimate a deterministic function of its inputs: + +```sh +wai-quantum estimate factoring # 2048-bit factoring, six reference qubits +wai-quantum estimate chemistry --qubit ns-e4 --frontier +wai-quantum estimate factoring --qubit ns-e4 --receipt | wai-quantum verify +``` + +The model is the planar-ISA one of arXiv:2211.07629: surface and Floquet-code +patches, lattice surgery, and 15-to-1 distillation factories chosen by +exhaustive search over unit types, distances and copy counts. The tests +reproduce that source's own tables: +- every code distance and factory of its per-qubit factoring table, down to the + copy counts; +- the distance, factory count, physical qubits and runtime of all twelve + factoring and chemistry rows of its summary table; +- all twelve factory counts of its dynamics rows. Factoring on the `(ns, 10⁻⁴)` + qubit comes out at `d = 13`, 18 factories, 8,716,258 physical qubits and + 17 h 43 min. + +Holding the model to those tables exposed places where the source's prose and +its numbers disagree, and one factory that is not the minimum of its own model. +Each is documented in the module and settled by the tables. + +A sealed estimate names its whole job (counts, qubit, code, budget, slowdown and +search limits, to the bit). Re-running it anywhere must reproduce every figure, +down to each factory round, and the bytes are the same natively and as a WASI +component. + ## Determinism Simulation is **byte-exact**: amplitudes are dyadic fixed-point, so a circuit diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index 59ef74f18..5a868bd7c 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -26,7 +26,8 @@ mod linalg; feature = "quantum_mps", feature = "quantum_vml", feature = "quantum_phasor", - feature = "quantum_kernel" + feature = "quantum_kernel", + feature = "quantum_resource" ))] mod repro; @@ -84,6 +85,12 @@ pub mod quantum_spd_receipt; #[cfg(feature = "quantum_bptn")] pub mod quantum_bptn; +#[cfg(feature = "quantum_resource")] +pub mod quantum_resource; + +#[cfg(feature = "quantum_resource_receipt")] +pub mod quantum_resource_receipt; + #[cfg(feature = "quantum_info")] pub mod quantum_info; diff --git a/wai/wai-quantum/src/quantum_resource.rs b/wai/wai-quantum/src/quantum_resource.rs new file mode 100644 index 000000000..3b4a00777 --- /dev/null +++ b/wai/wai-quantum/src/quantum_resource.rs @@ -0,0 +1,1024 @@ +//! Fault-tolerant RESOURCE ESTIMATION — `wai.quantum.resource`. +//! +//! Every serious claim about what a fault-tolerant quantum computer will do — +//! "this molecule in a month", "this key in a week", "a hundred thousand qubits +//! suffice" — is a resource estimate. It starts from a logical program (how +//! many qubits, how many T, Toffoli and rotation gates, in how many layers), then +//! adds an error-correcting code, a magic-state factory and a physical qubit. +//! It ends in two numbers: physical qubits and wall-clock time. Funders, auditors +//! and the people who build the machines all argue over those numbers. The +//! arithmetic between the inputs and the outputs is usually a spreadsheet or an +//! unpublished script, run once. +//! +//! This module makes that arithmetic a deterministic, re-runnable function. An +//! estimate is a pure function of its inputs, computed with `+ − × ÷` and +//! in-crate logarithms, so the same inputs give the same bytes on every machine. +//! It can be sealed into a signed receipt that anyone can re-run and check. +//! +//! # The model +//! +//! The planar-ISA model of arXiv:2211.07629 is implemented in full. +//! +//! **Logical layout and time.** The program is compiled by parallel-synthesis +//! sequential Pauli computation: +//! - `Q = 2·Q_alg + ⌈√(8·Q_alg)⌉ + 1` tiles; +//! - at least `C = M + R_U + R_T + T_U·D_U + 3·R_CCZ` logical steps; +//! - `R = T_U·R_U + 4·R_CCZ + R_T` T states, where a rotation costs +//! `T_U = ⌈0.53·log₂(R_U/ε_syn) + 5.3⌉` T gates. +//! +//! **QEC code.** A code has logical error `P(d) = a·(p/p*)^((d+1)/2)` per tile per +//! step, `n(d)` physical qubits per tile and step time `τ(d)`. The distance is the +//! smallest odd `d` with `Q·C·P(d) ≤ ε_log`. +//! +//! **Factories.** A factory is a sequence of 15-to-1 distillation rounds. Each +//! round is either the space-efficient unit or the Reed–Muller-preparation unit, +//! run on physical qubits (first round only) or on a code patch of non-decreasing +//! distance. A round has acceptance `1 − 15·P_T − 356·P` and output error +//! `35·P_T³ + 7.1·P`. Each round gets the fewest copies such that, at the +//! `1%/rounds` binomial quantile, it still feeds the next round. +//! +//! The search is exhaustive. Among all factories meeting the per-state error +//! target `ε_dis/R`, it chooses the one with the smallest space-time volume +//! `n(D)·τ(D)`. It then provisions `⌈R·τ(D)/(M(D)·t)⌉` copies, where `t` is the +//! algorithm's runtime. +//! +//! **Budget.** The error budget is split evenly between logical failures, +//! distillation and synthesis. A part the program does not use (no rotations, +//! no T states) is not reserved; its share goes to the other parts. +//! +//! # It reproduces the published tables +//! +//! The model is held to the source's own worked examples. +//! +//! - **Factory and distance tables:** every distance and factory is reproduced, +//! down to the copy counts. +//! - **Summary table, factoring and chemistry rows:** the code distance, factory +//! count, physical-qubit count and runtime are reproduced to the printed +//! precision. One example: 2048-bit factoring on the `(ns, 10⁻⁴)` qubit gives +//! `d = 13`, 18 factories of 5,760 qubits, 8.72 M physical qubits and 17 h 43 min. +//! +//! The tests in this module assert each row. Doing so exposed three places where +//! the source's prose and its tables disagree; in each case the tables win. +//! +//! - **Factoring Toffoli count.** It is printed as `3.73·10¹⁰`, but every +//! downstream number in the source (its `C`, `R`, runtime and factory count) +//! requires `3.73·10⁹`. +//! - **Logical Reed–Muller unit.** It is described as taking 13 logical steps, +//! but the source's factory durations (57.2 µs at `d = 13`, 128 µs on the +//! `(ns, 10⁻³)` qubit) require 11. +//! - **Physical distillation units.** They are described as taking 46 and 23 +//! measurement times, but the source's factory counts on measurement-based +//! qubits require 45 and 24. +//! - **Quantum-dynamics example.** Its algorithm-level counts as printed cannot +//! yield its own `C ≈ 1.5·10⁵` and `R ≈ 2.4·10⁶`. Its rows are therefore +//! reproduced from ISA-level counts ([`IsaCounts`]), which this module accepts +//! directly. Its twelve factory counts pin the unprinted step count to +//! `C ≈ 157,000`. At that count all twelve factory counts and code distances +//! match. +//! +//! Two rows of the source cannot both be right. Factoring on the `(ns, 10⁻⁶)` +//! measurement-based qubit is one. +//! - **The per-qubit table** prints a 16,416-qubit, 23 µs factory. +//! - **The summary table** has factories taking 1.2% of 6.2 M qubits, about +//! 5,700 each. +//! - **The search here** finds a 2,640-qubit, 21.9 µs factory that meets the +//! same target under the stated model, about a sixth of the printed volume. +//! It matches the summary row's distance, factory count, qubits and runtime. +//! +//! The source's chemistry runtime on the `(ns, 10⁻⁴)` measurement-based qubit is +//! printed as "24 mins". Its own step count and step time give 24 *days*. +//! +//! # Honest boundaries +//! +//! - **The source's assumptions, stated as assumptions.** These include +//! uniform, independent circuit noise; free transport of T states; and +//! smooth T consumption. +//! - **Planar codes only.** Surface and Floquet-honeycomb patches with lattice +//! surgery. Block (qLDPC) codes, magic-state cultivation and reaction-time- +//! limited scheduling need models of their own. +//! - **Rounding.** The source prints its times truncated (`163.8 ms` as +//! `163 ms`). This module reports exact values; the tests compare at the +//! printed precision. + +use crate::repro::ln; + +/// Which native entangling operations a physical qubit offers. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum InstructionSet { + /// Unitary two-qubit gates (CNOT/CZ) plus single-qubit measurement. + GateBased, + /// Joint Pauli measurements as the native entangling operation. + Measurement, +} + +/// A physical qubit, at the level the estimate needs: operation times and the +/// worst-case Clifford and T error rates. +#[derive(Clone, Debug, PartialEq)] +pub struct PhysicalQubit { + pub name: String, + pub instruction_set: InstructionSet, + /// Gate time in nanoseconds (for a measurement-based qubit, the physical + /// T-gate time). + pub gate_ns: f64, + /// Measurement time in nanoseconds. + pub meas_ns: f64, + /// Error rate of every Clifford operation, including idling and + /// measurement. + pub clifford_error: f64, + /// Error rate of a physical T state. + pub t_error: f64, +} + +impl PhysicalQubit { + fn preset(name: &str, set: InstructionSet, gate_ns: f64, meas_ns: f64, p: f64, pt: f64) -> PhysicalQubit { + PhysicalQubit { name: name.into(), instruction_set: set, gate_ns, meas_ns, clifford_error: p, t_error: pt } + } + /// Gate-based, 100 µs operations, Clifford error 10⁻³ (T 10⁻⁶). + pub fn us_e3() -> PhysicalQubit { + Self::preset("(us, 1e-3)", InstructionSet::GateBased, 100_000.0, 100_000.0, 1e-3, 1e-6) + } + /// Gate-based, 100 µs operations, Clifford error 10⁻⁴ (T 10⁻⁶). + pub fn us_e4() -> PhysicalQubit { + Self::preset("(us, 1e-4)", InstructionSet::GateBased, 100_000.0, 100_000.0, 1e-4, 1e-6) + } + /// Gate-based, 50 ns gates and 100 ns measurement, error 10⁻³. + pub fn ns_e3() -> PhysicalQubit { + Self::preset("(ns, 1e-3)", InstructionSet::GateBased, 50.0, 100.0, 1e-3, 1e-3) + } + /// Gate-based, 50 ns gates and 100 ns measurement, error 10⁻⁴. + pub fn ns_e4() -> PhysicalQubit { + Self::preset("(ns, 1e-4)", InstructionSet::GateBased, 50.0, 100.0, 1e-4, 1e-4) + } + /// Measurement-based, 100 ns operations, Clifford error 10⁻⁴, T error 5%. + pub fn meas_ns_e4() -> PhysicalQubit { + Self::preset("(ns, 1e-4)*", InstructionSet::Measurement, 100.0, 100.0, 1e-4, 0.05) + } + /// Measurement-based, 100 ns operations, Clifford error 10⁻⁶, T error 1%. + pub fn meas_ns_e6() -> PhysicalQubit { + Self::preset("(ns, 1e-6)*", InstructionSet::Measurement, 100.0, 100.0, 1e-6, 0.01) + } + /// The six reference qubits, in the source's table order. + pub fn presets() -> Vec { + vec![Self::us_e3(), Self::us_e4(), Self::ns_e3(), Self::ns_e4(), Self::meas_ns_e4(), Self::meas_ns_e6()] + } +} + +/// A planar QEC code model: `P(d) = a·(p/p*)^((d+1)/2)`, +/// `n(d) = q2·d² + q1·d + q0`, `τ(d) = (g·t_gate + m·t_meas)·d`. +#[derive(Clone, Debug, PartialEq)] +pub struct QecModel { + pub name: String, + pub instruction_set: InstructionSet, + /// Crossing prefactor `a`. + pub prefactor: f64, + /// Threshold `p*`. + pub threshold: f64, + /// `(q2, q1, q0)`. + pub qubits: (u64, u64, i64), + /// `(g, m)`: gate and measurement times per unit of distance. + pub cycle: (f64, f64), +} + +impl QecModel { + /// Rotated surface code on gate-based qubits. + pub fn surface_gate() -> QecModel { + QecModel { + name: "surface code".into(), + instruction_set: InstructionSet::GateBased, + prefactor: 0.03, + threshold: 0.01, + qubits: (2, 0, 0), + cycle: (4.0, 2.0), + } + } + /// Surface code on measurement-based qubits. + pub fn surface_measurement() -> QecModel { + QecModel { + name: "surface code".into(), + instruction_set: InstructionSet::Measurement, + prefactor: 0.08, + threshold: 0.0015, + qubits: (2, 0, 0), + cycle: (0.0, 20.0), + } + } + /// Floquet honeycomb code on measurement-based qubits (arXiv:2107.02194). + pub fn floquet() -> QecModel { + QecModel { + name: "floquet code".into(), + instruction_set: InstructionSet::Measurement, + prefactor: 0.07, + threshold: 0.01, + qubits: (4, 8, -8), + cycle: (0.0, 3.0), + } + } + /// The codes the source considers for a given qubit. + pub fn for_qubit(q: &PhysicalQubit) -> Vec { + match q.instruction_set { + InstructionSet::GateBased => vec![Self::surface_gate()], + InstructionSet::Measurement => vec![Self::surface_measurement(), Self::floquet()], + } + } + + /// Physical qubits per tile at distance `d`. + pub fn tile_qubits(&self, d: u32) -> u64 { + let d = d as i64; + let (a, b, c) = self.qubits; + (a as i64 * d * d + b as i64 * d + c) as u64 + } + + /// Duration of one logical step at distance `d`, in nanoseconds. + pub fn step_ns(&self, q: &PhysicalQubit, d: u32) -> f64 { + (self.cycle.0 * q.gate_ns + self.cycle.1 * q.meas_ns) * d as f64 + } + + /// Probability that one tile fails during one logical step. + pub fn logical_error(&self, p: f64, d: u32) -> f64 { + self.prefactor * powu(p / self.threshold, d.div_ceil(2)) + } +} + +/// `x^k` by repeated squaring in a fixed order: the same bits everywhere. +fn powu(x: f64, mut k: u32) -> f64 { + let (mut base, mut acc) = (x, 1.0); + while k > 0 { + if k & 1 == 1 { + acc *= base; + } + base *= base; + k >>= 1; + } + acc +} + +fn log2(x: f64) -> f64 { + ln(x) / core::f64::consts::LN_2 +} + +/// `⌈√x⌉` in integers. +fn ceil_sqrt(x: u64) -> u64 { + if x == 0 { + return 0; + } + // Newton from above, then fix up: exact for every u64 this module meets. + let mut r = (x as f64).sqrt() as u64; + while r.saturating_mul(r) < x { + r += 1; + } + while r > 0 && (r - 1).saturating_mul(r - 1) >= x { + r -= 1; + } + r +} + +/// What a logical program needs, before layout. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct LogicalCounts { + /// Algorithm qubits. + pub qubits: u64, + /// Pauli measurements. + pub measurements: u64, + /// Arbitrary-angle single-qubit rotations. + pub rotations: u64, + /// Non-Clifford layers containing at least one arbitrary-angle rotation. + pub rotation_depth: u64, + /// T gates. + pub t_gates: u64, + /// Toffoli or CCZ gates. + pub ccz: u64, +} + +/// What the logical layout needs: tiles, minimum steps, T states. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct IsaCounts { + /// Tiles `Q` (algorithm patches plus routing ancillas). + pub tiles: u64, + /// The fewest logical steps the program can run in, `C_min`. + pub min_steps: u64, + /// T states consumed, `R`. + pub t_states: u64, + /// T gates per synthesized rotation (0 when there are no rotations). + pub t_per_rotation: u64, + /// Whether rotations are synthesized (and so share the error budget). + pub synthesizes: bool, +} + +/// How the error budget is shared out. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Budget { + pub logical: f64, + pub distillation: f64, + pub synthesis: f64, +} + +impl Budget { + /// Even split over the parts the program uses. + pub fn split(total: f64, uses_t_states: bool, synthesizes: bool) -> Budget { + let parts = 1.0 + uses_t_states as u8 as f64 + synthesizes as u8 as f64; + let share = total / parts; + Budget { + logical: share, + distillation: if uses_t_states { share } else { 0.0 }, + synthesis: if synthesizes { share } else { 0.0 }, + } + } +} + +/// Lay a logical program out (parallel-synthesis sequential Pauli computation). +pub fn layout(c: &LogicalCounts, total_budget: f64) -> IsaCounts { + let tiles = 2 * c.qubits + ceil_sqrt(8 * c.qubits) + 1; + let synthesizes = c.rotations > 0; + let uses_t = synthesizes || c.t_gates > 0 || c.ccz > 0; + let b = Budget::split(total_budget, uses_t, synthesizes); + let t_per_rotation = if synthesizes { + let x = 0.53 * log2(c.rotations as f64 / b.synthesis) + 5.3; + ceil_pos(x) + } else { + 0 + }; + let min_steps = c.measurements + c.rotations + c.t_gates + t_per_rotation * c.rotation_depth + 3 * c.ccz; + let t_states = t_per_rotation * c.rotations + 4 * c.ccz + c.t_gates; + IsaCounts { tiles, min_steps, t_states, t_per_rotation, synthesizes } +} + +fn ceil_pos(x: f64) -> u64 { + let t = x as u64; + if (t as f64) < x { t + 1 } else { t } +} + +/// The two 15-to-1 distillation units. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum UnitKind { + /// 20 tiles (12 physical qubits), 13 steps (45 measurement times). + SpaceEfficient, + /// 31 tiles (31 physical qubits), 11 steps (24 measurement times). + ReedMuller, +} + +impl UnitKind { + const ALL: [UnitKind; 2] = [UnitKind::SpaceEfficient, UnitKind::ReedMuller]; + fn physical(self) -> (u64, f64) { + match self { + UnitKind::SpaceEfficient => (12, 45.0), + UnitKind::ReedMuller => (31, 24.0), + } + } + fn logical(self) -> (u64, u64) { + match self { + UnitKind::SpaceEfficient => (20, 13), + UnitKind::ReedMuller => (31, 11), + } + } +} + +/// One distillation round of a factory. +#[derive(Clone, Debug, PartialEq)] +pub struct Round { + pub kind: UnitKind, + /// Code distance; `1` means the round runs on bare physical qubits. + pub distance: u32, + /// Copies of the unit run in parallel. + pub copies: u64, + /// Physical qubits of one copy. + pub unit_qubits: u64, + /// Duration of the round, nanoseconds. + pub duration_ns: f64, + /// Probability that one copy is rejected. + pub failure: f64, + /// Error rate of the T states this round outputs. + pub output_error: f64, +} + +/// A T-state factory. +#[derive(Clone, Debug, PartialEq)] +pub struct Factory { + pub rounds: Vec, + /// Physical qubits: the largest round (rounds run one after another). + pub qubits: u64, + pub duration_ns: f64, + pub output_error: f64, + /// T states output in at least 99% of runs. + pub outputs: u64, +} + +/// Smallest `k` with `P(X ≤ k) ≥ alpha`, `X ~ Binomial(n, s)`. Weights are +/// grown outward from the mode in a fixed order, so the answer is the same +/// bits everywhere, and no factorial or power overflows. +fn binomial_quantile(n: u64, s: f64, alpha: f64) -> u64 { + if s >= 1.0 { + return n; + } + if s <= 0.0 { + return 0; + } + let q = 1.0 - s; + let mode = (((n + 1) as f64) * s) as u64; + let mode = mode.min(n); + // w[i - lo] ∝ P(X = i) for i in lo..=hi. + let mut down = Vec::new(); + let mut w = 1.0; + let mut i = mode; + while i > 0 { + w *= (i as f64) / ((n - i + 1) as f64) * (q / s); + if w < 1e-40 { + break; + } + i -= 1; + down.push(w); + } + let lo = mode - down.len() as u64; + let mut up = Vec::new(); + let mut w = 1.0; + let mut i = mode; + while i < n { + w *= ((n - i) as f64) / ((i + 1) as f64) * (s / q); + if w < 1e-40 { + break; + } + i += 1; + up.push(w); + } + let weights: Vec = down.iter().rev().copied().chain(core::iter::once(1.0)).chain(up.iter().copied()).collect(); + let total: f64 = weights.iter().sum(); + let mut cum = 0.0; + for (j, w) in weights.iter().enumerate() { + cum += w; + if cum >= alpha * total { + return lo + j as u64; + } + } + lo + weights.len() as u64 - 1 +} + +/// Fewest copies of a unit (one output, success `s`) such that at the +/// `alpha` quantile at least `needed` outputs survive. +fn copies_for(needed: u64, s: f64, alpha: f64) -> Option { + let ok = |n: u64| binomial_quantile(n, s, alpha) >= needed; + let mut n = needed.max(1); + while !ok(n) { + n *= 2; + if n >= 1_000_000_000 { + return None; + } + } + let (mut lo, mut hi) = (n / 2, n); + while lo < hi { + let mid = lo + (hi - lo) / 2; + if ok(mid) { hi = mid } else { lo = mid + 1 } + } + Some(hi) +} + +/// Limits for the exhaustive factory search. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct SearchLimits { + pub max_distance: u32, + pub max_rounds: usize, +} + +impl Default for SearchLimits { + fn default() -> Self { + SearchLimits { max_distance: 51, max_rounds: 3 } + } +} + +/// Build the factory for a fixed sequence of `(unit, distance)` rounds: +/// propagate error rates forward, then size copies backward. Returns `None` +/// for an invalid sequence, one that misses `target`, or one that cannot beat +/// `bound` in space-time volume even with a single copy per round. +fn build_factory( + seq: &[(UnitKind, u32)], + qubit: &PhysicalQubit, + qec: &QecModel, + target: f64, + bound: f64, +) -> Option { + let alpha = 0.01 / seq.len() as f64; + let mut input = qubit.t_error; + let mut rounds = Vec::with_capacity(seq.len()); + for &(kind, d) in seq { + let (clifford, unit_qubits, duration_ns) = if d == 1 { + let (nq, steps) = kind.physical(); + (qubit.clifford_error, nq, steps * qubit.meas_ns) + } else { + let (tiles, steps) = kind.logical(); + ( + qec.logical_error(qubit.clifford_error, d), + tiles * qec.tile_qubits(d), + steps as f64 * qec.step_ns(qubit, d), + ) + }; + let failure = 15.0 * input + 356.0 * clifford; + if failure <= 0.0 || failure >= 1.0 { + return None; + } + let output_error = 35.0 * input * input * input + 7.1 * clifford; + if output_error > input { + return None; + } + rounds.push(Round { kind, distance: d, copies: 1, unit_qubits, duration_ns, failure, output_error }); + input = output_error; + } + if input > target { + return None; + } + let floor_qubits = rounds.iter().map(|r| r.unit_qubits).max().unwrap_or(0); + let duration_ns: f64 = rounds.iter().map(|r| r.duration_ns).sum(); + if floor_qubits as f64 * duration_ns > bound { + return None; + } + let mut needed = 1; + for r in rounds.iter_mut().rev() { + r.copies = copies_for(needed, 1.0 - r.failure, alpha)?; + needed = 15 * r.copies; + } + let qubits = rounds.iter().map(|r| r.copies * r.unit_qubits).max().unwrap_or(0); + Some(Factory { rounds, qubits, duration_ns, output_error: input, outputs: 1 }) +} + +/// The factory with the smallest space-time volume whose T states meet +/// `target`, or `None` if no factory within `limits` does. +pub fn find_factory(qubit: &PhysicalQubit, qec: &QecModel, target: f64, limits: SearchLimits) -> Option { + let mut distances: Vec = vec![1]; + distances.extend((3..=limits.max_distance).step_by(2)); + let mut best: Option<(f64, Factory)> = None; + let mut seq: Vec<(UnitKind, u32)> = Vec::with_capacity(limits.max_rounds); + for rounds in 1..=limits.max_rounds { + search(&mut seq, rounds, &distances, qubit, qec, target, &mut best); + } + best.map(|(_, f)| f) +} + +fn search( + seq: &mut Vec<(UnitKind, u32)>, + rounds: usize, + distances: &[u32], + qubit: &PhysicalQubit, + qec: &QecModel, + target: f64, + best: &mut Option<(f64, Factory)>, +) { + if seq.len() == rounds { + let bound = best.as_ref().map_or(f64::INFINITY, |(v, _)| *v); + let Some(f) = build_factory(seq, qubit, qec, target, bound) else { return }; + let volume = f.qubits as f64 * f.duration_ns; + let better = match best { + None => true, + Some((v, b)) => volume < *v || (volume == *v && f.qubits < b.qubits), + }; + if better { + *best = Some((volume, f)); + } + return; + } + let floor = seq.last().map_or(1, |&(_, d)| d.max(3)); + for &d in distances { + // Physical rounds only first; distances never shrink. + if (d == 1 && !seq.is_empty()) || (d != 1 && d < floor) { + continue; + } + for kind in UnitKind::ALL { + seq.push((kind, d)); + search(seq, rounds, distances, qubit, qec, target, best); + seq.pop(); + } + } +} + +/// How fast to run: at the minimum step count, or slowed by a factor to +/// trade factories for time. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct EstimateConfig { + pub error_budget: f64, + /// Run in `slowdown × C_min` logical steps (`1` = as fast as possible). + pub slowdown: u64, + pub limits: SearchLimits, +} + +impl EstimateConfig { + pub fn new(error_budget: f64) -> EstimateConfig { + EstimateConfig { error_budget, slowdown: 1, limits: SearchLimits::default() } + } +} + +/// A complete physical estimate. +#[derive(Clone, Debug, PartialEq)] +pub struct Estimate { + pub qubit: String, + pub code: String, + pub isa: IsaCounts, + pub budget: Budget, + /// Code distance of the algorithm's tiles. + pub distance: u32, + pub tile_qubits: u64, + pub step_ns: f64, + /// Logical steps actually run. + pub steps: u64, + pub runtime_ns: f64, + /// Per-tile, per-step logical error at the chosen distance. + pub logical_error: f64, + /// Per-state error the T states must meet. + pub t_target: f64, + pub factory: Option, + pub factories: u64, + pub algorithm_qubits: u64, + pub factory_qubits: u64, + pub physical_qubits: u64, +} + +impl Estimate { + /// Share of physical qubits spent on T factories. + pub fn factory_fraction(&self) -> f64 { + self.factory_qubits as f64 / self.physical_qubits as f64 + } + + /// Space-time volume: physical qubits × runtime (qubit-nanoseconds). + pub fn volume(&self) -> f64 { + self.physical_qubits as f64 * self.runtime_ns + } +} + +#[derive(Clone, Debug, PartialEq, Eq)] +pub enum EstimateError { + /// The physical error rate is at or above the code's threshold. + AboveThreshold, + /// No odd distance up to the limit meets the logical budget. + DistanceExceeded, + /// No factory within the search limits meets the T-state target. + NoFactory, + /// An empty program, a zero budget or a zero slowdown. + Invalid, + /// The code and the qubit use different instruction sets. + Mismatch, + /// A sealing meter returned without running the estimate. + NotRun, +} + +/// Estimate a logical program on a qubit and code. +pub fn estimate(c: &LogicalCounts, qubit: &PhysicalQubit, qec: &QecModel, cfg: &EstimateConfig) -> Result { + if c.qubits == 0 || cfg.error_budget.is_nan() || cfg.error_budget <= 0.0 { + return Err(EstimateError::Invalid); + } + estimate_isa(&layout(c, cfg.error_budget), qubit, qec, cfg) +} + +/// Estimate from ISA-level counts directly. +pub fn estimate_isa(isa: &IsaCounts, qubit: &PhysicalQubit, qec: &QecModel, cfg: &EstimateConfig) -> Result { + if isa.tiles == 0 || isa.min_steps == 0 || cfg.slowdown == 0 || !(cfg.error_budget > 0.0 && cfg.error_budget < 1.0) { + return Err(EstimateError::Invalid); + } + if qec.instruction_set != qubit.instruction_set { + return Err(EstimateError::Mismatch); + } + if qubit.clifford_error >= qec.threshold { + return Err(EstimateError::AboveThreshold); + } + let budget = Budget::split(cfg.error_budget, isa.t_states > 0, isa.synthesizes); + let t_target = if isa.t_states > 0 { budget.distillation / isa.t_states as f64 } else { 0.0 }; + let factory = if isa.t_states == 0 { + None + } else if t_target >= qubit.t_error { + // Raw T states already meet the target: no distillation. + None + } else { + Some(find_factory(qubit, qec, t_target, cfg.limits).ok_or(EstimateError::NoFactory)?) + }; + let mut steps = isa.min_steps.saturating_mul(cfg.slowdown); + loop { + let distance = distance_for(qec, qubit, isa.tiles, steps, budget.logical, cfg.limits.max_distance)?; + let step_ns = qec.step_ns(qubit, distance); + let runtime_ns = step_ns * steps as f64; + // The program must outlast one factory run. + if let Some(f) = &factory + && runtime_ns < f.duration_ns + { + let need = ceil_pos(f.duration_ns / step_ns); + if need > steps { + steps = need; + continue; + } + } + let (factories, factory_qubits) = match &factory { + Some(f) => { + let n = ceil_pos(isa.t_states as f64 * f.duration_ns / (f.outputs as f64 * runtime_ns)); + (n, n * f.qubits) + } + None => (0, 0), + }; + let tile_qubits = qec.tile_qubits(distance); + let algorithm_qubits = isa.tiles * tile_qubits; + return Ok(Estimate { + qubit: qubit.name.clone(), + code: qec.name.clone(), + isa: *isa, + budget, + distance, + tile_qubits, + step_ns, + steps, + runtime_ns, + logical_error: qec.logical_error(qubit.clifford_error, distance), + t_target, + factory, + factories, + algorithm_qubits, + factory_qubits, + physical_qubits: algorithm_qubits + factory_qubits, + }); + } +} + +/// Smallest odd distance with `tiles · steps · P(d) ≤ budget`. +fn distance_for(qec: &QecModel, qubit: &PhysicalQubit, tiles: u64, steps: u64, budget: f64, max: u32) -> Result { + let volume = tiles as f64 * steps as f64; + let mut d = 1; + while d <= max { + if volume * qec.logical_error(qubit.clifford_error, d) <= budget { + return Ok(d); + } + d += 2; + } + Err(EstimateError::DistanceExceeded) +} + +/// The cheapest estimate over every code the qubit supports, by space-time +/// volume (physical qubits × runtime). +pub fn estimate_best(c: &LogicalCounts, qubit: &PhysicalQubit, cfg: &EstimateConfig) -> Result { + let isa = layout(c, cfg.error_budget); + estimate_best_isa(&isa, qubit, cfg) +} + +/// [`estimate_best`] from ISA-level counts. +pub fn estimate_best_isa(isa: &IsaCounts, qubit: &PhysicalQubit, cfg: &EstimateConfig) -> Result { + let mut best: Option = None; + let mut last_err = EstimateError::Mismatch; + for qec in QecModel::for_qubit(qubit) { + match estimate_isa(isa, qubit, &qec, cfg) { + Ok(e) => { + if best.as_ref().is_none_or(|b| e.volume() < b.volume()) { + best = Some(e); + } + } + Err(e) => last_err = e, + } + } + best.ok_or(last_err) +} + +/// The space-time frontier: estimates at slowdowns `1, 2, 4, …` up to +/// `max_slowdown`, keeping only those that use fewer qubits than every faster +/// one. +pub fn frontier(isa: &IsaCounts, qubit: &PhysicalQubit, qec: &QecModel, cfg: &EstimateConfig, max_slowdown: u64) -> Vec { + let mut out: Vec = Vec::new(); + let mut s = 1; + while s <= max_slowdown { + let c = EstimateConfig { slowdown: s, ..*cfg }; + if let Ok(e) = estimate_isa(isa, qubit, qec, &c) + && out.last().is_none_or(|p| e.physical_qubits < p.physical_qubits) + { + out.push(e); + } + s *= 2; + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + + const US: f64 = 1e3; + const MS: f64 = 1e6; + const S: f64 = 1e9; + const HOUR: f64 = 3600.0 * S; + const DAY: f64 = 24.0 * HOUR; + + /// The chemistry program (double-factorized qubitization), ε = 0.01. + fn chemistry() -> LogicalCounts { + 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, + } + } + + /// 2048-bit factoring, ε = 1/3, with the Toffoli count every downstream + /// number of the source requires (3.73·10⁹; it prints 3.73·10¹⁰). + fn factoring() -> LogicalCounts { + LogicalCounts { + qubits: 12_581, + measurements: 1_080_000_000, + rotations: 12, + rotation_depth: 12, + t_gates: 12, + ccz: 3_730_000_000, + } + } + + /// Whether `x` matches a value printed to `figures` significant figures. + /// The source rounds some figures and truncates others, so anything within + /// one unit of the last printed digit matches. + fn printed(x: f64, shown: f64, figures: i32) -> bool { + let unit = 10f64.powi(shown.abs().log10().floor() as i32 - figures + 1); + (x - shown).abs() < unit + } + + #[test] + fn integer_helpers_are_exact() { + for x in [1u64, 2, 3, 4, 8, 9, 10, 100_648, 10_544, u32::MAX as u64] { + let r = ceil_sqrt(x); + assert!(r * r >= x && (r - 1) * (r - 1) < x, "{x}"); + } + assert_eq!(powu(0.1, 0), 1.0); + assert_eq!(powu(2.0, 10), 1024.0); + assert!((powu(0.01, 5) - 1e-10).abs() < 1e-24); + // A Binomial(16, 1 − 0.002568) has its 0.5% quantile at 15, and a + // Binomial(15, ·) does not reach 15. + assert_eq!(binomial_quantile(16, 1.0 - 0.002568, 0.005), 15); + assert_eq!(binomial_quantile(15, 1.0 - 0.002568, 0.005), 14); + assert_eq!(copies_for(15, 1.0 - 0.002568, 0.005), Some(16)); + } + + #[test] + fn the_layout_reproduces_the_isa_counts() { + let c = layout(&chemistry(), 0.01); + assert_eq!(c.tiles, 2740); + assert_eq!(c.t_per_rotation, 25); + assert!(printed(c.min_steps as f64, 4.12e11, 3)); // printed 4.10e11, from unrounded inputs + assert!((c.t_states as f64 / 5.44e11 - 1.0).abs() < 0.005); // inputs printed to 3 s.f. + let f = layout(&factoring(), 1.0 / 3.0); + assert_eq!(f.tiles, 25_481); + assert_eq!(f.t_per_rotation, 9); + assert!(printed(f.min_steps as f64, 1.23e10, 3)); + assert!(printed(f.t_states as f64, 1.49e10, 3)); + } + + #[test] + fn the_worked_factories_are_reproduced_exactly() { + let q = PhysicalQubit::ns_e4(); + let s = QecModel::surface_gate(); + // One round, d = 9: 3240 qubits, 46.8 µs, 5.6·10⁻¹¹. + let f = build_factory(&[(UnitKind::SpaceEfficient, 9)], &q, &s, 1.0, f64::INFINITY).unwrap(); + assert_eq!((f.qubits, f.rounds[0].copies), (3240, 1)); + assert!((f.duration_ns - 46.8 * US).abs() < 1e-6); + assert!(printed(f.output_error, 5.6e-11, 2)); + // 16 × (d = 3) → 1 × (d = 11): 5760 qubits, 72.8 µs, 5.51·10⁻¹³. + let f = build_factory(&[(UnitKind::SpaceEfficient, 3), (UnitKind::SpaceEfficient, 11)], &q, &s, 1.0, f64::INFINITY).unwrap(); + assert_eq!(f.rounds.iter().map(|r| r.copies).collect::>(), vec![16, 1]); + assert_eq!(f.qubits, 5760); + assert!((f.duration_ns - 72.8 * US).abs() < 1e-6); + assert!(printed(f.output_error, 5.51e-13, 3)); + // 16 × (d = 5) → Reed–Muller (d = 13): 16000 qubits, 83.2 µs, 2.13·10⁻¹⁵. + let f = build_factory(&[(UnitKind::SpaceEfficient, 5), (UnitKind::ReedMuller, 13)], &q, &s, 1.0, f64::INFINITY).unwrap(); + assert_eq!(f.qubits, 16_000); + assert_eq!(f.rounds[1].unit_qubits, 10_478); + assert!((f.duration_ns - 83.2 * US).abs() < 1e-6); + assert!(printed(f.output_error, 2.13e-15, 3)); + } + + /// Distance and factory for factoring on all six qubits (the source's + /// per-qubit table): `P(d) ≤ 3.5·10⁻¹⁶`, `P_T ≤ 7.4·10⁻¹²`. + #[test] + fn the_per_qubit_factoring_table_is_reproduced() { + let isa = layout(&factoring(), 1.0 / 3.0); + let cfg = EstimateConfig::new(1.0 / 3.0); + // (code, d, n(d), τ(d) truncated, n(D), τ(D) truncated) + let rows: [(&str, u32, u64, f64, u64, f64); 6] = [ + ("surface code", 27, 1458, 16.0 * MS, 17_640, 163.0 * MS), + ("surface code", 13, 338, 7.0 * MS, 4840, 85.0 * MS), + ("surface code", 27, 1458, 10.0 * US, 33_320, 128.0 * US), + ("surface code", 13, 338, 5.0 * US, 5760, 72.0 * US), + ("floquet code", 15, 1012, 4.0 * US, 21_840, 52.0 * US), + ("floquet code", 7, 244, 2.0 * US, 16_416, 23.0 * US), + ]; + let unit = |x: f64| if x >= MS { MS } else { US }; + for (i, (q, row)) in PhysicalQubit::presets().iter().zip(rows).enumerate() { + let e = estimate_best_isa(&isa, q, &cfg).unwrap(); + let f = e.factory.as_ref().unwrap(); + if i == 5 { + // The printed factory is not the minimum of the stated model: + // the one found here meets the same target in a sixth of the + // space-time volume (see the module docs). + assert_eq!((e.code.as_str(), e.distance, e.tile_qubits), (row.0, row.1, row.2)); + assert!(f.output_error <= e.t_target); + assert!(f.qubits as f64 * f.duration_ns < row.4 as f64 * row.5 / 5.0); + continue; + } + let got = ( + e.code.as_str(), + e.distance, + e.tile_qubits, + (e.step_ns / unit(row.3)).floor() * unit(row.3), + f.qubits, + (f.duration_ns / unit(row.5)).floor() * unit(row.5), + ); + assert_eq!(got, row, "{}: {:?}", q.name, f.rounds); + } + } + + /// The summary rows for factoring and chemistry: distance and factory + /// count exactly; physical qubits and runtime to the printed figures. + #[test] + fn the_summary_rows_are_reproduced() { + const YEAR: f64 = 365.25 * DAY; + const MONTH: f64 = YEAR / 12.0; + // (d, factories, qubits in millions, runtime as printed, its unit) + let factoring_rows: [(u32, u64, f64, f64, f64); 6] = [ + (27, 13, 37.0, 6.2, YEAR), + (13, 14, 8.6, 3.0, YEAR), + (27, 15, 37.0, 1.5, DAY), + (13, 18, 8.7, 18.0, HOUR), + (15, 15, 26.0, 15.0, HOUR), + (7, 13, 6.2, 7.1, HOUR), + ]; + let chemistry_rows: [(u32, u64, f64, f64, f64); 6] = [ + (33, 15, 6.4, 260.0, YEAR), + (17, 14, 1.6, 130.0, YEAR), + (33, 17, 6.9, 2.0, MONTH), + (17, 17, 1.9, 1.0, MONTH), + (17, 19, 4.5, 24.0, DAY), // printed "24 mins"; its own C and τ(d) give days + (9, 19, 1.3, 12.0, DAY), + ]; + for (prog, eps, rows) in [(factoring(), 1.0 / 3.0, factoring_rows), (chemistry(), 0.01, chemistry_rows)] { + let cfg = EstimateConfig::new(eps); + for (q, (d, nf, mq, t, unit)) in PhysicalQubit::presets().iter().zip(rows) { + let e = estimate_best(&prog, q, &cfg).unwrap(); + assert_eq!((e.distance, e.factories), (d, nf), "{} {:?}", q.name, e.factory); + assert!(printed(e.physical_qubits as f64 / 1e6, mq, 2), "{}: {}", q.name, e.physical_qubits); + assert!(printed(e.runtime_ns / unit, t, 2), "{}: {} vs {t}", q.name, e.runtime_ns / unit); + } + } + } + + #[test] + fn the_factoring_example_matches_to_the_minute() { + let e = estimate(&factoring(), &PhysicalQubit::ns_e4(), &QecModel::surface_gate(), &EstimateConfig::new(1.0 / 3.0)).unwrap(); + assert_eq!(e.distance, 13); + assert_eq!(e.factories, 18); + assert_eq!(e.physical_qubits, 18 * 5760 + 25_481 * 338); + let minutes = (e.runtime_ns / (60.0 * S)) as u64; + assert_eq!((minutes / 60, minutes % 60), (17, 43)); + let f = e.factory.unwrap(); + assert_eq!(f.rounds.iter().map(|r| (r.kind, r.distance, r.copies)).collect::>(), + vec![(UnitKind::SpaceEfficient, 3, 16), (UnitKind::SpaceEfficient, 11, 1)]); + } + + /// The dynamics rows, from ISA-level counts: `Q = 230` and `R = 2.4·10⁶` as + /// printed, and `C = 157,000`, which the twelve printed factory counts pin + /// down (the source prints `C` only as `1.5·10⁵`). Fast and ×10 slowed: + /// distance, factory count and qubits. + #[test] + fn the_dynamics_rows_follow_from_the_isa_counts() { + let isa = IsaCounts { tiles: 230, min_steps: 157_000, t_states: 2_400_000, t_per_rotation: 20, synthesizes: true }; + let fast: [(u32, u64, f64); 6] = [(19, 199, 3.0), (9, 199, 0.68), (19, 242, 8.2), (9, 199, 0.68), (9, 260, 5.8), (5, 224, 0.62)]; + let slow: [(u32, u64, f64); 6] = [(21, 18, 0.46), (11, 17, 0.11), (21, 22, 0.94), (11, 17, 0.11), (11, 22, 0.61), (5, 23, 0.09)]; + for (slowdown, rows) in [(1, fast), (10, slow)] { + let cfg = EstimateConfig { slowdown, ..EstimateConfig::new(0.001) }; + for (q, (d, nf, mq)) in PhysicalQubit::presets().iter().zip(rows) { + let e = estimate_best_isa(&isa, q, &cfg).unwrap(); + assert_eq!((e.distance, e.factories), (d, nf), "{} ×{slowdown}", q.name); + // "0.09M" is printed to one figure, the rest to two. + let figures = if mq < 0.1 { 1 } else { 2 }; + assert!(printed(e.physical_qubits as f64 / 1e6, mq, figures), "{} ×{slowdown}: {}", q.name, e.physical_qubits); + } + } + } + + #[test] + fn the_frontier_trades_qubits_for_time() { + let isa = layout(&chemistry(), 0.01); + let q = PhysicalQubit::ns_e4(); + let fr = frontier(&isa, &q, &QecModel::surface_gate(), &EstimateConfig::new(0.01), 64); + assert!(fr.len() >= 2); + for w in fr.windows(2) { + assert!(w[1].physical_qubits < w[0].physical_qubits); + assert!(w[1].runtime_ns > w[0].runtime_ns); + } + } + + #[test] + fn bad_inputs_are_refused() { + let isa = layout(&factoring(), 1.0 / 3.0); + let mut bad = PhysicalQubit::ns_e3(); + bad.clifford_error = 0.02; + let cfg = EstimateConfig::new(1.0 / 3.0); + assert_eq!(estimate_isa(&isa, &bad, &QecModel::surface_gate(), &cfg), Err(EstimateError::AboveThreshold)); + assert_eq!(estimate_isa(&isa, &PhysicalQubit::ns_e3(), &QecModel::floquet(), &cfg), Err(EstimateError::Mismatch)); + assert_eq!(estimate_isa(&isa, &PhysicalQubit::ns_e3(), &QecModel::surface_gate(), &EstimateConfig { slowdown: 0, ..cfg }), Err(EstimateError::Invalid)); + assert_eq!(estimate(&LogicalCounts::default(), &PhysicalQubit::ns_e3(), &QecModel::surface_gate(), &cfg), Err(EstimateError::Invalid)); + } + + #[test] + fn a_clifford_only_program_needs_no_factory() { + let c = LogicalCounts { qubits: 100, measurements: 10_000, ..Default::default() }; + let e = estimate(&c, &PhysicalQubit::ns_e4(), &QecModel::surface_gate(), &EstimateConfig::new(0.01)).unwrap(); + assert!(e.factory.is_none()); + assert_eq!((e.factories, e.factory_qubits), (0, 0)); + assert_eq!(e.budget.logical, 0.01); + } +} + diff --git a/wai/wai-quantum/src/quantum_resource_receipt.rs b/wai/wai-quantum/src/quantum_resource_receipt.rs new file mode 100644 index 000000000..aa67682b4 --- /dev/null +++ b/wai/wai-quantum/src/quantum_resource_receipt.rs @@ -0,0 +1,584 @@ +//! Signed, re-checkable resource estimates (`wai.quantum.resource`). +//! +//! A resource estimate is the number a claim stands on: so many physical +//! qubits, so much time. The receipt binds that number to everything it came +//! from: +//! +//! - **the job**: the program's counts, the physical qubit, the code (or the +//! rule "cheapest code this qubit supports"), the error budget, the slowdown +//! and the search limits. Rates and times are hashed as their exact IEEE-754 +//! bits. +//! - **the estimate**: every figure in it, down to each factory round's copy +//! count and output error. +//! - **the energy** the sink spent producing it (`joules_micro`, `0` if +//! unmetered). +//! +//! The estimate is a pure function of the job, computed with `+ − × ÷` and +//! in-crate logarithms. [`ResourceReceipt::verify_computation`] therefore +//! re-runs the job anywhere and demands the identical estimate, bit for bit. +//! Change one input (an error rate by one ulp, the budget, the code) and the +//! check fails. +//! +//! The receipt is a JWP profile like the crate's others (Merkle root over leaf +//! hashes, Ed25519). It can carry an energy-class label +//! (`quantum_energy::Labelled`). + +use crate::merkle::merkle_root; +use crate::quantum_resource::{ + estimate_best, estimate_best_isa, estimate, estimate_isa, Budget, Estimate, EstimateConfig, EstimateError, + Factory, InstructionSet, IsaCounts, LogicalCounts, PhysicalQubit, QecModel, Round, SearchLimits, UnitKind, +}; +use ed25519_dalek::{Signature, Signer, SigningKey, Verifier, VerifyingKey}; +use serde_json::{json, Value}; + +const DOMAIN_RECEIPT: &[u8] = b"wai:resource-receipt\x01"; +const DOMAIN_RECEIPT_ID: &[u8] = b"wai:resource-receipt-id\x01"; +const DOMAIN_JOB: &[u8] = b"wai:resource-job\x01"; +const DOMAIN_RESULT: &[u8] = b"wai:resource-estimate\x01"; + +/// The program being estimated, at the level it is known. +#[derive(Clone, Debug, PartialEq)] +pub enum Program { + /// Algorithm-level counts, laid out by the estimator. + Logical(LogicalCounts), + /// Layout-level counts, given directly. + Isa(IsaCounts), +} + +/// Everything that determines an estimate. +#[derive(Clone, Debug, PartialEq)] +pub struct ResourceJob { + pub program: Program, + pub qubit: PhysicalQubit, + /// The code to use, or `None` for the cheapest the qubit supports. + pub code: Option, + pub config: EstimateConfig, +} + +impl ResourceJob { + pub fn run(&self) -> Result { + match (&self.program, &self.code) { + (Program::Logical(c), Some(q)) => estimate(c, &self.qubit, q, &self.config), + (Program::Logical(c), None) => estimate_best(c, &self.qubit, &self.config), + (Program::Isa(i), Some(q)) => estimate_isa(i, &self.qubit, q, &self.config), + (Program::Isa(i), None) => estimate_best_isa(i, &self.qubit, &self.config), + } + } + + /// BLAKE3 over the job's canonical encoding. + pub fn job_hash(&self) -> [u8; 32] { + let mut o = DOMAIN_JOB.to_vec(); + match &self.program { + Program::Logical(c) => { + o.push(1); + for x in [c.qubits, c.measurements, c.rotations, c.rotation_depth, c.t_gates, c.ccz] { + u(x, &mut o); + } + } + Program::Isa(i) => { + o.push(2); + isa(i, &mut o); + } + } + qubit(&self.qubit, &mut o); + match &self.code { + Some(q) => { + o.push(1); + code(q, &mut o); + } + None => o.push(0), + } + f(self.config.error_budget, &mut o); + u(self.config.slowdown, &mut o); + u(self.config.limits.max_distance as u64, &mut o); + u(self.config.limits.max_rounds as u64, &mut o); + *blake3::hash(&o).as_bytes() + } +} + +fn u(x: u64, o: &mut Vec) { + o.extend_from_slice(&x.to_be_bytes()); +} + +fn f(x: f64, o: &mut Vec) { + o.extend_from_slice(&x.to_bits().to_be_bytes()); +} + +fn text(s: &str, o: &mut Vec) { + u(s.len() as u64, o); + o.extend_from_slice(s.as_bytes()); +} + +fn set(s: InstructionSet) -> u8 { + match s { + InstructionSet::GateBased => 1, + InstructionSet::Measurement => 2, + } +} + +fn isa(i: &IsaCounts, o: &mut Vec) { + for x in [i.tiles, i.min_steps, i.t_states, i.t_per_rotation] { + u(x, o); + } + o.push(i.synthesizes as u8); +} + +fn qubit(q: &PhysicalQubit, o: &mut Vec) { + text(&q.name, o); + o.push(set(q.instruction_set)); + for x in [q.gate_ns, q.meas_ns, q.clifford_error, q.t_error] { + f(x, o); + } +} + +fn code(q: &QecModel, o: &mut Vec) { + text(&q.name, o); + o.push(set(q.instruction_set)); + f(q.prefactor, o); + f(q.threshold, o); + u(q.qubits.0, o); + u(q.qubits.1, o); + u(q.qubits.2 as u64, o); + f(q.cycle.0, o); + f(q.cycle.1, o); +} + +fn kind_tag(k: UnitKind) -> u8 { + match k { + UnitKind::SpaceEfficient => 1, + UnitKind::ReedMuller => 2, + } +} + +fn result_hash(e: &Estimate) -> [u8; 32] { + let mut o = DOMAIN_RESULT.to_vec(); + text(&e.qubit, &mut o); + text(&e.code, &mut o); + isa(&e.isa, &mut o); + for x in [e.budget.logical, e.budget.distillation, e.budget.synthesis] { + f(x, &mut o); + } + u(e.distance as u64, &mut o); + u(e.tile_qubits, &mut o); + f(e.step_ns, &mut o); + u(e.steps, &mut o); + f(e.runtime_ns, &mut o); + f(e.logical_error, &mut o); + f(e.t_target, &mut o); + match &e.factory { + Some(fa) => { + o.push(1); + u(fa.rounds.len() as u64, &mut o); + for r in &fa.rounds { + o.push(kind_tag(r.kind)); + u(r.distance as u64, &mut o); + u(r.copies, &mut o); + u(r.unit_qubits, &mut o); + f(r.duration_ns, &mut o); + f(r.failure, &mut o); + f(r.output_error, &mut o); + } + u(fa.qubits, &mut o); + f(fa.duration_ns, &mut o); + f(fa.output_error, &mut o); + u(fa.outputs, &mut o); + } + None => o.push(0), + } + for x in [e.factories, e.algorithm_qubits, e.factory_qubits, e.physical_qubits] { + u(x, &mut o); + } + *blake3::hash(&o).as_bytes() +} + +fn hx(b: &[u8]) -> String { + b.iter().map(|x| format!("{x:02x}")).collect() +} + +fn unhex(s: &str) -> Option<[u8; N]> { + if s.len() != 2 * N || !s.bytes().all(|b| b.is_ascii_hexdigit()) { + return None; + } + let mut out = [0u8; N]; + for (i, c) in s.as_bytes().chunks(2).enumerate() { + out[i] = u8::from_str_radix(std::str::from_utf8(c).ok()?, 16).ok()?; + } + Some(out) +} + +fn bits(x: f64) -> String { + format!("{:016x}", x.to_bits()) +} + +fn unbits(v: &Value) -> Option { + u64::from_str_radix(v.as_str()?, 16).ok().map(f64::from_bits) +} + +/// A signed record of one resource estimate and the energy spent on it. +#[derive(Clone, Debug, PartialEq)] +pub struct ResourceReceipt { + pub job_hash: [u8; 32], + pub estimate: Estimate, + pub result_hash: [u8; 32], + /// Merkle root over (job, estimate) hashes. + pub merkle_root: [u8; 32], + /// Measured microjoules the sink spent (`0` = unmetered). Attested by the + /// signature, not verifiable by anyone else. + pub joules_micro: u64, + pub parent_receipt_hash: Option<[u8; 32]>, + pub signer_pubkey: [u8; 32], + pub signer_id: String, + pub sig: [u8; 64], +} + +impl ResourceReceipt { + fn leaves(&self) -> Vec<[u8; 32]> { + vec![self.job_hash, self.result_hash] + } + + fn payload_with(&self, label: Option<&crate::quantum_energy::EnergyClass>) -> Vec { + let mut o = crate::quantum_energy::labelled_domain(DOMAIN_RECEIPT, label); + o.extend_from_slice(&self.job_hash); + o.extend_from_slice(&self.result_hash); + o.extend_from_slice(&self.merkle_root); + o.extend_from_slice(&self.joules_micro.to_be_bytes()); + crate::quantum_energy::write_label(label, &mut o); + match self.parent_receipt_hash { + Some(h) => { + o.push(1); + o.extend_from_slice(&h); + } + None => o.push(0), + } + o.extend_from_slice(&self.signer_pubkey); + o.extend_from_slice(self.signer_id.as_bytes()); + o + } + + /// Run `job` under `meter` and seal what it computed. `meter` is handed a + /// closure that performs the run; it must call it once, and returns the + /// microjoules it measured (`|run| { run(); 0 }` for an unmetered seal). + pub fn seal( + signer: &SigningKey, + signer_id: impl Into, + job: &ResourceJob, + parent_receipt_hash: Option<[u8; 32]>, + meter: impl FnOnce(&mut dyn FnMut()) -> u64, + ) -> Result { + let mut out = None; + let joules_micro = meter(&mut || out = Some(job.run())); + let estimate = out.ok_or(EstimateError::NotRun)??; + let mut r = ResourceReceipt { + job_hash: job.job_hash(), + result_hash: result_hash(&estimate), + estimate, + merkle_root: [0; 32], + joules_micro, + parent_receipt_hash, + signer_pubkey: signer.verifying_key().to_bytes(), + signer_id: signer_id.into(), + sig: [0; 64], + }; + r.merkle_root = merkle_root(&r.leaves()); + r.sig = signer.sign(&r.payload_with(None)).to_bytes(); + Ok(r) + } + + fn body_verifies(&self) -> bool { + result_hash(&self.estimate) == self.result_hash && merkle_root(&self.leaves()) == self.merkle_root + } + + /// The receipt is internally consistent and its signature is valid. Does + /// not re-run the job: see [`Self::verify_computation`]. + pub fn verify(&self) -> bool { + if !self.body_verifies() { + return false; + } + let Ok(k) = VerifyingKey::from_bytes(&self.signer_pubkey) else { + return false; + }; + k.verify(&self.payload_with(None), &Signature::from_bytes(&self.sig)).is_ok() + } + + /// The re-run half of [`Self::verify_computation`]. + pub(crate) fn computation_matches(&self, job: &ResourceJob) -> bool { + job.job_hash() == self.job_hash && job.run().is_ok_and(|e| result_hash(&e) == self.result_hash) + } + + /// The strong check: `job` is the job this receipt names, re-running it + /// here reproduces the sealed estimate bit for bit, and the receipt + /// verifies. + pub fn verify_computation(&self, job: &ResourceJob) -> bool { + self.computation_matches(job) && self.verify() + } + + /// Stable id, for chaining a derived receipt via `parent_receipt_hash`. + pub fn receipt_hash(&self) -> [u8; 32] { + let mut h = blake3::Hasher::new(); + h.update(DOMAIN_RECEIPT_ID); + h.update(&self.payload_with(None)); + h.update(&self.sig); + *h.finalize().as_bytes() + } + + /// Canonical JSON, `kind` first and then keys sorted. Every floating-point + /// figure appears as its exact bits (`*_bits`), which are what a verifier + /// reads back; the headline runtime also appears as a readable number. + pub fn to_json(&self) -> String { + let e = &self.estimate; + let factory = e.factory.as_ref().map(|fa| { + json!({ + "duration_ns_bits": bits(fa.duration_ns), + "output_error_bits": bits(fa.output_error), + "outputs": fa.outputs, + "qubits": fa.qubits, + "rounds": fa.rounds.iter().map(|r| json!({ + "copies": r.copies, + "distance": r.distance, + "duration_ns_bits": bits(r.duration_ns), + "failure_bits": bits(r.failure), + "kind": match r.kind { UnitKind::SpaceEfficient => "space-efficient", UnitKind::ReedMuller => "reed-muller" }, + "output_error_bits": bits(r.output_error), + "unit_qubits": r.unit_qubits, + })).collect::>(), + }) + }); + let v = json!({ + "kind": "quantum-resource", + "estimate": { + "algorithm_qubits": e.algorithm_qubits, + "budget_bits": [bits(e.budget.logical), bits(e.budget.distillation), bits(e.budget.synthesis)], + "code": e.code, + "distance": e.distance, + "factories": e.factories, + "factory": factory, + "factory_qubits": e.factory_qubits, + "isa": { + "min_steps": e.isa.min_steps, + "synthesizes": e.isa.synthesizes, + "t_per_rotation": e.isa.t_per_rotation, + "t_states": e.isa.t_states, + "tiles": e.isa.tiles, + }, + "logical_error_bits": bits(e.logical_error), + "physical_qubits": e.physical_qubits, + "qubit": e.qubit, + "runtime_ns_bits": bits(e.runtime_ns), + "runtime_s": e.runtime_ns / 1e9, + "step_ns_bits": bits(e.step_ns), + "steps": e.steps, + "t_target_bits": bits(e.t_target), + "tile_qubits": e.tile_qubits, + }, + "job_hash": hx(&self.job_hash), + "joules_micro": self.joules_micro, + "parent_receipt_hash": self.parent_receipt_hash.map(|h| hx(&h)), + "receipt_hash": hx(&self.receipt_hash()), + "result_hash": hx(&self.result_hash), + "root_hash": hx(&self.merkle_root), + "sig": hx(&self.sig), + "signer_id": self.signer_id, + "signer_pubkey": hx(&self.signer_pubkey), + }); + // The crate's serde_json keeps insertion order, and keys are written in + // sorted order after `kind`, so this is canonical. + v.to_string() + } + + pub fn from_json(s: &str) -> Option { + let v: Value = serde_json::from_str(s).ok()?; + if v.get("kind")?.as_str()? != "quantum-resource" { + return None; + } + let e = v.get("estimate")?; + let n = |o: &Value, k: &str| o.get(k).and_then(Value::as_u64); + let b = |o: &Value, k: &str| o.get(k).and_then(unbits); + let i = e.get("isa")?; + let budget = e.get("budget_bits")?.as_array()?; + let factory = match e.get("factory")? { + Value::Null => None, + fa => { + let mut rounds = Vec::new(); + for r in fa.get("rounds")?.as_array()? { + rounds.push(Round { + kind: match r.get("kind")?.as_str()? { + "space-efficient" => UnitKind::SpaceEfficient, + "reed-muller" => UnitKind::ReedMuller, + _ => return None, + }, + distance: u32::try_from(n(r, "distance")?).ok()?, + copies: n(r, "copies")?, + unit_qubits: n(r, "unit_qubits")?, + duration_ns: b(r, "duration_ns_bits")?, + failure: b(r, "failure_bits")?, + output_error: b(r, "output_error_bits")?, + }); + } + Some(Factory { + rounds, + qubits: n(fa, "qubits")?, + duration_ns: b(fa, "duration_ns_bits")?, + output_error: b(fa, "output_error_bits")?, + outputs: n(fa, "outputs")?, + }) + } + }; + let estimate = Estimate { + qubit: e.get("qubit")?.as_str()?.into(), + code: e.get("code")?.as_str()?.into(), + isa: IsaCounts { + tiles: n(i, "tiles")?, + min_steps: n(i, "min_steps")?, + t_states: n(i, "t_states")?, + t_per_rotation: n(i, "t_per_rotation")?, + synthesizes: i.get("synthesizes")?.as_bool()?, + }, + budget: Budget { + logical: unbits(budget.first()?)?, + distillation: unbits(budget.get(1)?)?, + synthesis: unbits(budget.get(2)?)?, + }, + distance: u32::try_from(n(e, "distance")?).ok()?, + tile_qubits: n(e, "tile_qubits")?, + step_ns: b(e, "step_ns_bits")?, + steps: n(e, "steps")?, + runtime_ns: b(e, "runtime_ns_bits")?, + logical_error: b(e, "logical_error_bits")?, + t_target: b(e, "t_target_bits")?, + factory, + factories: n(e, "factories")?, + algorithm_qubits: n(e, "algorithm_qubits")?, + factory_qubits: n(e, "factory_qubits")?, + physical_qubits: n(e, "physical_qubits")?, + }; + let s = |k: &str| v.get(k).and_then(Value::as_str); + let parent = match v.get("parent_receipt_hash") { + Some(Value::String(h)) => Some(unhex::<32>(h)?), + _ => None, + }; + Some(ResourceReceipt { + job_hash: unhex(s("job_hash")?)?, + estimate, + result_hash: unhex(s("result_hash")?)?, + merkle_root: unhex(s("root_hash")?)?, + joules_micro: v.get("joules_micro")?.as_u64()?, + parent_receipt_hash: parent, + signer_pubkey: unhex(s("signer_pubkey")?)?, + signer_id: s("signer_id")?.to_string(), + sig: unhex(s("sig")?)?, + }) + } +} + +impl Default for ResourceJob { + /// 2048-bit factoring on the `(ns, 10⁻⁴)` gate-based qubit, `ε = 1/3`. + fn default() -> Self { + ResourceJob { + program: Program::Logical(LogicalCounts { + qubits: 12_581, + measurements: 1_080_000_000, + rotations: 12, + rotation_depth: 12, + t_gates: 12, + ccz: 3_730_000_000, + }), + qubit: PhysicalQubit::ns_e4(), + code: None, + config: EstimateConfig { error_budget: 1.0 / 3.0, slowdown: 1, limits: SearchLimits::default() }, + } + } +} + +crate::quantum_energy::labellable!(ResourceReceipt, DOMAIN_RECEIPT_ID); + +#[cfg(test)] +mod tests { + use super::*; + + fn key(b: u8) -> SigningKey { + SigningKey::from_bytes(&[b; 32]) + } + + fn unmetered(run: &mut dyn FnMut()) -> u64 { + run(); + 0 + } + + #[test] + fn a_sealed_estimate_verifies_and_reproduces() { + let j = ResourceJob::default(); + let r = ResourceReceipt::seal(&key(1), "sink-a", &j, None, |run| { + run(); + 1_700 + }) + .unwrap(); + assert!(r.verify()); + assert!(r.verify_computation(&j)); + assert_eq!(r.estimate, j.run().unwrap()); + assert_eq!((r.estimate.distance, r.estimate.factories), (13, 18)); + assert_eq!(r.joules_micro, 1_700); + } + + #[test] + fn the_json_round_trips_to_the_same_bytes() { + // With a factory, and a Clifford-only program without one. + let clifford = ResourceJob { + program: Program::Logical(LogicalCounts { qubits: 50, measurements: 1_000, ..Default::default() }), + ..ResourceJob::default() + }; + for job in [ResourceJob::default(), clifford] { + let r = ResourceReceipt::seal(&key(2), "sink \"b\"", &job, Some([7; 32]), unmetered).unwrap(); + let back = ResourceReceipt::from_json(&r.to_json()).unwrap(); + assert_eq!(back, r); + assert_eq!(back.to_json(), r.to_json()); + assert!(back.verify()); + assert!(back.verify_computation(&job)); + } + } + + #[test] + fn tampering_is_caught() { + let j = ResourceJob::default(); + let r = ResourceReceipt::seal(&key(3), "sink", &j, None, unmetered).unwrap(); + let mut a = r.clone(); + a.joules_micro = 1; + assert!(!a.verify(), "the energy is signed"); + let mut b = r.clone(); + b.estimate.physical_qubits -= 1; + assert!(!b.verify(), "the headline is hashed"); + let mut c = r.clone(); + c.estimate.factory.as_mut().unwrap().rounds[0].copies = 15; + assert!(!c.verify(), "every factory round is hashed"); + // One ulp on the physical error rate is a different job. + let mut other = j.clone(); + other.qubit.clifford_error = f64::from_bits(other.qubit.clifford_error.to_bits() + 1); + assert!(!r.verify_computation(&other)); + let mut fixed = j.clone(); + fixed.code = Some(QecModel::surface_gate()); + assert!(!r.verify_computation(&fixed), "choosing the code is part of the job"); + let mut slower = j.clone(); + slower.config.slowdown = 2; + assert!(!r.verify_computation(&slower)); + } + + #[test] + fn a_meter_that_never_runs_the_job_cannot_seal() { + assert_eq!( + ResourceReceipt::seal(&key(4), "sink", &ResourceJob::default(), None, |_run| 99), + Err(EstimateError::NotRun) + ); + } + + #[test] + fn an_energy_class_label_seals_over_it() { + use crate::quantum_energy::{EnergyClass, Labelled}; + let r = ResourceReceipt::seal(&key(5), "sink", &ResourceJob::default(), None, |run| { + run(); + 1_000 + }) + .unwrap(); + let labelled = Labelled::seal(r, EnergyClass::ModelBased, &key(5)).unwrap(); + assert!(labelled.verify()); + assert!(!labelled.receipt.verify(), "the inner signature now covers the label"); + } +} diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs index 5b50dc8c1..99902a644 100644 --- a/wai/wai-quantum/src/repro.rs +++ b/wai/wai-quantum/src/repro.rs @@ -7,146 +7,51 @@ //! fdlibm's algorithm, as musl carries it, using only `+ − × ÷`, which Rust //! never fuses. -/// `sin` and `cos` of `x`, bit-identical on every IEEE-754 platform. The -/// algorithm is fdlibm's (as musl carries it): Cody–Waite reduction by π/2 in -/// three parts, then minimax polynomials on [−π/4, π/4]. Only `+ − × ÷` are -/// used, and Rust never fuses them, so the bits do not depend on the target. -/// Accurate to within an ulp of the true value. -/// -/// For `|x| > 2²⁰` (and for non-finite `x`) both are NaN: the reduction used -/// here is exact only below that, and a NaN is better than a value that -/// depends on the platform. -pub fn sin_cos(x: f64) -> (f64, f64) { - if !(x.abs() <= 1_048_576.0) { - return (f64::NAN, f64::NAN); - } - let ix = (x.to_bits() >> 32) as u32 & 0x7fff_ffff; - if ix <= 0x3fe9_21fb { - // |x| ≤ π/4 - if ix < 0x3e40_0000 { - // |x| < 2⁻²⁷: sin x = x and cos x = 1 to the last bit. - return (x, 1.0); - } - return (k_sin(x, 0.0, false), k_cos(x, 0.0)); - } - let (n, y0, y1) = rem_pio2(x, ix); - let (s, c) = (k_sin(y0, y1, true), k_cos(y0, y1)); - match n & 3 { - 0 => (s, c), - 1 => (c, -s), - 2 => (-s, -c), - _ => (-c, s), - } -} - -const S1: f64 = -1.666_666_666_666_663_243_48e-1; -const S2: f64 = 8.333_333_333_322_489_461_24e-3; -const S3: f64 = -1.984_126_982_985_794_931_34e-4; -const S4: f64 = 2.755_731_370_707_006_767_89e-6; -const S5: f64 = -2.505_076_025_340_686_341_95e-8; -const S6: f64 = 1.589_690_995_211_550_102_21e-10; - -fn k_sin(x: f64, y: f64, tail: bool) -> f64 { - let z = x * x; - let w = z * z; - let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6); - let v = z * x; - if !tail { - x + v * (S1 + z * r) - } else { - x - ((z * (0.5 * y - v * r) - y) - v * S1) - } -} - -const C1: f64 = 4.166_666_666_666_660_190_37e-2; -const C2: f64 = -1.388_888_888_887_410_957_49e-3; -const C3: f64 = 2.480_158_728_947_672_941_78e-5; -const C4: f64 = -2.755_731_435_139_066_330_35e-7; -const C5: f64 = 2.087_572_321_298_174_827_90e-9; -const C6: f64 = -1.135_964_755_778_819_482_65e-11; - -fn k_cos(x: f64, y: f64) -> f64 { - let z = x * x; - let w = z * z; - let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6)); - let hz = 0.5 * z; - let w = 1.0 - hz; - w + (((1.0 - w) - hz) + (z * r - x * y)) -} - -const TOINT: f64 = 1.5 / f64::EPSILON; -const PIO4: f64 = 7.853_981_633_974_482_789_99e-1; -const INVPIO2: f64 = 6.366_197_723_675_813_824_33e-1; -const PIO2_1: f64 = 1.570_796_326_734_125_614_17; -const PIO2_1T: f64 = 6.077_100_506_506_192_249_32e-11; -const PIO2_2: f64 = 6.077_100_506_303_965_976_60e-11; -const PIO2_2T: f64 = 2.022_266_248_795_950_631_54e-21; -const PIO2_3: f64 = 2.022_266_248_711_166_455_80e-21; -const PIO2_3T: f64 = 8.478_427_660_368_899_569_97e-32; - -/// `x = n·π/2 + (y0 + y1)` with `|y0 + y1| ≤ π/4`, for `|x| ≤ 2²⁰`. -fn rem_pio2(x: f64, ix: u32) -> (i64, f64, f64) { - let mut f = (x * INVPIO2 + TOINT) - TOINT; - let mut r = x - f * PIO2_1; - let mut w = f * PIO2_1T; - // Matters only under directed rounding; kept so the port is exact. - if r - w < -PIO4 { - f -= 1.0; - r = x - f * PIO2_1; - w = f * PIO2_1T; - } else if r - w > PIO4 { - f += 1.0; - r = x - f * PIO2_1; - w = f * PIO2_1T; - } - let mut y0 = r - w; - let ex = (ix >> 20) as i32; - let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32; - if ex - ey > 16 { - let t = r; - w = f * PIO2_2; - r = t - w; - w = f * PIO2_2T - ((t - r) - w); - y0 = r - w; - let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32; - if ex - ey > 49 { - let t = r; - w = f * PIO2_3; - r = t - w; - w = f * PIO2_3T - ((t - r) - w); - y0 = r - w; - } - } - let y1 = (r - y0) - w; - (f as i64, y0, y1) -} +#[cfg(any( + feature = "quantum_spd", + feature = "quantum_pauli", + feature = "quantum_mps", + feature = "quantum_vml", + feature = "quantum_phasor", + feature = "quantum_kernel" +))] +mod trig; +#[cfg(any( + feature = "quantum_spd", + feature = "quantum_pauli", + feature = "quantum_mps", + feature = "quantum_vml", + feature = "quantum_phasor", + feature = "quantum_kernel" +))] +pub use trig::sin_cos; // --------------------------------------------------------------------------- // Natural logarithm // --------------------------------------------------------------------------- -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG1: f64 = 6.666_666_666_666_735_130e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG2: f64 = 3.999_999_999_940_941_908e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG3: f64 = 2.857_142_874_366_239_149e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG4: f64 = 2.222_219_843_214_978_396e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG5: f64 = 1.818_357_216_161_805_012e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG6: f64 = 1.531_383_769_920_937_332e-1; -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] const LG7: f64 = 1.479_819_860_511_658_591e-1; /// The natural logarithm, within an ulp of the true value: `−∞` at zero, NaN /// below it. -#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] pub fn ln(x: f64) -> f64 { let mut bits = x.to_bits(); let mut hx = (bits >> 32) as u32; @@ -357,7 +262,7 @@ mod tests { }) } - #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] + #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource"))] #[test] fn ln_is_within_an_ulp_of_the_platform() { for x in inputs().filter(|&x| x > 0.0) { diff --git a/wai/wai-quantum/src/repro/trig.rs b/wai/wai-quantum/src/repro/trig.rs new file mode 100644 index 000000000..94a8d2391 --- /dev/null +++ b/wai/wai-quantum/src/repro/trig.rs @@ -0,0 +1,116 @@ +//! `sin` and `cos`, for the layers that rotate: fdlibm's algorithm as musl +//! carries it, using only `+ − × ÷`. + +/// `sin` and `cos` of `x`, bit-identical on every IEEE-754 platform. The +/// algorithm is fdlibm's (as musl carries it): Cody–Waite reduction by π/2 in +/// three parts, then minimax polynomials on [−π/4, π/4]. Only `+ − × ÷` are +/// used, and Rust never fuses them, so the bits do not depend on the target. +/// Accurate to within an ulp of the true value. +/// +/// For `|x| > 2²⁰` (and for non-finite `x`) both are NaN: the reduction used +/// here is exact only below that, and a NaN is better than a value that +/// depends on the platform. +pub fn sin_cos(x: f64) -> (f64, f64) { + if !(x.abs() <= 1_048_576.0) { + return (f64::NAN, f64::NAN); + } + let ix = (x.to_bits() >> 32) as u32 & 0x7fff_ffff; + if ix <= 0x3fe9_21fb { + // |x| ≤ π/4 + if ix < 0x3e40_0000 { + // |x| < 2⁻²⁷: sin x = x and cos x = 1 to the last bit. + return (x, 1.0); + } + return (k_sin(x, 0.0, false), k_cos(x, 0.0)); + } + let (n, y0, y1) = rem_pio2(x, ix); + let (s, c) = (k_sin(y0, y1, true), k_cos(y0, y1)); + match n & 3 { + 0 => (s, c), + 1 => (c, -s), + 2 => (-s, -c), + _ => (-c, s), + } +} + +const S1: f64 = -1.666_666_666_666_663_243_48e-1; +const S2: f64 = 8.333_333_333_322_489_461_24e-3; +const S3: f64 = -1.984_126_982_985_794_931_34e-4; +const S4: f64 = 2.755_731_370_707_006_767_89e-6; +const S5: f64 = -2.505_076_025_340_686_341_95e-8; +const S6: f64 = 1.589_690_995_211_550_102_21e-10; + +fn k_sin(x: f64, y: f64, tail: bool) -> f64 { + let z = x * x; + let w = z * z; + let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6); + let v = z * x; + if !tail { + x + v * (S1 + z * r) + } else { + x - ((z * (0.5 * y - v * r) - y) - v * S1) + } +} + +const C1: f64 = 4.166_666_666_666_660_190_37e-2; +const C2: f64 = -1.388_888_888_887_410_957_49e-3; +const C3: f64 = 2.480_158_728_947_672_941_78e-5; +const C4: f64 = -2.755_731_435_139_066_330_35e-7; +const C5: f64 = 2.087_572_321_298_174_827_90e-9; +const C6: f64 = -1.135_964_755_778_819_482_65e-11; + +fn k_cos(x: f64, y: f64) -> f64 { + let z = x * x; + let w = z * z; + let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6)); + let hz = 0.5 * z; + let w = 1.0 - hz; + w + (((1.0 - w) - hz) + (z * r - x * y)) +} + +const TOINT: f64 = 1.5 / f64::EPSILON; +const PIO4: f64 = 7.853_981_633_974_482_789_99e-1; +const INVPIO2: f64 = 6.366_197_723_675_813_824_33e-1; +const PIO2_1: f64 = 1.570_796_326_734_125_614_17; +const PIO2_1T: f64 = 6.077_100_506_506_192_249_32e-11; +const PIO2_2: f64 = 6.077_100_506_303_965_976_60e-11; +const PIO2_2T: f64 = 2.022_266_248_795_950_631_54e-21; +const PIO2_3: f64 = 2.022_266_248_711_166_455_80e-21; +const PIO2_3T: f64 = 8.478_427_660_368_899_569_97e-32; + +/// `x = n·π/2 + (y0 + y1)` with `|y0 + y1| ≤ π/4`, for `|x| ≤ 2²⁰`. +fn rem_pio2(x: f64, ix: u32) -> (i64, f64, f64) { + let mut f = (x * INVPIO2 + TOINT) - TOINT; + let mut r = x - f * PIO2_1; + let mut w = f * PIO2_1T; + // Matters only under directed rounding; kept so the port is exact. + if r - w < -PIO4 { + f -= 1.0; + r = x - f * PIO2_1; + w = f * PIO2_1T; + } else if r - w > PIO4 { + f += 1.0; + r = x - f * PIO2_1; + w = f * PIO2_1T; + } + let mut y0 = r - w; + let ex = (ix >> 20) as i32; + let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32; + if ex - ey > 16 { + let t = r; + w = f * PIO2_2; + r = t - w; + w = f * PIO2_2T - ((t - r) - w); + y0 = r - w; + let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32; + if ex - ey > 49 { + let t = r; + w = f * PIO2_3; + r = t - w; + w = f * PIO2_3T - ((t - r) - w); + y0 = r - w; + } + } + let y1 = (r - y0) - w; + (f as i64, y0, y1) +} diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs index 20d2f74c4..5edf8663a 100644 --- a/wai/wai-quantum/tests/bit_identity.rs +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -14,7 +14,8 @@ feature = "quantum_mps", feature = "quantum_pauli", feature = "quantum_spd", - feature = "quantum_bptn" + feature = "quantum_bptn", + feature = "quantum_resource" ))] use wai_quantum::*; @@ -121,3 +122,28 @@ fn tensor_network_evolution() { v.push(t.discarded()); assert_eq!(digest(&v), 0x0e51_33ff_7987_aa26); } + +/// Every figure of the factoring and chemistry estimates on all six reference +/// qubits: distances, factory rounds, qubit counts, times and error rates. +/// The rotation T-count goes through the in-crate logarithm. +#[test] +fn resource_estimates() { + use quantum_resource::*; + let programs = [ + (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), + (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), + ]; + let mut v = Vec::new(); + for (c, eps) in programs { + for q in PhysicalQubit::presets() { + let e = estimate_best(&c, &q, &EstimateConfig::new(eps)).unwrap(); + v.extend([e.distance as f64, e.steps as f64, e.runtime_ns, e.logical_error, e.t_target]); + v.extend([e.factories as f64, e.physical_qubits as f64, e.isa.t_per_rotation as f64]); + for r in &e.factory.unwrap().rounds { + v.extend([r.distance as f64, r.copies as f64, r.duration_ns, r.failure, r.output_error]); + } + } + } + assert_eq!(digest(&v), 0xce9d_43b2_d498_f379); +} + diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index f575aca24..283b86ca3 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -3900,7 +3900,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.23" +version = "0.3.24" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index f7496f134..d83397713 100644 --- a/wai/wai-rs/Cargo.toml +++ b/wai/wai-rs/Cargo.toml @@ -341,6 +341,10 @@ quantum_spd = ["quantum", "wai-quantum/quantum_spd"] quantum_spd_receipt = ["quantum_spd", "wai-quantum/quantum_spd_receipt"] # Tensor-network states on the hardware graph (simple update, belief propagation). quantum_bptn = ["quantum_spd", "wai-quantum/quantum_bptn"] +# Fault-tolerant resource estimation (layout, distance, T factories, frontier), +# and its signed receipts. +quantum_resource = ["dep:wai-quantum", "wai-quantum/quantum_resource"] +quantum_resource_receipt = ["quantum_resource", "wai-quantum/quantum_resource_receipt"] quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"] # Variational quantum machine learning (extensions/quantum-sim § VML, # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot @@ -468,6 +472,7 @@ required-features = [ "quantum_qudit", "quantum_spd_receipt", "quantum_bptn", + "quantum_resource_receipt", ] [[bin]] diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index 1036a78b9..059a29314 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -36,6 +36,17 @@ //! wai-quantum 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 estimate factoring # 2048-bit factoring, all six qubits +//! wai-quantum estimate chemistry --qubit ns-e4 --frontier +//! wai-quantum estimate --qubits 100 --ccz 1e6 --budget 0.01 +//! wai-quantum estimate factoring --qubit ns-e4 --receipt | wai-quantum verify +//! ``` +//! //! 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 @@ -51,6 +62,8 @@ use wai::quantum_info::{chsh, DensityMatrix}; use wai::quantum_qudit::QuditState; use wai::quantum_energy::{EnergyClass, MaybeLabelled}; use wai::quantum_receipt::QuantumReceipt; +use wai::quantum_resource::{self as qr, EstimateConfig, LogicalCounts, PhysicalQubit, QecModel}; +use wai::quantum_resource_receipt::{Program, ResourceJob, ResourceReceipt}; use wai::quantum_spd_receipt::{SpdJob, SpdReceipt}; use wai::quantum_source::{holevo_bound, schumacher_limit, von_neumann_entropy}; use wai::quantum_toolchain as qt; @@ -128,6 +141,11 @@ fn usage() -> ! { UTILITY SCALE:\n \ wai-quantum utility [--theta T] [--steps N] [--threshold D] [-q QUBIT] [--json] [--receipt]\n \ wai-quantum utility --method bptn [--chi X] [--theta T] [--steps N] [-q QUBIT] [--json]\n\n\ + FAULT-TOLERANT RESOURCES:\n \ + wai-quantum estimate [factoring|chemistry] [--qubit NAME|all] [--code surface|surface-meas|floquet]\n \ + [--budget E] [--slowdown N] [--frontier] [--json] [--receipt]\n \ + wai-quantum 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.\n\n\ Channels: identity, bitflip, dephasing, depolarizing, damping.\n\ Run `wai-quantum list` for the algorithm catalog." ); @@ -185,6 +203,126 @@ fn describe(c: &Circuit) -> (String, String, u64, String) { ) } +/// 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. +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 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, + ) +} + fn main() { let argv: Vec = std::env::args().collect(); if argv.len() < 2 { @@ -270,6 +408,25 @@ fn main() { // carries its acquisition class inside the signature // (`quantum_energy::Labelled`, as `wai_quantum_meter` emits). let src = read_stdin(); + 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 {} uJ", 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(); @@ -388,6 +545,112 @@ fn main() { } } + // 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)")), + }); + 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 }; + // The same deterministic demo key as `receipt`: the estimate is + // bit-identical everywhere, so is its receipt. + let key = SigningKey::from_bytes(&[42u8; 32]); + let r = ResourceReceipt::seal(&key, "did:key:wai-quantum-cli", &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} {}", "slowdown", "d", "factories", "physical", "runtime"); + 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} {:<13} {:>3} {:>9} {:>8} {:>14} {}", + "qubit", "code", "d", "factories", "share", "physical", "runtime" + ); + for (q, e) in &estimates { + match e { + Ok(e) => println!( + "{:<12} {:<13} {:>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), + } + } + } + "entropy" => { let c = circuit_from(args, algo); let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}"))); diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index d1e8661f7..ff8928a8e 100644 --- a/wai/wai-rs/src/lib.rs +++ b/wai/wai-rs/src/lib.rs @@ -590,6 +590,12 @@ pub use wai_quantum::quantum_spd_receipt; #[cfg(feature = "quantum_bptn")] pub use wai_quantum::quantum_bptn; + +#[cfg(feature = "quantum_resource")] +pub use wai_quantum::quantum_resource; + +#[cfg(feature = "quantum_resource_receipt")] +pub use wai_quantum::quantum_resource_receipt; #[cfg(feature = "quantum_qir")] pub use wai_quantum::quantum_qir; diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 4720f02a9..0da1aeda9 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -734,7 +734,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.23" +version = "0.3.24" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157)