qsimlab.shor

qsimlab.shor

Shor's algorithm at gate level: factoring runs, oracle circuits, resource counts, cost laws.

Provisional (phase 2, see python/API.md §1 and §8).

Every run simulates the semiclassical order-finding circuit exactly (one recycled control qubit, Griffiths–Niu phase corrections) with a gate-level oracle built from X/CNOT/Toffoli gates, on the bit-sliced branch engine of research/shor/shor.md (or the dense/sparse engines for the permutation and Beauregard oracles):

>>> import qsimlab.shor as shor
>>> r = shor.factor(1_005_973, seed=1)            # 20-bit N, 88 qubits
>>> r.factors, r.base, r.order
((997, 1009), 980062, 41832)
>>> r.qubits, r.toffoli_gates
(88, 70720)
>>> bool((r.runs[0].support_trace <= r.runs[0].predicted_support).all())
True

The cost of a run is set by the multiplicative order r of the base, not by N (the support law, research/theory/theory-shor.md T1): before round i the state holds at most B_i = min(2^i, r / gcd(r, 2^(t-i))) branches. factor predicts the peak memory from it before running and refuses runs over budget with ResourceLimitError. The order used for that prediction is computed classically (Pollard–Brent factoring); the simulated algorithm never sees it. This is exact simulation, not a factoring speed-up: for a generic semiprime and a random base, r ≈ N/c and the work grows like N · n³.

Oracles (oracle=, see ORACLES):

• "windowed-opt" (default): Gidney's windowed table-lookup multiplier with the superoptimised blocks of research/shor/superopt.md; 4n + 4 + w qubits;
• "windowed-mbu-lookup" / "windowed-mbu": measurement-based uncomputation (research/shor/mbu-shor.md): fewer Toffolis, X-basis measurements with classical fix-ups;
• "windowed", "ripple": the round-4 windowed oracle and the Cuccaro ripple-carry oracle (3n + 4 qubits);
• "beauregard": Draper/Beauregard QFT arithmetic, 2n + 3 qubits, dense state (≈ 10-bit N);
• "permutation": a classical lookup table as the oracle, n + 1 qubits (measurement statistics only, not a compilable circuit);
• "ge" / "eh": Gidney–Ekerå exponent windowing (research/shor/ge-shor.md) for Shor's order finding, or the Ekerå–Håstad short-discrete-log schedule (1.5n exponent bits, lattice post-processing; balanced semiprimes).

ShorRun dataclass

ShorRun(
    base: int,
    measured: Union[int, Tuple[int, ...]],
    order: Optional[int],
    factors: Optional[Tuple[int, int]],
    qubits: int,
    total_gates: int,
    toffoli_gates: int,
    measurements: int,
    peak_support: int,
    peak_amplitude_bytes: int,
    gate_branch_ops: int,
    true_order: int,
    predicted_peak_support: int,
    predicted_bytes: int,
    engine: str,
    wall_time: float,
    support_trace: Optional[ndarray] = None,
    p1_trace: Optional[ndarray] = None,
)

One order-finding (or Ekerå–Håstad) run.

• base: the base a; measured: the measured integer (bit i = round i), a tuple (k, j) for "eh" (one integer per exponent register);
• order / factors: what the classical post-processing recovered (None if it failed — that is the classical part of Shor failing, e.g. an odd order or a^(r/2) ≡ −1);
• qubits, total_gates (every gate of the run, control gates included), toffoli_gates, measurements (mid-circuit X-basis measurements of the MBU oracles);
• peak_support: largest number of stored branches; peak_amplitude_bytes: amplitude memory of dense/sparse engines, peak branch count of the GE engine; gate_branch_ops: gates × branches evaluated (the work of the sliced engine);
• true_order: the multiplicative order computed classically (for the memory guard); predicted_peak_support / predicted_bytes: the guard's prediction;
• support_trace: branches stored before every round (sliced engine, trace=True), p1_trace: P(control = 1) of every round.

predicted_support property

predicted_support: Optional[ndarray]

B_i for the rounds in support_trace (None without a trace).

FactorResult dataclass

FactorResult(
    engine: str,
    components: List[Tuple[int, int, str]],
    seed: int,
    precision: str,
    wall_time: float,
    explanation: Optional[Explanation] = None,
    N: int = 0,
    oracle: str = "",
    factors: Optional[Tuple[int, int]] = None,
    base: Optional[int] = None,
    order: Optional[int] = None,
    measured: Union[int, Tuple[int, ...], None] = None,
    qubits: int = 0,
    total_gates: int = 0,
    toffoli_gates: int = 0,
    peak_support: int = 0,
    budget: int = 0,
    runs: List[ShorRun] = list(),
)

Bases: Result

Result of factor (a qsimlab.sim.Result).

factors is (p, q) with p ≤ q or None; base, order, measured and the resource fields describe the last (successful) run; runs lists every attempt (ShorRun). components has one (qubits, gates, engine) entry per run.

ResourceCounts dataclass

ResourceCounts(
    N: int,
    base: int,
    oracle: str,
    qubits: int,
    rounds: int,
    gate_level: bool,
    oracle_gates: int,
    toffoli: int,
    cnot: int,
    x: int,
    measurements: int,
    fixups: int,
    slice_steps: Optional[int] = None,
    per_round_gates: Optional[ndarray] = None,
    per_round_toffoli: Optional[ndarray] = None,
)

Whole-run counts of the oracle blocks of one order-finding (or EH) run.

oracle_gates counts every operation of the rounds controlled-U blocks (X, CNOT, Toffoli, X-basis measurements and their Z/CZ fix-ups); toffoli, cnot, x, measurements, fixups break it down. The semiclassical control adds at most 2 H, 1 phase and 1 recycling X per round (control_gates_max). For the MBU and GE oracles the counts depend on the (fixed, seeded) outcome stream of the engine. slice_steps: GE only. gate_level is False for the permutation oracle (no gates).

OracleCircuit dataclass

OracleCircuit(
    circuit: Circuit,
    N: int,
    a: int,
    oracle: str,
    control: int,
    work: List[int],
    num_qubits: int,
    gates: int,
    toffolis: int,
)

A controlled-U_a block: |c⟩|x⟩|0…0⟩ → |c⟩|a^c·x mod N⟩|0…0⟩ for x < N.

control is qubit 0, work[k] holds bit k of x, every other qubit is an ancilla that starts and ends in |0⟩.

SupportPrediction dataclass

SupportPrediction(
    order: int,
    nu: int,
    rounds: int,
    peak: int,
    sum: int,
    engine: str,
    predicted_bytes: int,
    bounds: ndarray,
)

The support law for (N, a): order r, nu = ν₂(r), rounds t = 2n, bounds (B_i), peak (max_i B_i), sum (Σ B_i; the sliced engine's work is ≈ 2 Ḡ Σ B_i gate·branch steps for Ḡ gates per round), and the engine and peak bytes factor would predict.

NoisySuccess dataclass

NoisySuccess(
    engine: str,
    components: List[Tuple[int, int, str]],
    seed: int,
    precision: str,
    wall_time: float,
    explanation: Optional[Explanation] = None,
    N: int = 0,
    base: int = 0,
    p: float = 0.0,
    noise: str = "",
    trajectories: int = 0,
    locations: int = 0,
    order: int = 0,
    success: float = 0.0,
    ci: Tuple[float, float] = (0.0, 1.0),
    order_rate: float = 0.0,
    peak_rate: float = 0.0,
    capped: int = 0,
    mean_faults: float = 0.0,
    faults: Optional[ndarray] = None,
    factor_ok: Optional[ndarray] = None,
    order_ok: Optional[ndarray] = None,
    peak_ok: Optional[ndarray] = None,
    measured: Optional[List[Optional[int]]] = None,
)

Bases: Result

Result of noisy_success (a qsimlab.sim.Result).

• success: fraction of trajectories whose post-processing returned a factor, with ci its Wilson 95 % interval;
• order_rate: the true order was recovered; peak_rate: the measured y is a "good" outcome, |y/2^t − s/r| < 1/(2r²) (the metric of research/shor/shor-noise.md);
• locations: fault locations L of the circuit, mean_faults;
• capped: trajectories abandoned because their support exceeded cap (counted as failures);
• per-trajectory arrays faults, factor_ok, order_ok, peak_ok.

support_bounds

support_bounds(order: int, rounds: int) -> ndarray

The support law (research/theory/theory-shor.md T1): B_i = min(2^i, r / gcd(r, 2^(t-i))) for i = 0 … t-1 — an upper bound on the number of branches before round i, attained except with probability ≤ 4/r_odd per round.

>>> support_bounds(12, 8).tolist()        # r = 12 = 4 · 3, t = 8 rounds
[1, 2, 3, 3, 3, 3, 3, 6]

factor

factor(
    N: int,
    *,
    oracle: str = "windowed-opt",
    window: Optional[int] = None,
    exponent_window: Optional[int] = None,
    base: Optional[int] = None,
    precision: str = "f64",
    seed: Optional[int] = None,
    tries: int = 10,
    budget: Union[Budget, int, str, None] = None,
    engine: str = "auto",
    trace: bool = True,
    threads: Optional[int] = None,
) -> FactorResult

Factor N by simulating Shor's algorithm at gate level.

Up to tries runs; each draws a random base a (or uses base), simulates the whole semiclassical order-finding circuit and post-processes the measured integer (continued fractions, small multiples), stopping at the first run that yields a factor. With the same seed the bases and measured integers are those of qsim run shor --seed <seed> (the CLI of the Rust crate).

• window: lookup window w of the windowed oracles (default 4) or the multiplicand window w_m of "ge"/"eh" (default 3); exponent_window: w_e of "ge"/"eh" (default 2).
• precision: amplitudes in "f64" or "f32" (≈ 35 % less memory; distributions within 1e-7 of f64, research/shor/shor.md).
• budget: refuse (ResourceLimitError, with needed/limit) any run whose predicted peak memory exceeds it. Default: min(32 GiB, half the physical memory).
• engine: "auto" (sliced branches for the gate-level X/CNOT/CCX oracles, the cheaper of fused dense/sparse for "permutation", dense for "beauregard"), or "sliced", "dense", "sparse".

N must be odd, composite, not a prime power, and below 2^62.

>>> r = factor(143, oracle="ripple", seed=3)
>>> r.factors
(11, 13)
>>> pow(r.base, r.order, 143)
1

resource_counts

resource_counts(
    N: int,
    oracle: str = "windowed-opt",
    *,
    window: Optional[int] = None,
    exponent_window: Optional[int] = None,
    base: Optional[int] = None,
    per_round: bool = False,
    threads: Optional[int] = None,
) -> ResourceCounts

Qubits and whole-run gate counts of the circuit, without simulating.

base defaults to the smallest base coprime to N (counts depend on it only weakly, through the constants of the lookup tables).

>>> c = resource_counts(1_005_973, "windowed-opt", base=980_062)
>>> c.qubits, c.oracle_gates, c.toffoli          # research/data/mbu-shor/counts.txt
(88, 271380, 70720)

oracle_circuit

oracle_circuit(
    N: int,
    a: int,
    oracle: str = "windowed-opt",
    *,
    window: Optional[int] = None,
) -> OracleCircuit

The controlled-multiplication-by-a block of oracle as a Circuit.

Supported: "windowed-opt", "windowed", "ripple" (X/CNOT/CCX only) and "beauregard" (QFT arithmetic). The MBU oracles need multi-bit classical feed-forward that a Circuit cannot express (UnsupportedOperationError); count them with resource_counts.

>>> o = oracle_circuit(15, 7, "ripple")
>>> o.num_qubits, o.work, o.toffolis > 0
(16, [1, 2, 3, 4], True)

shor_circuit

shor_circuit(
    N: int,
    a: int,
    oracle: str = "ripple",
    *,
    window: Optional[int] = None,
) -> Circuit

The whole semiclassical order-finding circuit as a Circuit.

Qubit 0 is the recycled control, 1 … n the work register (prepared in |1⟩). Round i: recycle (X conditioned on measurement i−1), H, controlled U^(2^(t−1−i)), one conditioned phase per earlier bit, H, measure. Measurement i is bit i of the measured integer, so simulate(c, samples(k)) samples the same distribution as factor.

>>> c = shor_circuit(15, 7, "ripple")
>>> c.num_qubits, c.num_measurements
(16, 8)

exact_distribution

exact_distribution(
    N: int,
    a: int,
    oracle: str = "permutation",
    *,
    window: Optional[int] = None,
    prune: float = 0.0,
    threads: Optional[int] = None,
) -> ndarray

Exact probabilities of the 2n-bit measured integer (float64[2**(2n)]), by walking the whole measurement tree of the semiclassical circuit (n ≤ 8 bits; 10 for "permutation", 6 for "beauregard").

>>> p = exact_distribution(15, 7)
>>> [int(y) for y in np.flatnonzero(p > 1e-12)], round(float(p[64]), 6)
([0, 64, 128, 192], 0.25)

predict_support

predict_support(
    N: int,
    a: int,
    oracle: str = "windowed-opt",
    *,
    window: Optional[int] = None,
    exponent_window: Optional[int] = None,
    precision: str = "f64",
) -> SupportPrediction

Predict the support trace, peak memory and work of one run, without running it.

>>> p = predict_support(1_005_973, 980_062)
>>> p.order, p.nu, p.peak, p.engine
(41832, 3, 20916, 'sliced')

multiplicative_order

multiplicative_order(a: int, N: int) -> int

The multiplicative order of a modulo N (classical; Pollard–Brent + Carmichael).

>>> multiplicative_order(7, 15), multiplicative_order(980_062, 1_005_973)
(4, 41832)

carmichael

carmichael(N: int) -> int

Carmichael's function λ(N) (the largest order of any unit mod N).

>>> carmichael(1_537_596_787)
256252500

noisy_success

noisy_success(
    N: int,
    a: int,
    p: float,
    *,
    noise: str = "depolarizing",
    trajectories: int = 200,
    oracle: str = "windowed",
    window: Optional[int] = None,
    faults: Optional[int] = None,
    seed: Optional[int] = None,
    reset_ancillas: bool = False,
    cap: int = 1 << 26,
    threads: Optional[int] = None,
) -> NoisySuccess

Success probability of the gate-level circuit under circuit-level Pauli noise.

Exact Monte-Carlo trajectories (research/shor/shor-noise.md): a Pauli fault at rate p after every oracle gate on each of its qubits, and on the control (preparation, after each H and the phase correction, readout flip). noise: "depolarizing", "bitflip" or "phaseflip". faults=k instead conditions every trajectory on exactly k faults at uniformly random locations (stratified estimates). reset_ancillas: an ideal reset of every ancilla after each round. Oracles: "windowed", "windowed-opt", "ripple" (≤ 129 qubits). Trajectory i uses a stream derived from (seed, i): results do not depend on threads.

>>> r = noisy_success(143, 2, 0.0, trajectories=20, seed=1)
>>> r.success, r.peak_rate, r.mean_faults
(1.0, 1.0, 0.0)