wai-quantum 0.4.0

//! Architecture-level ESTIMATION — `wai.quantum.arch`.
//!
//! Since 2023 the projected cost of a fault-tolerant computation has fallen by
//! an order of magnitude without a better qubit. The savings come from the
//! architecture around the code:
//!
//! - **Yoked storage** (arXiv:2312.04522). Idle logical qubits sit in
//!   surface-code patches bound together by row (and column) parity checks.
//!   This takes about a third of the qubits of plain patches at the same
//!   error.
//! - **Magic-state cultivation** (arXiv:2409.17595). A T state is grown inside
//!   one patch and kept only if every check passes, instead of being distilled
//!   from fifteen.
//! - **Reaction-limited scheduling.** A step runs at the pace of the slowest of
//!   three things: lattice surgery, the magic-state supply, or the classical
//!   control loop.
//!
//! This module prices all three and assembles them into a whole-machine
//! estimate. As in the planar estimator (`quantum_resource`), every number is a
//! deterministic function of its inputs: `+ − × ÷` and in-crate logarithms,
//! with the same bits on every machine.
//!
//! # Storage
//!
//! A patch of distance `d` covers `2(d+1)²` physical qubits; the margin leaves
//! room for lattice surgery. Storage errors follow path-counting fits. Each
//! fit gives the logical error per patch per round:
//!
//! ```text
//! p = r^(a−1) · n^(b−1) / (C · Λ^d)
//! ```
//!
//! Here `r` is the number of rounds between yoke checks and `n` the number of
//! patches per yoked block. The exponents are `(a, b) = (1, 1)` for plain
//! patches, `(2, 2)` for row yokes and `(4, 2)` for grid yokes. This is
//! multiplied by `n / k` to give the error per *logical* qubit per round.
//!
//! - **Cold storage.** Blocks share one hallway that walks through them to
//!   measure the yokes. A row-yoked block of `g` groups is checked every
//!   `(8g + 2)·d` rounds. A `w × w` grid block is checked every `(25w + 4)·d`
//!   rounds and stores `w² − 4w + 2` logical qubits in `(w+1)²` patches.
//! - **Hot storage.** Every patch faces an access hallway. Row yokes are
//!   checked every `50d` rounds.
//!
//! The search for the cheapest layout per logical qubit is held to the source's
//! own footprint tool. There are three fit sets ([`StorageModel::si1000`],
//! [`StorageModel::uniform`], [`StorageModel::uniform_refit`]). For each, the
//! layout is checked at 61 targets from 10⁻¹⁸ to 10⁻³, for plain, row-yoked
//! and grid-yoked storage, cold and hot. Every distance, block shape and
//! qubit count is identical.
//!
//! # Magic states
//!
//! [`cultivate`] reads the cheapest T state meeting an error target off the
//! cultivation curves. Those curves come from the published end-to-end
//! simulations: 10⁹–10¹² shots per point, released under Apache-2.0. A curve
//! gives the expected qubit·rounds per kept state (retries included) at each
//! postselection cutoff, for injection distance 3 or 5 and noise 5·10⁻⁴,
//! 10⁻³ or 2·10⁻³.
//!
//! The simulations swap the T gate for an S gate, so they can run on a
//! stabilizer simulator. As the source does, the error is then doubled to
//! stand in for T.
//!
//! [`CczFactory`] turns eight T states into one CCZ state, with output error
//! `28·p_T²` (arXiv:1212.5069). The factory is a 3 × 4 block of patches running
//! six layers of temporally encoded lattice surgery.
//!
//! # Time
//!
//! [`Timing`] sets the clock. A step takes the longest of three times:
//! - the lattice-surgery period, `d` rounds;
//! - the CCZ supply period, rounds per CCZ divided by the number of factories;
//! - the reaction time of the control system.
//!
//! [`Timing::bound`] says which of the three binds.
//!
//! On that clock:
//! - an `n`-bit addition consumes `n − 1` CCZ states and then uncomputes, for
//!   `2(n−1)` steps;
//! - a lookup over `w` address bits takes `2ʷ − 1` steps;
//! - a phaseup takes half a lookup.
//!
//! # The 2025 factoring estimate, re-run
//!
//! arXiv:2505.15917 estimates that a 2048-bit RSA modulus can be factored in
//! under a week on under a million noisy qubits. Its assumptions are:
//! - a square grid;
//! - a uniform error rate of 10⁻³;
//! - a 1 µs cycle;
//! - a 10 µs reaction time.
//!
//! [`FactoringParams`] carries its algorithm parameters. [`factoring_estimate`]
//! re-runs the whole estimate from them.
//!
//! **Reproduced.**
//! - **Physical qubits.** The total of 897,864 is reproduced exactly from the
//!   source's layout:
//!   - 1280 cold logical qubits at 430 qubits each;
//!   - 131 hot patches at distance 25;
//!   - a 7 × 18 compute region holding six factories.
//! - **Cold-storage density.** The 430 is the published fit evaluated. Grid
//!   blocks 16 patches wide at distance 11 give 429.03 qubits per logical
//!   qubit at 10⁻¹⁵ per round.
//! - **Logical counts.** The source's Toffoli count, `6.5·10⁹` per factoring,
//!   and its 1399 logical qubits follow from its operation tallies under its
//!   counting rule.
//!
//! **Where the source disagrees with itself.** Re-running it surfaced six
//! places. None of them moves the conclusion: under a million qubits and under
//! a week.
//! - **Runtime per shot.** The source's tallies, at its stated durations (2 ms
//!   per addition and per lookup, 1 ms per phaseup), give 10.62 h per shot.
//!   It prints 12.07 h, so its own estimate is conservative by 14%.
//! - **Hot storage.** The layout keeps 131 logical qubits hot, but the source
//!   gives three other counts.
//!   - Its peak formula `m + 3f + 2ℓ + len m` evaluates to 1432 at loop 4,
//!     which is 152 hot.
//!   - The text prints 1409 beside that formula.
//!   - Its table gives 1399, which is 119 hot.
//!
//!   At 152 hot the machine takes 926,256 physical qubits, still under a
//!   million.
//! - **Phaseup Toffolis.** The source's reference counting code has an exponent
//!   sign slip. It evaluates a phaseup with a 6-bit address to 0.25 Toffolis;
//!   the text says 8. Counting as the text states adds 1.7% to the total.
//! - **Loop 4.** The table lists 1.5 additions and 2.5 lookups per iteration.
//!   The reference code performs 2.5 additions and 1.5 lookups. Both take
//!   2 ms, so the runtime is unchanged.
//! - **Cultivation.** The source reads 30,000 qubit·rounds for a 10⁻⁷ T state
//!   off a figure. The data behind that figure give about 22,600. This makes
//!   the factory estimate conservative as well.
//! - **Block granularity.** 1280 logical qubits do not fill whole 194-qubit
//!   grid blocks. Seven whole blocks cost 582,624 physical qubits, not
//!   550,400.
//!
//! [`factoring_estimate`] applies every correction. It finds 958,480 physical
//! qubits:
//! - seven cold grid blocks, 582,624;
//! - 152 hot patches, 205,504;
//! - the compute region, 170,352.
//!
//! At the clock, with each operation rounded up to a whole millisecond, a shot
//! takes 10.58 h. A factoring takes 4.3 days, shots lost to logical errors
//! included.
//!
//! # Honest boundaries
//!
//! - **Fits are fits.** The storage formulas are path-counting extrapolations
//!   from simulations at one noise strength, 10⁻³. Their authors call them
//!   conservative. Outside that regime they are guesses.
//! - **Cultivation is tabulated, not modelled.** Only the three simulated noise
//!   strengths are supported; [`cultivate`] refuses any other. Between
//!   tabulated cutoffs the volume is interpolated log-linearly. Points backed
//!   by few failures (see [`CultivationPoint::failures`]) are uncertain by a
//!   factor of a few.
//! - **The factoring workload is the source's.** The operation tallies and the
//!   deviation rate come from it, and are reproduced, not re-derived. The
//!   residue-system search behind them is a conjecture of the source
//!   (its Assumption 1).
//! - **Lattice surgery is priced, not routed.** Operation times are lower
//!   bounds from the CCZ count and the layer count. There is no placement or
//!   routing.

use crate::repro::{exp, ln};

/// Physical qubits a distance-`d` patch covers, with room for lattice surgery:
/// `2(d+1)²`.
pub fn patch_qubits(d: u32) -> u64 {
    2 * u64::from(d + 1) * u64::from(d + 1)
}

/// A path-counting fit for the logical error of stored patches.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct StorageFit {
    /// The error suppression factor `Λ` per unit of distance.
    pub suppression: f64,
    /// The divisor `C`.
    pub divisor: f64,
    /// `a`: the error grows as `r^(a−1)` in the rounds between yoke checks.
    pub round_exponent: u32,
    /// `b`: the error grows as `n^(b−1)` in the patches per block.
    pub patch_exponent: u32,
}

impl StorageFit {
    /// Logical error per patch per round at distance `d`, with `rounds`
    /// between yoke checks and `patches` per block.
    pub fn per_patch_round(&self, d: u32, rounds: u64, patches: u64) -> f64 {
        let mut scale = self.divisor;
        for _ in 0..d {
            scale *= self.suppression;
        }
        let mut x = 1.0;
        for _ in 1..self.round_exponent {
            x *= rounds as f64;
        }
        for _ in 1..self.patch_exponent {
            x *= patches as f64;
        }
        x / scale
    }

    /// The plain-patch fit equivalent to a threshold model
    /// `P(d) = a·(p/p*)^((d+1)/2)` per round, at physical error `p`. For
    /// example, the fit `quantum_match` measures from its own memory
    /// experiments.
    pub fn flat_from_threshold(prefactor: f64, threshold: f64, p: f64) -> StorageFit {
        let ratio_sqrt = exp(0.5 * ln(p / threshold));
        StorageFit {
            suppression: 1.0 / ratio_sqrt,
            divisor: 1.0 / (prefactor * ratio_sqrt),
            round_exponent: 1,
            patch_exponent: 1,
        }
    }
}

/// Fits for plain, row-yoked and grid-yoked patches.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct StorageModel {
    pub flat: StorageFit,
    pub rows: StorageFit,
    pub grid: StorageFit,
}

const fn fit(suppression: f64, divisor: f64, round_exponent: u32, patch_exponent: u32) -> StorageFit {
    StorageFit { suppression, divisor, round_exponent, patch_exponent }
}

impl StorageModel {
    /// arXiv:2312.04522, main text. SI1000 circuit noise at `p = 10⁻³`, fitted
    /// to one significant figure.
    pub fn si1000() -> StorageModel {
        StorageModel { flat: fit(3.0, 20.0, 1, 1), rows: fit(8.0, 500.0, 2, 2), grid: fit(50.0, 200_000.0, 4, 2) }
    }

    /// arXiv:2312.04522, supplementary note 1. Uniform depolarizing circuit
    /// noise at `p = 10⁻³`, fitted to one significant figure.
    pub fn uniform() -> StorageModel {
        StorageModel { flat: fit(4.0, 10.0, 1, 1), rows: fit(10.0, 2000.0, 2, 2), grid: fit(200.0, 20_000.0, 4, 2) }
    }

    /// The uniform-noise fits refitted to one and a half significant figures,
    /// as used by arXiv:2505.15917 (its storage figure). More pessimistic than
    /// [`StorageModel::uniform`].
    pub fn uniform_refit() -> StorageModel {
        StorageModel { flat: fit(3.5, 40.0, 1, 1), rows: fit(15.0, 100.0, 2, 2), grid: fit(150.0, 50_000.0, 4, 2) }
    }
}

/// Where a block sits: cold (dense, shared check hallway) or hot (every patch
/// faces an access hallway).
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Storage {
    Cold,
    Hot,
}

/// How a block's patches are bound together.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Yoke {
    /// Plain patches.
    None,
    /// One parity check along each row of patches (distance-2 outer code).
    Rows,
    /// Parity checks along rows and columns of a square (distance-4 outer
    /// code). Cold storage only.
    Grid,
}

/// A stored block: its shape, its distance and how well it keeps.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct StorageBlock {
    pub storage: Storage,
    pub yoke: Yoke,
    /// Patch distance.
    pub distance: u32,
    /// The distance a plain patch would need for the same target; sets the
    /// height of hot-storage access hallways.
    pub hallway_distance: u32,
    /// Patches per yoked group (a row), or `w²` for a grid.
    pub patches_per_group: u32,
    /// Groups (rows) in the block.
    pub groups: u32,
    /// Footprint, in patches, including hallways and workspace.
    pub rows: u32,
    pub cols: u32,
    /// Logical qubits the block stores.
    pub logical: u32,
    /// Rounds between yoke checks.
    pub rounds_between_checks: u64,
    /// Logical error per logical qubit per round.
    pub error_per_logical_round: f64,
}

impl StorageBlock {
    /// Physical qubits the block covers. A hot block's access hallways are
    /// `d0 + 1` patches tall, where `d0` is the plain-patch distance.
    pub fn physical_qubits(&self) -> u64 {
        let patch = patch_qubits(self.distance);
        match self.storage {
            Storage::Cold => patch * u64::from(self.rows) * u64::from(self.cols),
            Storage::Hot => {
                let hallways = u64::from(self.groups.div_ceil(2));
                let code_rows = u64::from(self.rows) - hallways;
                let hallway_patch = 2 * u64::from(self.distance + 1) * u64::from(self.hallway_distance + 1);
                patch * code_rows * u64::from(self.cols) + hallway_patch * hallways * u64::from(self.cols)
            }
        }
    }

    /// Physical qubits per stored logical qubit.
    pub fn qubits_per_logical(&self) -> f64 {
        self.physical_qubits() as f64 / f64::from(self.logical)
    }
}

const MAX_DISTANCE: u32 = 200;

/// The smallest plain-patch distance (at least 3) whose per-round error meets
/// `target`; `None` past distance 200.
pub fn flat_distance(model: &StorageModel, target: f64) -> Option<u32> {
    (3..=MAX_DISTANCE).find(|&d| model.flat.per_patch_round(d, 1, 1) <= target)
}

/// A rectangular block: `groups` rows of `patches_per_group` patches, plain or
/// row-yoked, at the smallest distance meeting `target` per logical qubit per
/// round. A row-yoked group must have an even number of patches, two of which
/// hold the yokes. `None` if no distance up to 200 meets the target, or the
/// shape stores nothing.
pub fn rectangular_block(model: &StorageModel, storage: Storage, yoke: Yoke, patches_per_group: u32, groups: u32, target: f64) -> Option<StorageBlock> {
    let yokes = match yoke {
        Yoke::None => 0,
        Yoke::Rows => 2,
        Yoke::Grid => return None,
    };
    if groups == 0 || (yokes > 0 && patches_per_group % 2 == 1) || patches_per_group <= yokes {
        return None;
    }
    let fit = if yokes == 0 { &model.flat } else { &model.rows };
    let per_group = patches_per_group - yokes;
    let mut d = 3;
    let (rounds, error) = loop {
        let rounds = match storage {
            Storage::Hot => 50 * u64::from(d),
            Storage::Cold => (u64::from(groups) * 8 + 2) * u64::from(d),
        };
        let error = fit.per_patch_round(d, rounds, u64::from(patches_per_group)) * f64::from(patches_per_group) / f64::from(per_group);
        if error <= target {
            break (rounds, error);
        }
        if d >= MAX_DISTANCE {
            return None;
        }
        d += 1;
    };
    let (mut rows, mut cols) = (groups, patches_per_group);
    match storage {
        Storage::Hot => {
            cols += 1;
            rows += groups.div_ceil(2);
        }
        Storage::Cold => {
            if yokes > 0 {
                rows += 1;
            }
        }
    }
    Some(StorageBlock {
        storage,
        yoke,
        distance: d,
        hallway_distance: flat_distance(model, target)?,
        patches_per_group,
        groups,
        rows,
        cols,
        logical: per_group * groups,
        rounds_between_checks: rounds,
        error_per_logical_round: error,
    })
}

/// A cold `w × w` grid-yoked block (`w` a positive multiple of 4), at the
/// smallest distance meeting `target` per logical qubit per round. It stores
/// `w² − 4w + 2` logical qubits in `(w+1)²` patches (one row and one column of
/// workspace).
pub fn grid_block(model: &StorageModel, width: u32, target: f64) -> Option<StorageBlock> {
    if width == 0 || !width.is_multiple_of(4) {
        return None;
    }
    let patches = width * width;
    let logical = patches - 4 * width + 2;
    let mut d = 3;
    let (rounds, error) = loop {
        let rounds = u64::from(d) * u64::from(width) * 25 + u64::from(d) * 4;
        let error = model.grid.per_patch_round(d, rounds, u64::from(patches)) * f64::from(patches) / f64::from(logical);
        if error <= target {
            break (rounds, error);
        }
        if d >= MAX_DISTANCE {
            return None;
        }
        d += 1;
    };
    Some(StorageBlock {
        storage: Storage::Cold,
        yoke: Yoke::Grid,
        distance: d,
        hallway_distance: flat_distance(model, target)?,
        patches_per_group: patches,
        groups: 1,
        rows: width + 1,
        cols: width + 1,
        logical,
        rounds_between_checks: rounds,
        error_per_logical_round: error,
    })
}

/// The cheapest block per logical qubit storing at most `max_logical` logical
/// qubits. Rectangular blocks try 4 to 200 patches per group and every group
/// count. Grid blocks try every width that fits. The first minimum wins, in
/// order of patches per group and then groups (or width).
pub fn best_block(model: &StorageModel, storage: Storage, yoke: Yoke, target: f64, max_logical: u32) -> Option<StorageBlock> {
    let mut best: Option<(f64, StorageBlock)> = None;
    let mut consider = |b: StorageBlock| {
        let rate = b.qubits_per_logical();
        if best.as_ref().is_none_or(|(r, _)| rate < *r) {
            best = Some((rate, b));
        }
    };
    match yoke {
        Yoke::Grid => {
            if storage == Storage::Hot {
                return None;
            }
            let mut w = 4;
            while w * w - 4 * w + 2 <= max_logical {
                if let Some(b) = grid_block(model, w, target) {
                    consider(b);
                }
                w += 4;
            }
        }
        Yoke::None | Yoke::Rows => {
            let yokes = if yoke == Yoke::Rows { 2 } else { 0 };
            for n in 4..=200u32 {
                if n <= yokes {
                    continue;
                }
                let per_group = n - yokes;
                for groups in 1..=max_logical.div_ceil(per_group) {
                    if per_group * groups > max_logical {
                        continue;
                    }
                    if let Some(b) = rectangular_block(model, storage, yoke, n, groups, target) {
                        consider(b);
                    }
                }
            }
        }
    }
    best.map(|(_, b)| b)
}

/// How a number of logical qubits is stored: whole blocks.
#[derive(Clone, Debug, PartialEq)]
pub struct StoragePlan {
    /// The blocks, with how many of each.
    pub blocks: Vec<(StorageBlock, u64)>,
    /// Logical qubits asked for.
    pub logical: u64,
}

impl StoragePlan {
    pub fn physical_qubits(&self) -> u64 {
        self.blocks.iter().map(|(b, k)| b.physical_qubits() * k).sum()
    }
    /// Logical qubits the blocks can hold (at least [`StoragePlan::logical`]).
    pub fn capacity(&self) -> u64 {
        self.blocks.iter().map(|(b, k)| u64::from(b.logical) * k).sum()
    }
}

/// Store `logical` qubits in whole cold blocks of one kind of yoke. A plan is
/// copies of one block shape, plus at most one block of another shape for the
/// remainder. Every pair of shapes holding at most `max_block` logical qubits
/// is tried, and the plan with the fewest physical qubits is kept. The shapes
/// are grid widths from 4, and rectangles of up to 200 patches per group and
/// 16 groups.
///
/// The fits were drawn from blocks of up to 256 patches. The source's own
/// layouts cap a block at 250 logical qubits. Larger blocks extrapolate the
/// fits, and they look cheaper than they have been shown to be.
pub fn store(model: &StorageModel, yoke: Yoke, target: f64, logical: u64, max_block: u32) -> Option<StoragePlan> {
    if logical == 0 {
        return Some(StoragePlan { blocks: Vec::new(), logical });
    }
    let mut shapes: Vec<StorageBlock> = Vec::new();
    match yoke {
        Yoke::Grid => {
            let mut w = 4;
            while w * w - 4 * w + 2 <= max_block {
                shapes.extend(grid_block(model, w, target));
                w += 4;
            }
        }
        Yoke::None | Yoke::Rows => {
            for n in 4..=200u32 {
                for groups in 1..=16 {
                    shapes.extend(rectangular_block(model, Storage::Cold, yoke, n, groups, target).filter(|b| b.logical <= max_block));
                }
            }
        }
    }
    let mut best: Option<(u64, StoragePlan)> = None;
    let mut consider = |blocks: Vec<(StorageBlock, u64)>| {
        let plan = StoragePlan { blocks, logical };
        let cost = plan.physical_qubits();
        if best.as_ref().is_none_or(|(c, _)| cost < *c) {
            best = Some((cost, plan));
        }
    };
    for main in &shapes {
        let cap = u64::from(main.logical);
        consider(vec![(*main, logical.div_ceil(cap))]);
        for rest in &shapes {
            let rest_cap = u64::from(rest.logical);
            let copies = logical.saturating_sub(rest_cap).div_ceil(cap);
            let mut blocks = vec![(*rest, 1)];
            if copies > 0 {
                blocks.insert(0, (*main, copies));
            }
            consider(blocks);
        }
    }
    best.map(|(_, plan)| plan)
}

/// One point on a published cultivation curve: a postselection cutoff, with
/// what it keeps and what it costs.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CultivationPoint {
    /// Injection (colour-code) distance, 3 or 5.
    pub d1: u32,
    /// Uniform circuit noise strength.
    pub p: f64,
    /// Kept shots and failures among them, as simulated.
    pub kept: u64,
    pub failures: u64,
    /// Expected qubit·rounds per kept state, retries included.
    pub volume: f64,
}

impl CultivationPoint {
    /// Logical error of a kept T state: twice the simulated S-state rate.
    pub fn error(&self) -> f64 {
        2.0 * self.failures as f64 / self.kept as f64
    }
}

/// `(d1, p, kept, failures, qubit·rounds)`. These are the representative
/// cutoffs of the end-to-end cultivation simulations of arXiv:2409.17595
/// (data: doi:10.5281/zenodo.13777072, Apache-2.0). Each curve escapes into a
/// distance-15 surface code. The volumes are the source's own: the circuit's
/// growth-weighted qubit·rounds divided by the keep rate. The error falls
/// down each curve as the volume rises.
#[rustfmt::skip]
const CULTIVATION: [(u32, f64, u64, u64, u32); 78] = [
    (3, 5e-4, 507_113_042, 76, 4434), (3, 5e-4, 567_320_249, 112, 3963), (3, 5e-4, 630_762_204, 260, 3564),
    (3, 5e-4, 702_651_030, 4541, 3200),
    (3, 1e-3, 225_419_907, 299, 8071), (3, 1e-3, 257_427_648, 383, 7067), (3, 1e-3, 289_504_256, 486, 6284),
    (3, 1e-3, 326_502_293, 677, 5572), (3, 1e-3, 367_643_011, 1251, 4948), (3, 1e-3, 423_967_169, 6213, 4291),
    (3, 1e-3, 471_154_637, 15137, 3861), (3, 1e-3, 524_756_158, 169_010, 3467),
    (3, 2e-3, 32_114_403, 400, 37284), (3, 2e-3, 40_189_645, 555, 29793), (3, 2e-3, 44_762_157, 671, 26749),
    (3, 2e-3, 53_846_749, 921, 22236), (3, 2e-3, 64_036_036, 1342, 18698), (3, 2e-3, 73_311_524, 1827, 16332),
    (3, 2e-3, 82_408_196, 2507, 14530), (3, 2e-3, 92_020_117, 3438, 13012), (3, 2e-3, 107_411_997, 5695, 11147),
    (3, 2e-3, 124_485_101, 11808, 9618), (3, 2e-3, 144_091_705, 20247, 8310), (3, 2e-3, 165_174_154, 33151, 7249),
    (3, 2e-3, 187_192_142, 58772, 6396), (3, 2e-3, 208_384_384, 107_491, 5746), (3, 2e-3, 237_794_538, 258_511, 5035),
    (3, 2e-3, 264_661_056, 825_206, 4524),
    (5, 5e-4, 92_480_789_479, 2, 18963), (5, 5e-4, 119_289_578_411, 3, 14702), (5, 5e-4, 133_864_311_235, 4, 13101),
    (5, 5e-4, 149_019_473_596, 16, 11769), (5, 5e-4, 168_245_418_061, 42, 10424), (5, 5e-4, 190_063_821_114, 101, 9227),
    (5, 5e-4, 213_505_407_573, 1225, 8214), (5, 5e-4, 237_773_048_111, 3798, 7376),
    (5, 1e-3, 5_072_794_335, 2, 127_164), (5, 1e-3, 6_326_363_510, 3, 101_967), (5, 1e-3, 7_316_858_650, 5, 88163),
    (5, 1e-3, 8_528_845_303, 7, 75635), (5, 1e-3, 11_303_894_346, 12, 57067), (5, 1e-3, 12_580_794_262, 17, 51275),
    (5, 1e-3, 14_509_395_372, 43, 44459), (5, 1e-3, 17_056_858_658, 85, 37819), (5, 1e-3, 19_856_355_181, 130, 32487),
    (5, 1e-3, 22_392_999_123, 361, 28807), (5, 1e-3, 24_893_295_331, 855, 25914), (5, 1e-3, 27_896_134_842, 1269, 23124),
    (5, 1e-3, 31_155_995_215, 2210, 20705), (5, 1e-3, 34_799_549_714, 3586, 18537), (5, 1e-3, 38_815_144_303, 5781, 16619),
    (5, 1e-3, 43_978_628_618, 15089, 14668), (5, 1e-3, 49_076_765_044, 88245, 13144),
    (5, 2e-3, 19_158_113, 1, 2_657_934), (5, 2e-3, 26_331_515, 2, 1_933_842), (5, 2e-3, 31_195_093, 3, 1_632_340),
    (5, 2e-3, 37_354_668, 6, 1_363_176), (5, 2e-3, 43_645_077, 10, 1_166_706), (5, 2e-3, 50_852_147, 13, 1_001_354),
    (5, 2e-3, 58_670_182, 18, 867_919), (5, 2e-3, 65_338_141, 27, 779_346), (5, 2e-3, 73_001_601, 36, 697_533),
    (5, 2e-3, 86_267_386, 63, 590_269), (5, 2e-3, 111_432_443, 95, 456_967), (5, 2e-3, 130_018_899, 134, 391_643),
    (5, 2e-3, 156_035_212, 195, 326_343), (5, 2e-3, 186_603_798, 300, 272_883), (5, 2e-3, 217_837_248, 492, 233_757),
    (5, 2e-3, 247_509_021, 780, 205_734), (5, 2e-3, 286_129_637, 1187, 177_965), (5, 2e-3, 338_075_961, 1890, 150_620),
    (5, 2e-3, 389_809_192, 3249, 130_631), (5, 2e-3, 438_097_321, 5457, 116_232), (5, 2e-3, 499_133_825, 10420, 102_019),
    (5, 2e-3, 555_534_346, 19226, 91661), (5, 2e-3, 628_514_709, 51577, 81018), (5, 2e-3, 700_358_340, 175_735, 72707),
    (5, 2e-3, 778_533_166, 1_971_026, 65406),
];

/// Every tabulated cultivation point, ordered by `(d1, p)` and then by rising
/// error.
pub fn cultivation_points() -> Vec<CultivationPoint> {
    CULTIVATION
        .iter()
        .map(|&(d1, p, kept, failures, volume)| CultivationPoint { d1, p, kept, failures, volume: f64::from(volume) })
        .collect()
}

/// A cultivated T state: its error, its expected cost and how it was grown.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CultivatedT {
    pub error: f64,
    /// Expected qubit·rounds per kept state.
    pub volume: f64,
    pub d1: u32,
}

/// The cheapest cultivated T state with error at most `target` at noise `p`.
/// The volume is interpolated log-linearly in `(ln error, ln volume)` between
/// the tabulated cutoffs that bracket the target. `None` if `p` is not one of
/// the simulated strengths (5·10⁻⁴, 10⁻³, 2·10⁻³) or no curve reaches the
/// target.
pub fn cultivate(p: f64, target: f64) -> Option<CultivatedT> {
    let points = cultivation_points();
    let mut best: Option<CultivatedT> = None;
    for d1 in [3, 5] {
        let curve: Vec<&CultivationPoint> = points.iter().filter(|q| q.d1 == d1 && q.p == p).collect();
        let Some(first) = curve.first() else { continue };
        if first.error() > target {
            continue;
        }
        // Points rise in error and fall in volume. Take the last point at or
        // below the target, and interpolate towards the next one up.
        let k = curve.iter().rposition(|q| q.error() <= target).unwrap_or(0);
        let lo = curve[k];
        let found = match curve.get(k + 1) {
            Some(hi) if target > lo.error() => {
                let t = (ln(target) - ln(lo.error())) / (ln(hi.error()) - ln(lo.error()));
                CultivatedT { error: target, volume: exp(ln(lo.volume) + t * (ln(hi.volume) - ln(lo.volume))), d1 }
            }
            _ => CultivatedT { error: lo.error(), volume: lo.volume, d1 },
        };
        if best.is_none_or(|b| found.volume < b.volume) {
            best = Some(found);
        }
    }
    best
}

/// A CCZ factory fed by cultivated T states: `t_states` T states in, one CCZ
/// state out.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct CczFactory {
    /// Footprint, in patches.
    pub width: u32,
    pub height: u32,
    /// Layers of lattice surgery per CCZ state.
    pub layers: u32,
    /// Rounds per layer, as a fraction of `d` (2/3 with temporally encoded
    /// lattice surgery).
    pub layer_fraction: f64,
    /// T states consumed per CCZ state.
    pub t_states: u32,
    /// Output error is `suppression · p_T²`.
    pub suppression: f64,
}

impl CczFactory {
    /// 8T-to-CCZ distillation (arXiv:1212.5069), in a 3 × 4 patch layout
    /// (arXiv:2409.17595, figure 24) of six temporally encoded layers.
    pub fn eight_t() -> CczFactory {
        CczFactory { width: 4, height: 3, layers: 6, layer_fraction: 2.0 / 3.0, t_states: 8, suppression: 28.0 }
    }

    /// Physical qubits at patch distance `d`.
    pub fn physical_qubits(&self, d: u32) -> u64 {
        u64::from(self.width * self.height) * patch_qubits(d)
    }

    /// Rounds per CCZ state. The T states are cultivated inside the
    /// footprint, each costing `t_volume` qubit·rounds. The lattice surgery
    /// follows.
    pub fn rounds_per_ccz(&self, d: u32, t_volume: f64) -> f64 {
        let cultivation = f64::from(self.t_states) * t_volume / self.physical_qubits(d) as f64;
        cultivation + f64::from(self.layers) * self.layer_fraction * f64::from(d)
    }

    /// Error of a CCZ state built from T states of error `t_error`.
    pub fn error(&self, t_error: f64) -> f64 {
        self.suppression * t_error * t_error
    }
}

/// What sets the pace of a step.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Bound {
    /// Lattice surgery: `d` rounds per step.
    Surgery,
    /// The CCZ supply.
    MagicStates,
    /// The classical control loop.
    Reaction,
}

/// The machine's clock.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Timing {
    /// Surface-code cycle, µs.
    pub cycle_us: f64,
    /// Control-system reaction time, µs.
    pub reaction_us: f64,
    /// Patch distance in the compute region.
    pub distance: u32,
    /// Rounds per CCZ state per factory.
    pub ccz_rounds: f64,
    pub factories: u32,
}

impl Timing {
    /// One lattice-surgery step: `d` cycles.
    pub fn surgery_us(&self) -> f64 {
        f64::from(self.distance) * self.cycle_us
    }
    /// Time between CCZ states from all factories together.
    pub fn ccz_period_us(&self) -> f64 {
        self.ccz_rounds * self.cycle_us / f64::from(self.factories)
    }
    /// A step: the slowest of surgery, CCZ supply and reaction.
    pub fn step_us(&self) -> f64 {
        self.surgery_us().max(self.ccz_period_us()).max(self.reaction_us)
    }
    /// Which of the three sets [`Timing::step_us`] (surgery on a tie).
    pub fn bound(&self) -> Bound {
        let step = self.step_us();
        if self.surgery_us() == step {
            Bound::Surgery
        } else if self.ccz_period_us() == step {
            Bound::MagicStates
        } else {
            Bound::Reaction
        }
    }
    /// An `n`-bit addition: `n − 1` CCZ states, then an uncompute of the same
    /// length.
    pub fn add_us(&self, n: u32) -> f64 {
        2.0 * f64::from(n.saturating_sub(1)) * self.step_us()
    }
    /// A table lookup over `w` address bits: `2ʷ − 1` layers, enough for its
    /// `2ʷ − w − 1` CCZ states.
    pub fn lookup_us(&self, w: u32) -> f64 {
        ((1u64 << w) - 1) as f64 * self.step_us()
    }
    /// A phaseup over `w` address bits: at most half a lookup.
    pub fn phaseup_us(&self, w: u32) -> f64 {
        self.lookup_us(w) / 2.0
    }
}

/// The parameters of the approximate-residue factoring algorithm of
/// arXiv:2505.15917 (its table 4).
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct FactoringParams {
    /// Bits in the modulus.
    pub n: u32,
    /// Ekerå–Håstad parameter.
    pub s: u32,
    /// Bits per prime of the residue system.
    pub prime_bits: u32,
    /// Window sizes of loops 1, 3 and 4.
    pub w1: u32,
    pub w3: u32,
    pub w4: u32,
    /// Bits of the truncated accumulator.
    pub accumulator_bits: u32,
    /// Shot failure rate from the approximation (simulated by the source).
    pub deviant: f64,
}

impl FactoringParams {
    /// The source's highlighted parameters for `n` ∈ {1024, 1536, 2048, 3072,
    /// 4096, 6144, 8192}.
    pub fn published(n: u32) -> Option<FactoringParams> {
        let (s, prime_bits, w1, w3, w4, accumulator_bits, deviant) = match n {
            1024 => (8, 18, 6, 3, 6, 28, 0.0287),
            1536 => (8, 21, 6, 3, 5, 31, 0.0183),
            2048 => (8, 21, 6, 3, 5, 33, 0.0125),
            3072 => (8, 21, 6, 3, 5, 35, 0.0091),
            4096 => (8, 24, 6, 3, 5, 36, 0.0080),
            6144 => (8, 24, 6, 3, 5, 39, 0.0042),
            8192 => (8, 24, 6, 3, 5, 40, 0.0040),
            _ => return None,
        };
        Some(FactoringParams { n, s, prime_bits, w1, w3, w4, accumulator_bits, deviant })
    }

    /// Input qubits `m = ⌈n(1/2 + 1/(2s))⌉ + ⌈n/(2s)⌉`.
    pub fn input_qubits(&self) -> u64 {
        let n = u64::from(self.n);
        let s = u64::from(self.s);
        (n * (s + 1)).div_ceil(2 * s) + n.div_ceil(2 * s)
    }

    /// Primes in the residue system: enough `ℓ`-bit primes to cover `n` bits
    /// per multiplication, `|P| = ⌈n·⌈m/w₁⌉/ℓ⌉`.
    pub fn primes(&self) -> u64 {
        (u64::from(self.n) * self.input_qubits().div_ceil(u64::from(self.w1))).div_ceil(u64::from(self.prime_bits))
    }

    /// Expected shots per factoring: `(s+1)/(1 − P_deviant)/0.99`.
    pub fn expected_shots(&self) -> f64 {
        f64::from(self.s + 1) / (1.0 - self.deviant) / 0.99
    }

    /// The operation tallies of one shot.
    pub fn tallies(&self) -> Vec<Subroutine> {
        let p = self.primes();
        let m = self.input_qubits();
        let len_m = u64::from(64 - m.leading_zeros());
        let (l, f) = (u64::from(self.prime_bits), u64::from(self.accumulator_bits));
        let win1 = m.div_ceil(u64::from(self.w1));
        let win3 = l.div_ceil(u64::from(self.w3));
        let win4 = l.div_ceil(u64::from(self.w4));
        let (w1, w3, w4) = (self.w1, self.w3, self.w4);
        let sub = |name, family, iterations: u64, halves: [u64; 3], add_bits, lookup_bits, phaseup_bits| Subroutine {
            name,
            family,
            iterations,
            half_adds: halves[0],
            half_lookups: halves[1],
            half_phaseups: halves[2],
            add_bits: add_bits as u32,
            lookup_bits,
            phaseup_bits,
        };
        vec![
            sub("loop1", 1, (p + 1) * win1, [2, 2, 0], l + len_m, w1, 0),
            sub("loop2", 2, p * len_m, [4, 0, 0], l + len_m, 0, 0),
            sub("loop3 startup", 3, p, [0, 2, 0], 0, 2 * w3, 0),
            sub("loop3 body", 3, p * (win3 - 2) * win3, [4, 2, 0], l + 1, 2 * w3, 0),
            sub("loop4", 4, p * win4, [5, 3, 2], f + 1, w4, w4),
            sub("unloop3 body", 5, p * (win3 - 2) * 2 * win3, [5, 3, 2], l + 1, 2 * w3, 2 * w3),
            sub("unloop3 cleanup", 5, p, [0, 0, 2], 0, 0, 2 * w3),
            sub("unloop2", 6, p * len_m, [4, 0, 0], l + len_m, 0, 0),
        ]
    }

    /// Peak logical qubits by the source's qubit tally (its table 5): loop 4
    /// holds `m + 3f + 2ℓ + len m`.
    pub fn peak_logical_qubits(&self) -> u64 {
        let m = self.input_qubits();
        let len_m = u64::from(64 - m.leading_zeros());
        m + 3 * u64::from(self.accumulator_bits) + 2 * u64::from(self.prime_bits) + len_m
    }

    /// Logical qubits by the source's reference code (its table 4 column):
    /// `m + (ℓ + len m) + ℓ + f + max(ℓ, f, ℓ + len m)`.
    pub fn table_logical_qubits(&self) -> u64 {
        let m = self.input_qubits();
        let len_m = u64::from(64 - m.leading_zeros());
        let (l, f) = (u64::from(self.prime_bits), u64::from(self.accumulator_bits));
        m + (l + len_m) + l + f + l.max(f).max(l + len_m)
    }

    /// Toffolis per shot, operation by operation:
    /// - an `n`-bit addition is `n − 1`;
    /// - a lookup over `w` bits is `2ʷ − w − 1`;
    /// - a phaseup over `w` bits is `⌈√(2ʷ)⌉`, so 8 at 6 bits, as the source's
    ///   text states.
    pub fn toffolis_per_shot(&self) -> f64 {
        self.tallies()
            .iter()
            .map(|s| {
                let add = f64::from(s.add_bits.saturating_sub(1));
                let lookup = ((1u64 << s.lookup_bits) - u64::from(s.lookup_bits) - 1) as f64;
                let phaseup = isqrt_ceil(1u64 << s.phaseup_bits) as f64;
                s.adds() * add + s.lookups() * lookup + s.phaseups() * phaseup
            })
            .sum()
    }

    /// Toffolis per shot by the source's reference counting rule, the rule
    /// behind its published count:
    /// - every operation in a subroutine is priced at that subroutine's
    ///   largest size;
    /// - an addition costs its full register width;
    /// - a lookup over `N` entries costs `N − bitlen(N) − 1`.
    ///
    /// The rule prices a phaseup over `N` entries at
    /// `2^(⌊N₁/2⌋−N₁−1) + 2^(⌊N₂/2⌋−N₂−1)`, a quarter of a Toffoli at 6 address
    /// bits. [`PhaseupCount::AsStated`] substitutes the text's `⌈√N⌉`.
    pub fn toffolis_per_shot_source_rule(&self, phaseups: PhaseupCount) -> f64 {
        let tallies = self.tallies();
        (1..=6)
            .map(|family| {
                let rows: Vec<&Subroutine> = tallies.iter().filter(|s| s.family == family).collect();
                let add_bits = rows.iter().map(|s| s.add_bits).max().unwrap_or(0);
                let lookup_bits = rows.iter().filter(|s| s.half_lookups > 0).map(|s| s.lookup_bits).max().unwrap_or(0);
                let phaseup_bits = rows.iter().filter(|s| s.half_phaseups > 0).map(|s| s.phaseup_bits).max().unwrap_or(0);
                let (adds, lookups, flips) = rows.iter().fold((0.0, 0.0, 0.0), |(a, l, p), s| (a + s.adds(), l + s.lookups(), p + s.phaseups()));
                let mut total = adds * f64::from(add_bits);
                if lookups > 0.0 {
                    let big_n = 1u64 << lookup_bits;
                    total += lookups * (big_n - u64::from(64 - big_n.leading_zeros()) - 1) as f64;
                }
                if flips > 0.0 {
                    total += flips
                        * match phaseups {
                            PhaseupCount::AsCoded => {
                                let bits = 64 - (1u64 << phaseup_bits).leading_zeros() as i32;
                                let (n1, n2) = (bits / 2, bits - bits / 2);
                                pow2(n1 / 2 - n1 - 1) + pow2(n2 / 2 - n2 - 1)
                            }
                            PhaseupCount::AsStated => isqrt_ceil(1u64 << phaseup_bits) as f64,
                        };
                }
                total
            })
            .sum()
    }
}

/// How the source rule prices a phaseup.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum PhaseupCount {
    /// As the reference counting code evaluates it: a quarter of a Toffoli at
    /// 6 address bits.
    AsCoded,
    /// As the text states: `⌈√N⌉`, so 8 at 6 address bits.
    AsStated,
}

fn pow2(e: i32) -> f64 {
    let mut x = 1.0;
    for _ in 0..e.unsigned_abs() {
        x = if e < 0 { x / 2.0 } else { x * 2.0 };
    }
    x
}

fn isqrt_ceil(x: u64) -> u64 {
    let mut r = (x as f64).sqrt() as u64;
    while r * r > x {
        r -= 1;
    }
    while r * r < x {
        r += 1;
    }
    r
}

/// One subroutine of the factoring algorithm: how often it runs and what each
/// iteration does. Operation counts are kept in halves, because loop 4 and
/// unloop 3 average two and a half additions per iteration.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Subroutine {
    pub name: &'static str,
    /// Which of the six loops it belongs to (1 to 6: loop 1, 2, 3, 4,
    /// unloop 3, unloop 2).
    pub family: u32,
    pub iterations: u64,
    pub half_adds: u64,
    pub half_lookups: u64,
    pub half_phaseups: u64,
    /// Largest addition, in bits.
    pub add_bits: u32,
    /// Largest lookup and phaseup, in address bits.
    pub lookup_bits: u32,
    pub phaseup_bits: u32,
}

impl Subroutine {
    pub fn adds(&self) -> f64 {
        (self.iterations * self.half_adds) as f64 / 2.0
    }
    pub fn lookups(&self) -> f64 {
        (self.iterations * self.half_lookups) as f64 / 2.0
    }
    pub fn phaseups(&self) -> f64 {
        (self.iterations * self.half_phaseups) as f64 / 2.0
    }
}

/// The physical assumptions of a whole-machine estimate.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Machine {
    /// Uniform circuit noise strength.
    pub p: f64,
    pub cycle_us: f64,
    pub reaction_us: f64,
    pub storage: StorageModel,
    /// Target logical error per logical qubit per round, for every stored
    /// qubit.
    pub target: f64,
    pub factory: CczFactory,
    pub factories: u32,
    /// Columns of lattice-surgery workspace beside the factories.
    pub workspace_cols: u32,
    /// Error target for each T state fed to a factory.
    pub t_target: f64,
    /// Largest cold block, in logical qubits (see [`store`]).
    pub max_block: u32,
}

impl Machine {
    /// The source's machine:
    /// - `p = 10⁻³`, a 1 µs cycle and a 10 µs reaction time;
    /// - refitted storage at 10⁻¹⁵ per logical qubit per round;
    /// - six 8T-to-CCZ factories with three columns of workspace;
    /// - 10⁻⁷ T states.
    pub fn published() -> Machine {
        Machine {
            p: 1e-3,
            cycle_us: 1.0,
            reaction_us: 10.0,
            storage: StorageModel::uniform_refit(),
            target: 1e-15,
            factory: CczFactory::eight_t(),
            factories: 6,
            workspace_cols: 3,
            t_target: 1e-7,
            max_block: 250,
        }
    }
}

/// A whole-machine estimate of the factoring algorithm.
#[derive(Clone, Debug, PartialEq)]
pub struct FactoringEstimate {
    pub params: FactoringParams,
    /// Logical qubits: the input register stays cold, the rest is hot.
    pub cold_logical: u64,
    pub hot_logical: u64,
    /// Hot-patch distance and the cold storage plan.
    pub hot_distance: u32,
    pub cold: StoragePlan,
    /// Physical qubits in each region.
    pub cold_qubits: u64,
    pub hot_qubits: u64,
    pub compute_qubits: u64,
    /// The T state each factory consumes, and its CCZ output error.
    pub t_state: CultivatedT,
    pub ccz_error: f64,
    /// Rounds per CCZ per factory.
    pub ccz_rounds: f64,
    pub timing: Timing,
    /// Hours per shot at the clock's lower bound.
    pub hours_per_shot: f64,
    /// Hours per shot with every operation rounded up to a whole millisecond,
    /// the source's allowance for moving qubits in and out.
    pub hours_per_shot_rounded: f64,
    /// Expected shots, the chance a shot has no logical error, and expected
    /// days per factoring (shots lost to logical errors included, at the
    /// rounded durations).
    pub expected_shots: f64,
    pub shot_survival: f64,
    pub days: f64,
    pub toffolis_per_shot: f64,
}

impl FactoringEstimate {
    pub fn physical_qubits(&self) -> u64 {
        self.cold_qubits + self.hot_qubits + self.compute_qubits
    }
}

/// Why a whole-machine estimate could not be made.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum ArchError {
    /// No plain-patch distance up to 200 meets the storage target.
    Distance,
    /// No yoked layout meets the storage target.
    Storage,
    /// No cultivation curve at this noise strength reaches the T target.
    Cultivation,
}

impl core::fmt::Display for ArchError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        f.write_str(match self {
            ArchError::Distance => "no patch distance up to 200 meets the storage target",
            ArchError::Storage => "no yoked storage layout meets the storage target",
            ArchError::Cultivation => "no cultivation curve at this noise strength reaches the T target",
        })
    }
}

impl std::error::Error for ArchError {}

/// Re-run the factoring estimate on `machine`.
///
/// - **Cold storage.** The `m` input qubits sit in whole grid-yoked blocks,
///   via [`store`].
/// - **Hot storage.** The rest of the peak sits in plain patches.
/// - **Compute region.** `factories` factories are stacked beside the
///   workspace columns.
/// - **Time.** Every operation in the tallies runs on the [`Timing`] clock,
///   which is set by the cultivated T states.
/// - **Survival.** A shot survives when no stored logical qubit fails for its
///   whole duration.
pub fn factoring_estimate(params: &FactoringParams, machine: &Machine) -> Result<FactoringEstimate, ArchError> {
    let d = flat_distance(&machine.storage, machine.target).ok_or(ArchError::Distance)?;
    let cold_logical = params.input_qubits();
    let hot_logical = params.peak_logical_qubits() - cold_logical;
    let cold = store(&machine.storage, Yoke::Grid, machine.target, cold_logical, machine.max_block).ok_or(ArchError::Storage)?;
    let t_state = cultivate(machine.p, machine.t_target).ok_or(ArchError::Cultivation)?;
    let ccz_rounds = machine.factory.rounds_per_ccz(d, t_state.volume);
    let timing = Timing { cycle_us: machine.cycle_us, reaction_us: machine.reaction_us, distance: d, ccz_rounds, factories: machine.factories };
    let whole_ms = |us: f64| (us / 1000.0).ceil() * 1000.0;
    let (mut micros, mut micros_rounded) = (0.0, 0.0);
    for s in params.tallies() {
        let ops = [(s.adds(), timing.add_us(s.add_bits)), (s.lookups(), timing.lookup_us(s.lookup_bits)), (s.phaseups(), timing.phaseup_us(s.phaseup_bits))];
        for (count, us) in ops {
            micros += count * us;
            micros_rounded += count * whole_ms(us);
        }
    }
    let compute_patches = u64::from((machine.factory.width + machine.workspace_cols) * machine.factory.height * machine.factories);
    let (cold_qubits, hot_qubits, compute_qubits) = (cold.physical_qubits(), hot_logical * patch_qubits(d), compute_patches * patch_qubits(d));
    let stored = (cold_logical + hot_logical + compute_patches) as f64;
    let rounds = micros_rounded / machine.cycle_us;
    let shot_survival = exp(-(stored * rounds * machine.target));
    let expected_shots = params.expected_shots();
    let (hours_per_shot, hours_per_shot_rounded) = (micros / 3.6e9, micros_rounded / 3.6e9);
    Ok(FactoringEstimate {
        params: *params,
        cold_logical,
        hot_logical,
        hot_distance: d,
        cold,
        cold_qubits,
        hot_qubits,
        compute_qubits,
        t_state,
        ccz_error: machine.factory.error(t_state.error),
        ccz_rounds,
        timing,
        hours_per_shot,
        hours_per_shot_rounded,
        expected_shots,
        shot_survival,
        days: hours_per_shot_rounded * expected_shots / shot_survival / 24.0,
        toffolis_per_shot: params.toffolis_per_shot(),
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    fn rsa() -> FactoringParams {
        FactoringParams::published(2048).unwrap()
    }

    #[test]
    fn refitted_storage_gives_the_published_densities() {
        let m = StorageModel::uniform_refit();
        assert_eq!(flat_distance(&m, 1e-15), Some(25));
        assert_eq!(patch_qubits(25), 1352);
        let g = grid_block(&m, 16, 1e-15).unwrap();
        assert_eq!((g.distance, g.logical, g.rows, g.cols), (11, 194, 17, 17));
        assert_eq!(g.physical_qubits(), 83_232);
        assert!((g.qubits_per_logical() - 429.03).abs() < 0.01, "{}", g.qubits_per_logical());
        // Its stated rounding: 430 per logical qubit.
        assert_eq!(1280 * 430 + 131 * 1352 + 7 * 18 * 1352, 897_864);
    }

    #[test]
    fn uniform_grids_reach_the_supplements_density() {
        // Supplementary note 1: about 350 physical qubits per logical qubit
        // at about 10⁻¹² per round.
        let g = best_block(&StorageModel::uniform(), Storage::Cold, Yoke::Grid, 1e-12, 250).unwrap();
        assert!((340.0..370.0).contains(&g.qubits_per_logical()), "{}", g.qubits_per_logical());
    }

    #[test]
    fn the_layout_search_matches_the_source_tool_everywhere() {
        let data = include_str!("../tests/data/yoked_footprint_oracle.csv");
        let mut checked = 0;
        for line in data.lines().skip(1) {
            let f: Vec<&str> = line.split(',').collect();
            let model = match f[0] {
                "si1000" => StorageModel::si1000(),
                "uniform" => StorageModel::uniform(),
                "uniform_refit" => StorageModel::uniform_refit(),
                other => panic!("{other}"),
            };
            let storage = if f[1] == "hot" { Storage::Hot } else { Storage::Cold };
            let target: f64 = f[3].parse().unwrap();
            let got = match f[2] {
                "0" => best_block(&model, storage, Yoke::None, target, 250),
                "2" => best_block(&model, storage, Yoke::Rows, target, 250),
                "64" => grid_block(&model, 16, target),
                other => panic!("{other}"),
            }
            .unwrap_or_else(|| panic!("no layout for {line}"));
            let want: Vec<u64> = f[4..11].iter().map(|x| x.parse().unwrap()).collect();
            let have = [got.distance, got.groups, got.patches_per_group, got.rows, got.cols, got.logical].map(u64::from);
            assert_eq!(&have[..], &want[..6], "{line}");
            assert_eq!(got.rounds_between_checks, want[6], "{line}");
            // The tool leaves the hallway distance unset for grids, which
            // have no hallways.
            if f[11] != "None" {
                assert_eq!(got.hallway_distance.to_string(), f[11], "{line}");
            }
            let rate: f64 = f[12].parse().unwrap();
            assert_eq!(got.qubits_per_logical(), rate, "{line}");
            let lerp: f64 = f[13].parse().unwrap();
            assert!((got.error_per_logical_round / lerp - 1.0).abs() < 1e-12, "{line}");
            checked += 1;
        }
        assert_eq!(checked, 915);
    }

    #[test]
    fn whole_blocks_cost_more_than_the_rate_suggests() {
        let m = StorageModel::uniform_refit();
        let plan = store(&m, Yoke::Grid, 1e-15, 1280, 250).unwrap();
        assert!(plan.capacity() >= 1280);
        let widths: Vec<(u32, u64)> = plan.blocks.iter().map(|(b, k)| (b.rows - 1, *k)).collect();
        assert_eq!(widths, vec![(16, 7)]);
        assert_eq!(plan.physical_qubits(), 582_624);
        // Past the simulated block sizes the fits promise more: a 32-wide grid
        // at 410 qubits per logical qubit.
        let wide = store(&m, Yoke::Grid, 1e-15, 1280, 1000).unwrap();
        assert!(wide.physical_qubits() < 1280 * 430);
        assert!(wide.blocks.iter().any(|(b, _)| b.rows - 1 == 32));
    }

    #[test]
    fn cultivation_reads_the_published_curves() {
        let t = cultivate(1e-3, 1e-7).unwrap();
        assert_eq!(t.d1, 5);
        assert!((22_000.0..23_200.0).contains(&t.volume), "{}", t.volume);
        // Well above the distance-5 curve's reach, distance 3 is cheaper.
        let cheap = cultivate(1e-3, 1e-4).unwrap();
        assert_eq!(cheap.d1, 3);
        // Below every curve: nothing.
        assert_eq!(cultivate(1e-3, 1e-12), None);
        // Unsimulated noise: refused.
        assert_eq!(cultivate(7e-4, 1e-6), None);
        // Each curve rises in error as it falls in volume.
        let pts = cultivation_points();
        for w in pts.windows(2) {
            if (w[0].d1, w[0].p) == (w[1].d1, w[1].p) {
                assert!(w[1].error() > w[0].error() && w[1].volume < w[0].volume);
            }
        }
    }

    #[test]
    fn the_factory_and_clock_reproduce_the_sources_arithmetic() {
        let f = CczFactory::eight_t();
        assert_eq!(f.physical_qubits(25), 16_224);
        // 8 × 30,000 / 16,224 + 6 × (2/3) × 25 = 114.79 rounds ("114.7").
        assert!((f.rounds_per_ccz(25, 30_000.0) - 114.79).abs() < 0.01);
        assert!(f.error(1e-7) < 1e-12);
        // Rounded up to 150 for slack: six factories match the 25 µs surgery
        // period.
        let t = Timing { cycle_us: 1.0, reaction_us: 10.0, distance: 25, ccz_rounds: 150.0, factories: 6 };
        assert_eq!((t.surgery_us(), t.ccz_period_us(), t.step_us()), (25.0, 25.0, 25.0));
        assert_eq!(t.bound(), Bound::Surgery);
        assert_eq!(t.add_us(33), 1600.0);
        assert_eq!(t.lookup_us(6), 1575.0);
        let slow = Timing { reaction_us: 40.0, ..t };
        assert_eq!(slow.bound(), Bound::Reaction);
        let starved = Timing { factories: 2, ..t };
        assert_eq!(starved.bound(), Bound::MagicStates);
    }

    #[test]
    fn the_logical_counts_reproduce_the_published_table() {
        let p = rsa();
        assert_eq!(p.input_qubits(), 1280);
        assert_eq!(p.primes(), 20_871);
        assert_eq!(p.table_logical_qubits(), 1399);
        assert_eq!(p.peak_logical_qubits(), 1432);
        // The published count: 6.5·10⁹ Toffolis per factoring.
        let per_shot = p.toffolis_per_shot_source_rule(PhaseupCount::AsCoded);
        let published = per_shot * f64::from(p.s + 1) / (1.0 - p.deviant);
        assert!((6.45e9..6.55e9).contains(&published), "{published:e}");
        // The reference tally rounds loop 4's half operations to integers; held
        // in halves here, it is 29.5 Toffolis lower.
        assert_eq!(per_shot, 712_920_014.5);
        // Every published row. The reference implementation's counts and keep
        // rates (its simulated `1 − P_deviant`, printed rounded in the table)
        // agree to its rounding of half operations.
        for (n, keep, tofs, qubits) in [
            (1024, 0.971_311_569_213_867_2, 1_083_979_446.0, 742),
            (1536, 0.981_724_202_632_904, 3_117_105_798.0, 1074),
            (2048, 0.987_522_467_970_848_1, 6_497_351_306.0, 1399),
            (3072, 0.990_862_101_316_452, 18_500_632_076.0, 2043),
            (4096, 0.992_021_542_042_493_8, 40_261_413_896.0, 2692),
            (6144, 0.995_758_056_640_625, 119_251_229_120.0, 3978),
            (8192, 0.996_010_778_006_166_2, 266_770_779_692.0, 5261),
        ] {
            let q = FactoringParams::published(n).unwrap();
            assert_eq!(q.table_logical_qubits(), qubits, "n={n}");
            assert!((1.0 - q.deviant - keep).abs() < 5e-5, "n={n}");
            let t = q.toffolis_per_shot_source_rule(PhaseupCount::AsCoded) * f64::from(q.s + 1) / keep;
            assert!((t / tofs - 1.0).abs() < 1e-7, "n={n}: {t} vs {tofs}");
        }
    }

    #[test]
    fn the_phaseup_rule_undercounts() {
        let p = rsa();
        let coded = p.toffolis_per_shot_source_rule(PhaseupCount::AsCoded);
        let stated = p.toffolis_per_shot_source_rule(PhaseupCount::AsStated);
        // 1,481,841 six-bit phaseups at 8 rather than 0.25, and 104,355
        // five-bit ones at 6.
        assert_eq!(stated - coded, 1_481_841.0 * 7.75 + 104_355.0 * 5.75);
        assert!((0.016..0.018).contains(&(stated / coded - 1.0)), "{}", stated / coded - 1.0);
        // Counted operation by operation, at full sizes, the total lands
        // between the two.
        let own = p.toffolis_per_shot();
        assert!(own < stated, "{own} {stated}");
    }

    #[test]
    fn the_stated_durations_give_ten_point_six_hours_not_twelve() {
        // 2 ms per addition and per lookup, 1 ms per phaseup.
        let ms: f64 = rsa().tallies().iter().map(|s| 2.0 * (s.adds() + s.lookups()) + s.phaseups()).sum();
        let hours = ms / 3.6e6;
        assert!((hours - 10.621).abs() < 0.001, "{hours}");
        assert!(hours < 12.07 / 1.13);
    }

    #[test]
    fn the_whole_machine_stays_under_a_million_qubits_and_a_week() {
        let e = factoring_estimate(&rsa(), &Machine::published()).unwrap();
        assert_eq!(e.hot_distance, 25);
        assert_eq!(e.hot_logical, 152);
        assert_eq!(e.compute_qubits, 170_352);
        assert_eq!(e.cold_qubits, 582_624);
        assert_eq!(e.physical_qubits(), 958_480);
        assert!(e.physical_qubits() < 1_000_000);
        assert_eq!(e.t_state.d1, 5);
        assert!(e.ccz_error < 1e-12);
        assert_eq!(e.timing.bound(), Bound::Surgery);
        assert!(e.days < 7.0, "{}", e.days);
        // At the source's own hot count the total is 926,256 rather than 897,864.
        assert_eq!(1280 * 430 + 152 * 1352 + 7 * 18 * 1352, 926_256);
    }

    #[test]
    fn a_measured_fit_converts_to_a_storage_fit() {
        // a = 0.023, p* = 0.01 at p = 10⁻³: Λ = √10 and C = 1/(a·√0.1).
        let f = StorageFit::flat_from_threshold(0.023, 0.01, 1e-3);
        assert!((f.suppression - 10f64.sqrt()).abs() < 1e-12);
        assert!((f.per_patch_round(9, 1, 1) - 0.023 * 1e-5).abs() < 1e-18);
        let m = StorageModel { flat: f, ..StorageModel::uniform_refit() };
        assert_eq!(flat_distance(&m, 1e-15), Some(26));
    }
}