wai-quantum 0.6.4: Pauli frames, the GPU state vector, the density matrix to 16 qubits, an OpenQASM writer that keeps every gate; the service for agents and browsers #320
34 changed files with 7354 additions and 108 deletions
11
.github/workflows/wai-quantum-service.yml
vendored
11
.github/workflows/wai-quantum-service.yml
vendored
|
|
@ -101,6 +101,17 @@ jobs:
|
|||
set -euo pipefail
|
||||
nice -n 10 cargo test --locked --release 2>&1 | tee "$RUNNER_TEMP/test.log"
|
||||
grep -qE '^test result: ok\. [1-9][0-9]* passed; 0 failed' "$RUNNER_TEMP/test.log" || { echo "::error::no test ran"; exit 1; }
|
||||
# The single-precision state vector on this host's GPU, against the
|
||||
# double-precision one; the GPU is declared, so its absence fails.
|
||||
- name: The GPU state vector
|
||||
shell: bash
|
||||
working-directory: wai/wai-quantum
|
||||
env:
|
||||
WAI_QUANTUM_SV_METAL_EXPECT: "1"
|
||||
run: |
|
||||
set -euo pipefail
|
||||
nice -n 10 cargo test --locked --release --features quantum_sv_metal --lib quantum_sv_metal 2>&1 | tee "$RUNNER_TEMP/gpu.log"
|
||||
grep -qE '^test result: ok\. 4 passed; 0 failed' "$RUNNER_TEMP/gpu.log" || { echo "::error::not every GPU test ran and passed"; exit 1; }
|
||||
- name: The smoke checks over HTTP
|
||||
shell: bash
|
||||
working-directory: wai/wai-quantum-service
|
||||
|
|
|
|||
110
wai/CHANGELOG.md
110
wai/CHANGELOG.md
|
|
@ -6,6 +6,116 @@ wire format may change and such changes are recorded here.
|
|||
|
||||
## v0.4.0 — current
|
||||
|
||||
### wai-quantum 0.6.4: noisy Clifford circuits by Pauli frames, the state vector on the GPU, and the density matrix by superoperators
|
||||
|
||||
Additive. Every pinned digest is unchanged.
|
||||
|
||||
- **`DensityMatrix::evolve`** (`quantum_noisy`): the exact noisy state on up to 16 qubits
|
||||
(`MAX_EVOLVE_QUBITS`), against 12 for `DensityMatrix::run`, which stays as the referee.
|
||||
- Each gate goes on `ρ`'s rows, its conjugate on the columns. Each channel becomes one
|
||||
superoperator `Σᵢ wᵢ Kᵢ ⊗ Kᵢ*` on its qubits' row and column bits, in `run`'s order.
|
||||
- A gate and its noise fuse into one gate on at most four of the vector's qubits, applied
|
||||
in the state vector's tiled passes. A gate therefore costs one pass and no copy of `ρ`.
|
||||
- **Checked:** against `run`, every entry within 10⁻¹² on circuits with every channel and
|
||||
a Toffoli. At 14 qubits, past `run`'s size, against branching trajectories within four
|
||||
standard errors.
|
||||
- **Measured** on a 16-core desktop processor: 12 qubits in 0.53 s against 158 s for
|
||||
`run`, and 10 qubits in 39 ms against 4.1 s. 15 qubits with 132 noisy gates take 49 s.
|
||||
- **The service's `density_matrix`** runs on it, up to 16 qubits within a job's memory and
|
||||
a work cap. `auto`, with noise, still tries it first.
|
||||
- **`quantum_spd::expectation_capped`**: the sparse-Pauli engine held to a number of terms.
|
||||
Whenever a rotation leaves more, the threshold doubles and the sum is truncated again
|
||||
until it fits, and everything dropped that way joins the bound. A run that never reaches
|
||||
the cap is the uncapped run, bit for bit.
|
||||
- **Checked:** on the 127-qubit kicked-Ising circuit, a cap above the peak changes no bit.
|
||||
A cap a tenth of the peak stays within its bound of the uncapped value. On a small
|
||||
register, the bound covers the exact error.
|
||||
- **The service's `pauli`** holds as many terms as the job's memory allows (`max_terms`
|
||||
can lower it). A sum that would have outgrown the job now finishes at a coarser
|
||||
threshold, reported as `threshold_reached`, with its bound.
|
||||
- **The tensor-network executor blocked and threaded** (`quantum_tn`). Each pairwise
|
||||
contraction is a complex matrix product blocked for the cache, its rows shared between
|
||||
threads. Slices running in parallel keep their contractions on one thread each. Every
|
||||
output adds its terms in the same order, so the result is the same bits (a test holds it
|
||||
to the plain loop's, and the pinned digests hold).
|
||||
- **Measured:** large products went from 2.0 to about 25 billion complex multiply-adds a
|
||||
second on a 16-core desktop processor (1024 × 1024 × 1024: 0.52 s → 0.043 s).
|
||||
- **Faster state-vector kernels** (`quantum_sv`). Gates on two to five qubits run at
|
||||
compile-time sizes, several groups side by side, with groups of low strides gathered.
|
||||
Every output is the same sum in the same order, so the state is the same bits.
|
||||
- Measured in cache on one core of a desktop machine: two-qubit gates 280–345 → 400–460
|
||||
million amplitudes a second, and dense three- and four-qubit gates 10–25% faster.
|
||||
- On a 16-core desktop processor, dense three- to five-qubit gates are 35–50% faster.
|
||||
|
||||
- **`quantum_sv_metal`** (new feature, macOS): the dense state vector in single precision
|
||||
on the platform's GPU, in one buffer the GPU and the host share. Half the bytes of the
|
||||
double-precision state, so one more qubit fits.
|
||||
- Gates are fused as `quantum_sv` fuses them, then grouped into tile passes. Each
|
||||
threadgroup copies 4096 amplitudes into threadgroup memory, applies every gate of the
|
||||
pass there between barriers, and copies them back.
|
||||
- A kernel is built per group size, so each group stays in registers. The kernels are
|
||||
compiled with the platform's safe math mode. Norms, expectations and samples are summed
|
||||
on the host in double precision.
|
||||
- **Checked:**
|
||||
- Against the double-precision state vector on random circuits of 6 to 21 qubits at
|
||||
every fusion width: amplitudes within 2·10⁻⁶, expectations within 10⁻⁵.
|
||||
- Tiled and untiled runs are the same bits, and repeated runs are the same bits.
|
||||
- **Measured** on a desktop machine's GPU, with the machine heavily loaded by other work:
|
||||
a 30-qubit, depth-20 brickwork of 910 gates runs in about 4.2 s at fusion width 3 with
|
||||
tiles, against about 8 s gate by gate.
|
||||
- **The service** offers it as `statevector_single`, accuracy class `single`, where a GPU
|
||||
is present. `auto` tries it after the double-precision `statevector`, so it takes the
|
||||
circuits that one cannot hold: up to 33 qubits within a job's memory.
|
||||
- **Measured** through the public service, on a host under heavy load:
|
||||
- a 33-qubit job (a 64 GiB state, 97 gates, chosen by `auto`) in 70 s;
|
||||
- 31 qubits in 23.5 s;
|
||||
- 26 qubits in 2.1 s, against 10.1 s for the double-precision method. The expectation
|
||||
values agree to 10⁻⁸.
|
||||
|
||||
- **`quantum_tableau::frame_counts` and `frame_expectations`.** Clifford circuits under Pauli
|
||||
noise (`FrameStep`: depolarizing on one or two qubits, dephasing; readout flips), on up
|
||||
to 32768 qubits.
|
||||
- One noiseless outcome from the tableau is the reference. Each shot carries a Pauli
|
||||
frame, the difference its noise makes, 64 shots to a machine word.
|
||||
- Frames start with random Z bits. They leave `|0…0⟩` unchanged, and carried through the
|
||||
circuit they spread each shot uniformly over the noiseless outcomes.
|
||||
- Errors are placed by geometric gaps, with the crate's own logarithm.
|
||||
- An expectation value is the noiseless value times the frames' mean sign against the
|
||||
string, with its standard error.
|
||||
- Both readouts are the same bits at any thread count.
|
||||
- **Checked:** on random Clifford circuits of 3 to 5 qubits with every channel, the counts
|
||||
(200,000 shots, with readout flips) and the expectation values agree with the exact
|
||||
density matrix within five standard errors. Two mutations are caught: dropping one qubit
|
||||
of the two-qubit channel, and dropping the random Z bits.
|
||||
- **Measured** on a 16-core desktop processor: 100,000 full-width shots take 44 ms at 1,000
|
||||
qubits with 10,000 noisy gates, and 1.6 s at 10,000 qubits with 100,000.
|
||||
- **The service's `stabilizer` method takes Pauli noise.** Each QASM gate gets its noise on
|
||||
the qubits it touched, as the other noisy methods apply it. Expectations are estimated
|
||||
from at least 10,000 shots unless `trajectories` says otherwise. Amplitude damping is not
|
||||
a Pauli channel and is refused. With noise, `auto` tries `density_matrix` (12 qubits,
|
||||
exact), then `stabilizer`, then `trajectories`.
|
||||
- **The OpenQASM writer keeps every gate** (`quantum_toolchain::to_qasm`). A gate with
|
||||
three or more controls, such as Grover's oracle and diffusion past three qubits, was
|
||||
written as a comment, so the program written out was not the circuit. Every controlled
|
||||
gate without a standard name is now written in OpenQASM 3's modifier form,
|
||||
`ctrl(n) @ g`, and `from_qasm` reads `ctrl`, `ctrl(n)`, `negctrl` and `negctrl(n)` back.
|
||||
- **Checked:** each base gate under zero to five controls is read back op for op, and
|
||||
Grover's search on five and seven qubits keeps its state.
|
||||
- **The service** (`wai-quantum-service`, not published):
|
||||
- It reads the control modifiers on `x y z h s sdg t tdg p` and on `cx cy cz ch cp ccx
|
||||
ccz`. A gate with a standard name is read as that gate, so a Clifford program keeps
|
||||
the stabilizer method. A phase on m qubits is written exactly with 2^m − 1 `cx` and
|
||||
as many `p`, its angles dyadic when the phase is. **Checked:** on a 6-qubit state,
|
||||
every amplitude against the gate applied directly where the controls hold, within
|
||||
10⁻¹², each test failing when the expansion's signs or its releases are broken.
|
||||
- Receipts sign the RFC 8785 canonical text of the request and the result, so a result
|
||||
verifies after any JSON reader has carried it, a browser's included. Its numbers are
|
||||
checked against a JavaScript engine on RFC 8785's appendix B and on 2 million random
|
||||
doubles. Receipts made before it still verify.
|
||||
- It reads numbers correctly rounded (serde_json's `float_roundtrip`). Before, a result
|
||||
could change in its last bit on its way from the job's process to the reply.
|
||||
- A job carries its request as submitted, so whoever holds its id can check its receipt.
|
||||
|
||||
### `wai.image.jpegai`'s hyper decoders, context stages and synthesis through the integer operator set
|
||||
|
||||
- **§13.6.2 is normative.** [integer-operators](extensions/integer-operators.md) §13.6.2
|
||||
|
|
|
|||
2
wai/tools/quantum-meter/Cargo.lock
generated
2
wai/tools/quantum-meter/Cargo.lock
generated
|
|
@ -1017,7 +1017,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/tools/quantum-wgsl-check/Cargo.lock
generated
2
wai/tools/quantum-wgsl-check/Cargo.lock
generated
|
|
@ -599,7 +599,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-pack/Cargo.lock
generated
2
wai/wai-pack/Cargo.lock
generated
|
|
@ -713,7 +713,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-quantum-qpu/Cargo.lock
generated
2
wai/wai-quantum-qpu/Cargo.lock
generated
|
|
@ -422,7 +422,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-quantum-serve/Cargo.lock
generated
2
wai/wai-quantum-serve/Cargo.lock
generated
|
|
@ -599,7 +599,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
47
wai/wai-quantum-service/Cargo.lock
generated
47
wai/wai-quantum-service/Cargo.lock
generated
|
|
@ -72,6 +72,12 @@ version = "1.8.3"
|
|||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "2af50177e190e07a26ab74f8b1efbfe2ef87da2116221318cb1c2e82baf7de06"
|
||||
|
||||
[[package]]
|
||||
name = "bitflags"
|
||||
version = "2.13.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06"
|
||||
|
||||
[[package]]
|
||||
name = "blake3"
|
||||
version = "1.8.7"
|
||||
|
|
@ -468,6 +474,42 @@ dependencies = [
|
|||
"windows-sys",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "objc2"
|
||||
version = "0.6.5"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "08849bbd4767dfae9457696856ae1c84fe4e0281bbe4a7abff2d0e06fb7981f8"
|
||||
dependencies = [
|
||||
"objc2-encode",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "objc2-encode"
|
||||
version = "4.1.0"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "ef25abbcd74fb2609453eb695bd2f860d389e457f67dc17cafc8b8cbc89d0c33"
|
||||
|
||||
[[package]]
|
||||
name = "objc2-foundation"
|
||||
version = "0.3.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "e3e0adef53c21f888deb4fa59fc59f7eb17404926ee8a6f59f5df0fd7f9f3272"
|
||||
dependencies = [
|
||||
"bitflags",
|
||||
"objc2",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "objc2-metal"
|
||||
version = "0.3.2"
|
||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
||||
checksum = "a0125f776a10d00af4152d74616409f0d4a2053a6f57fa5b7d6aa2854ac04794"
|
||||
dependencies = [
|
||||
"bitflags",
|
||||
"objc2",
|
||||
"objc2-foundation",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "once_cell"
|
||||
version = "1.21.4"
|
||||
|
|
@ -815,10 +857,13 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
"objc2",
|
||||
"objc2-foundation",
|
||||
"objc2-metal",
|
||||
"serde_json",
|
||||
]
|
||||
|
||||
|
|
|
|||
|
|
@ -12,10 +12,12 @@ publish = false
|
|||
axum = "0.8"
|
||||
tokio = { version = "1", features = ["rt-multi-thread", "macros", "process", "time", "sync", "net", "io-util"] }
|
||||
serde = { version = "1", features = ["derive"] }
|
||||
serde_json = "1"
|
||||
# float_roundtrip: a number a client sends back reads as the double it was
|
||||
# (the default reader rounds twice), so its receipt still verifies.
|
||||
serde_json = { version = "1", features = ["float_roundtrip"] }
|
||||
blake3 = "1"
|
||||
getrandom = "0.4"
|
||||
wai-quantum = { path = "../wai-quantum", default-features = false, features = ["quantum", "quantum_toolchain", "quantum_receipt", "quantum_sv", "quantum_tn", "quantum_tableau", "quantum_mpsim", "quantum_spd", "quantum_noisy"] }
|
||||
wai-quantum = { path = "../wai-quantum", default-features = false, features = ["quantum", "quantum_toolchain", "quantum_receipt", "quantum_sv", "quantum_tn", "quantum_tableau", "quantum_mpsim", "quantum_spd", "quantum_noisy", "quantum_sv_metal"] }
|
||||
|
||||
[profile.release]
|
||||
lto = "thin"
|
||||
|
|
|
|||
|
|
@ -28,31 +28,36 @@ A job is JSON:
|
|||
|
||||
Counts are keyed by the program's classical bits, the highest bit first. With no measurement in the program, every qubit is measured into the bit of its own index.
|
||||
|
||||
The program may use the standard gates, user `gate` definitions, and the control modifiers `ctrl @`, `ctrl(n) @`, `negctrl @` and `negctrl(n) @` (on `x y z h s sdg t tdg p` and `cx cy cz ch cp ccx ccz`). Classical control, `inv @`, `pow @` and gates after a qubit's measurement are refused, each with its line and what to write instead (`src/qasm.rs` lists everything read).
|
||||
|
||||
## Methods
|
||||
|
||||
| method | takes | answers | accuracy |
|
||||
|---|---|---|---|
|
||||
| `exact` | angles that are dyadic multiples of 2π, up to 26 qubits | counts, expectations, amplitudes, a state hash | bit-exact |
|
||||
| `statevector` | any circuit, up to 32 qubits | counts, expectations, amplitudes | exact |
|
||||
| `stabilizer` | Clifford circuits (rotations by whole quarter turns), up to 32768 qubits | counts, expectations | exact |
|
||||
| `statevector_single` | any circuit with gates on at most four qubits, up to 33 qubits, where the host has a GPU | counts, expectations, amplitudes | single |
|
||||
| `stabilizer` | Clifford circuits (rotations by whole quarter turns), up to 32768 qubits, optionally with depolarizing, dephasing and readout noise | counts, expectations | exact; with noise, an estimate |
|
||||
| `pauli` | any circuit, up to 1024 qubits | expectations | bounded |
|
||||
| `mps` | any circuit, up to 4096 qubits, at bond dimension up to 256; best when entanglement stays low | counts, expectations, amplitudes | exact when nothing is truncated, otherwise heuristic |
|
||||
| `tensor_network` | any circuit | amplitudes | exact |
|
||||
| `density_matrix` | noisy circuits, up to 12 qubits | counts, expectations | exact |
|
||||
| `density_matrix` | noisy circuits, up to 16 qubits within a job's memory | counts, expectations | exact |
|
||||
| `trajectories` | noisy circuits, up to 32 qubits | counts, expectations | estimate |
|
||||
|
||||
How the methods run:
|
||||
|
||||
- **statevector** groups consecutive gates into passes over the state. Each pass applies its gates tile by tile, in a buffer that fits a core's cache. The grouping never changes a bit of the result.
|
||||
- **pauli** reports each value with its truncation `bound`. Unless `ladder` is `false`, it also reports a `ladder`: the value and bound at up to three thresholds a decade above the one asked for, so you can see the value settle as the threshold drops.
|
||||
- **pauli** reports each value with its truncation `bound`. It holds at most as many Pauli terms as the job's memory allows; a sum that would outgrow it raises its threshold until it fits (`threshold_reached`), and the bound covers everything that dropped. Unless `ladder` is `false`, it also reports a `ladder`: the value and bound at up to three thresholds a decade above the one asked for, so you can see the value settle as the threshold drops.
|
||||
- **stabilizer**, with noise, samples by Pauli frames against one noiseless outcome of the tableau, 64 shots to a machine word. Its counts have the law of independent noisy shots. Its expectation values are estimates with standard errors, from at least 10,000 shots unless `trajectories` sets the number. Amplitude damping is not a Pauli channel and is refused.
|
||||
- **trajectories** treats every shot as an independent noisy trajectory in law. Shots that took the same errors share one state, so the error-free part of a circuit, and the error-free majority of the shots, run once. `histories` in the result counts the distinct error histories that ran.
|
||||
|
||||
`auto` picks the first method that takes the job. With noise, it tries `density_matrix` then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it.
|
||||
`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it.
|
||||
|
||||
The accuracy classes:
|
||||
|
||||
- **bit-exact**: the same bits on every machine. Amplitudes are held in 30-bit fixed point, and every output is an exact function of them.
|
||||
- **exact**: no approximation beyond double-precision rounding.
|
||||
- **single**: no approximation beyond single-precision rounding, about 10⁻⁷ relative per gate; sums are taken in double precision.
|
||||
- **bounded**: an approximation, with a rigorous bound on the error (`bound`) beside each value.
|
||||
- **estimate**: a statistical estimate, with its standard error (`stderr`).
|
||||
- **heuristic**: an approximation, reported with the standard estimate of its error but not a bound. A matrix product state that truncated reports its discarded weight and its fidelity estimate, and a `warning` when the estimate is below 0.9.
|
||||
|
|
@ -66,7 +71,7 @@ Counts are always seeded samples.
|
|||
| `GET /v1` | the service: its methods, limits, load and receipt key |
|
||||
| `POST /v1/run` | submit a job and wait up to 25 s: the finished job, or `202` and the job if it is still running |
|
||||
| `POST /v1/jobs` | submit a job: `202` and the job, with its `id` |
|
||||
| `GET /v1/jobs/{id}` | the job: its status, and its result when it is done |
|
||||
| `GET /v1/jobs/{id}` | the job: its status, the request as submitted, and its result when it is done |
|
||||
| `DELETE /v1/jobs/{id}` | cancel a queued or running job |
|
||||
| `POST /v1/verify` | `{"request": ..., "result": ...}`: whether the result's receipt seals that request and that result |
|
||||
|
||||
|
|
@ -74,7 +79,12 @@ Every refusal says why: `400` for a program or job that does not read, `422` for
|
|||
|
||||
## Receipts
|
||||
|
||||
Every finished result carries a `receipt`: an Ed25519 signature, by the key `GET /v1` publishes, over the BLAKE3 digest of the request as submitted and the result as returned. `POST /v1/verify` checks one, and so can any Ed25519 implementation (see `src/receipt.rs`).
|
||||
Every finished result carries a `receipt`: an Ed25519 signature, by the key `GET /v1` publishes, over the BLAKE3 digest of the request as submitted and the result as returned. `POST /v1/verify` checks one, and so can any Ed25519 implementation:
|
||||
|
||||
- digest = BLAKE3(`wai:quantum-service-receipt` ‖ 0x02 ‖ len(R) ‖ R ‖ len(S) ‖ S), each len a big-endian u64;
|
||||
- R and S are the RFC 8785 canonical texts (`"canon": "rfc8785"`) of the request and of the result without its `receipt`: keys sorted, no whitespace, every number written as ECMAScript writes the double it reads.
|
||||
|
||||
So a result verifies after any JSON reader has carried it, a browser's included. Whole numbers in a request must be within ±(2^53 − 1), so that every reader reads them alike. Receipts issued before the canonical form carry no `canon` and still verify (see `src/receipt.rs`).
|
||||
|
||||
## Running it
|
||||
|
||||
|
|
|
|||
|
|
@ -49,6 +49,24 @@ json.dump({"request": req, "result": res}, open(sys.argv[4],"w"))' "$tmp/ghz.jso
|
|||
post /v1/verify "$tmp/v.json" | py 'import json,sys; d=json.load(sys.stdin); assert d["valid"] is True, d; print("receipt: verifies")'
|
||||
post /v1/verify "$tmp/t.json" | py 'import json,sys; d=json.load(sys.stdin); assert d["valid"] is False, d; print("receipt: a changed result does not verify")'
|
||||
|
||||
# The job carries its request as submitted, so whoever holds only the id can
|
||||
# check the receipt: verify from the service's answer alone.
|
||||
py 'import json,sys
|
||||
req=json.load(open(sys.argv[1])); job=json.load(open(sys.argv[2]))
|
||||
assert job["request"]==req, job.get("request")
|
||||
json.dump({"request": job["request"], "result": job["result"]}, open(sys.argv[3],"w"))' "$tmp/ghz.json" "$tmp/ghz.out" "$tmp/vj.json"
|
||||
post /v1/verify "$tmp/vj.json" | py 'import json,sys; d=json.load(sys.stdin); assert d["valid"] is True, d; print("receipt: verifies from the job alone")'
|
||||
|
||||
# As a browser sends it back: whole doubles without a fraction, keys reversed.
|
||||
py 'import json,sys
|
||||
def js(v):
|
||||
if isinstance(v, float) and v.is_integer(): return int(v)
|
||||
if isinstance(v, dict): return {k: js(v[k]) for k in reversed(list(v))}
|
||||
if isinstance(v, list): return [js(x) for x in v]
|
||||
return v
|
||||
v=json.load(open(sys.argv[1])); json.dump(js(v), open(sys.argv[2],"w"))' "$tmp/vj.json" "$tmp/vjs.json"
|
||||
post /v1/verify "$tmp/vjs.json" | py 'import json,sys; d=json.load(sys.stdin); assert d["valid"] is True, d; print("receipt: verifies after a browser-style round trip")'
|
||||
|
||||
# Continuous angles on a small circuit: the dense state vector, and the
|
||||
# probabilities sum to the shots.
|
||||
py 'import json; print(json.dumps({"qasm": "qubit[3] q; rx(0.3) q[0]; ry(1.2) q[1]; cx q[0], q[2]; crz(0.7) q[1], q[2];", "shots": 500, "observables": ["Z0", "Z1 Z2"]}))' > "$tmp/rot.json"
|
||||
|
|
@ -96,6 +114,91 @@ assert sum(r["counts"].values())==2000, r
|
|||
h=r["accuracy"]["histories"]; assert 1 <= h < 2000, r["accuracy"]
|
||||
print("noisy-14: ok,", r["method"], h, "histories for 2000 shots")' "$tmp/noisy.out"
|
||||
|
||||
# A wide noisy Clifford circuit: Pauli frames on the stabilizer tableau.
|
||||
py 'import json
|
||||
n=500
|
||||
q="qubit[%d] q;\n" % n + "h q[0];\n" + "".join("cx q[%d], q[%d];\n" % (i,i+1) for i in range(n-1))
|
||||
print(json.dumps({"qasm": q, "shots": 1000, "noise": {"depolarizing_2q": 0.0005, "readout": 0.001}, "observables": ["Z0 Z499"]}))' > "$tmp/nghz.json"
|
||||
post /v1/run "$tmp/nghz.json" > "$tmp/nghz.out"
|
||||
py 'import json,sys
|
||||
d=json.load(open(sys.argv[1])); assert d.get("status")=="done", d; r=d["result"]
|
||||
assert r["method"]=="stabilizer" and r["accuracy"]["class"]=="estimate", d
|
||||
assert sum(r["counts"].values())==1000, r
|
||||
e=r["expectations"][0]; assert 0.5 < e["value"] < 1.0 and e["stderr"] > 0, e
|
||||
print("noisy-ghz-500: ok,", r["method"], "<Z0 Z499> =", round(e["value"],4), "±", round(e["stderr"],4))' "$tmp/nghz.out"
|
||||
|
||||
# The single-precision state vector on the GPU, where the service has one,
|
||||
# against the double-precision one.
|
||||
if py 'import json,sys; d=json.load(open(sys.argv[1])); sys.exit(0 if any(m["method"]=="statevector_single" and m.get("available") for m in d["methods"]) else 1)' "$tmp/describe.json"; then
|
||||
py 'import json; print(json.dumps({"qasm": "qubit[12] q; h q; rx(0.3) q[0]; cx q[0], q[11]; crz(0.7) q[3], q[9]; ry(1.1) q[5];", "observables": ["Z0 Z11", "X5"], "method": "statevector_single"}))' > "$tmp/gpu.json"
|
||||
py 'import json,sys; d=json.load(open(sys.argv[1])); d["method"]="statevector"; json.dump(d, open(sys.argv[2],"w"))' "$tmp/gpu.json" "$tmp/cpu.json"
|
||||
post /v1/run "$tmp/gpu.json" > "$tmp/gpu.out"
|
||||
post /v1/run "$tmp/cpu.json" > "$tmp/cpu.out"
|
||||
py 'import json,sys
|
||||
g=json.load(open(sys.argv[1])); c=json.load(open(sys.argv[2])); assert g.get("status")=="done" and c.get("status")=="done", (g, c)
|
||||
gv=[e["value"] for e in g["result"]["expectations"]]; cv=[e["value"] for e in c["result"]["expectations"]]
|
||||
assert g["result"]["accuracy"]["class"]=="single" and all(abs(a-b)<1e-5 for a,b in zip(gv,cv)), (gv, cv)
|
||||
print("gpu: ok,", g["result"]["accuracy"]["device"], gv)' "$tmp/gpu.out" "$tmp/cpu.out"
|
||||
else
|
||||
echo "gpu: this instance has none"
|
||||
fi
|
||||
|
||||
# What a program finds first: the index, the examples, the OpenAPI document.
|
||||
curl -sS --max-time 30 -H 'accept: application/json' "$U/" > "$tmp/index.json"
|
||||
py 'import json,sys; d=json.load(open(sys.argv[1])); assert d["mcp"].endswith("/mcp") and d["openapi"].endswith("/v1/openapi.json"), d; print("index: ok,", d["mcp"])' "$tmp/index.json"
|
||||
curl -sS --max-time 30 "$U/v1/examples" > "$tmp/examples.json"
|
||||
curl -sS --max-time 30 "$U/v1/openapi.json" > "$tmp/openapi.json"
|
||||
py 'import json,sys
|
||||
e=json.load(open(sys.argv[1]))["examples"]; o=json.load(open(sys.argv[2]))
|
||||
assert len(e)>=5 and all("job" in x and "qasm" in x["job"] for x in e), e
|
||||
assert o["openapi"].startswith("3.1") and "/v1/run" in o["paths"] and "/mcp" in o["paths"], list(o["paths"])
|
||||
print("examples and openapi: ok,", len(e), "examples,", len(o["paths"]), "paths")' "$tmp/examples.json" "$tmp/openapi.json"
|
||||
|
||||
# An agent over the Model Context Protocol: initialize, list the tools, run
|
||||
# the first example, and check its receipt.
|
||||
mcp() { curl -sS --max-time 120 -X POST "$U/mcp" -H 'content-type: application/json' -H 'accept: application/json, text/event-stream' --data-binary "$1"; }
|
||||
mcp '{"jsonrpc":"2.0","id":1,"method":"initialize","params":{"protocolVersion":"2025-06-18","capabilities":{},"clientInfo":{"name":"smoke","version":"1"}}}' > "$tmp/mcp-init.json"
|
||||
code=$(curl -sS -o /dev/null -w '%{http_code}' --max-time 30 -X POST "$U/mcp" -H 'content-type: application/json' --data-binary '{"jsonrpc":"2.0","method":"notifications/initialized"}')
|
||||
mcp '{"jsonrpc":"2.0","id":2,"method":"tools/list"}' > "$tmp/mcp-tools.json"
|
||||
py 'import json,sys
|
||||
i=json.load(open(sys.argv[1]))["result"]; t=json.load(open(sys.argv[2]))["result"]["tools"]
|
||||
assert i["protocolVersion"]=="2025-06-18" and "tools" in i["capabilities"], i
|
||||
names={x["name"] for x in t}; assert {"run_quantum_circuit","get_quantum_job","verify_quantum_receipt","quantum_service_info","quantum_circuit_examples"} <= names, names
|
||||
assert sys.argv[3]=="202", sys.argv[3]
|
||||
print("mcp: ok, initialized,", len(t), "tools")' "$tmp/mcp-init.json" "$tmp/mcp-tools.json" "$code"
|
||||
py 'import json,sys
|
||||
job=json.load(open(sys.argv[1]))["examples"][0]["job"]
|
||||
print(json.dumps({"jsonrpc":"2.0","id":3,"method":"tools/call","params":{"name":"run_quantum_circuit","arguments":job}}))' "$tmp/examples.json" > "$tmp/mcp-call.json"
|
||||
mcp "@$tmp/mcp-call.json" > "$tmp/mcp-run.json"
|
||||
py 'import json,sys
|
||||
r=json.load(open(sys.argv[1]))["result"]; job=r["structuredContent"]; text=r["content"][0]["text"]
|
||||
assert not r["isError"] and job["status"]=="done" and "counts" in job["result"], r
|
||||
assert "method:" in text and "counts" in text, text
|
||||
req=json.load(open(sys.argv[2]))["examples"][0]["job"]
|
||||
print(json.dumps({"jsonrpc":"2.0","id":4,"method":"tools/call","params":{"name":"verify_quantum_receipt","arguments":{"request":req,"result":job["result"]}}}))' "$tmp/mcp-run.json" "$tmp/examples.json" > "$tmp/mcp-verify.json"
|
||||
curl -sS --max-time 60 -X POST "$U/mcp" -H 'content-type: application/json' --data-binary "@$tmp/mcp-verify.json" > "$tmp/mcp-verified.json"
|
||||
py 'import json,sys
|
||||
v=json.load(open(sys.argv[1]))["result"]; assert v["structuredContent"]["valid"] is True, v
|
||||
print("mcp run + verify: ok,", json.load(open(sys.argv[2]))["result"]["content"][0]["text"].splitlines()[1])' "$tmp/mcp-verified.json" "$tmp/mcp-run.json"
|
||||
|
||||
# Grover's search on five qubits as the IDE writes it: its oracle and its
|
||||
# diffusion are each a Z under four controls. Four rounds find the item.
|
||||
py 'import json
|
||||
n, marked = 5, 0b10110
|
||||
def mcz(): return "ctrl(4) @ z q[0],q[1],q[2],q[3],q[4];\n"
|
||||
q = "OPENQASM 3.0;\ninclude \"stdgates.inc\";\nqubit[5] q;\n" + "".join("h q[%d];\n" % i for i in range(n))
|
||||
for _ in range(4):
|
||||
flips = "".join("x q[%d];\n" % i for i in range(n) if not marked >> i & 1)
|
||||
q += flips + mcz() + flips
|
||||
q += "".join("h q[%d];\nx q[%d];\n" % (i, i) for i in range(n)) + mcz() + "".join("x q[%d];\nh q[%d];\n" % (i, i) for i in range(n))
|
||||
print(json.dumps({"qasm": q, "shots": 1000}))' > "$tmp/grover.json"
|
||||
post /v1/run "$tmp/grover.json" > "$tmp/grover.out"
|
||||
py 'import json,sys
|
||||
d=json.load(open(sys.argv[1])); r=d["result"]
|
||||
hits=r["counts"].get("10110", 0)
|
||||
assert d["status"]=="done" and hits > 980, d
|
||||
print("grover-5 (ctrl modifiers): ok,", r["method"], hits, "of 1000 on the marked item")' "$tmp/grover.out"
|
||||
|
||||
# A refusal says why.
|
||||
py 'import json; print(json.dumps({"qasm": "qubit[2] q; if (c) x q;"}))' > "$tmp/bad.json"
|
||||
code=$(curl -sS -o "$tmp/bad.out" -w '%{http_code}' -X POST "$U/v1/jobs" -H 'content-type: application/json' --data-binary "@$tmp/bad.json")
|
||||
|
|
|
|||
|
|
@ -9,17 +9,19 @@
|
|||
//! |---|---|---|---|
|
||||
//! | `exact` | dyadic angles, ≤ 26 qubits | counts, expectations, amplitudes, state hash | bit-exact |
|
||||
//! | `statevector` | any, up to the memory a job has (32 qubits) | counts, expectations, amplitudes | exact |
|
||||
//! | `stabilizer` | Clifford, ≤ 32768 qubits | counts, expectations | exact |
|
||||
//! | `statevector_single` | gates on ≤ 4 qubits, ≤ 33 qubits, on a GPU | counts, expectations, amplitudes | single |
|
||||
//! | `stabilizer` | Clifford, ≤ 32768 qubits, optionally with Pauli noise | counts, expectations | exact; an estimate with noise |
|
||||
//! | `pauli` | any, ≤ 1024 qubits | expectations | bounded |
|
||||
//! | `mps` | any, ≤ 4096 qubits, best on low-entanglement circuits | counts, expectations, amplitudes | exact, or heuristic with its fidelity estimate |
|
||||
//! | `tensor_network` | any | amplitudes | exact |
|
||||
//! | `density_matrix` | any, noisy, ≤ 12 qubits | counts, expectations | exact |
|
||||
//! | `density_matrix` | any, noisy, ≤ 16 qubits in a job's memory | counts, expectations | exact |
|
||||
//! | `trajectories` | any, noisy, ≤ 32 qubits | counts, expectations | estimate |
|
||||
//!
|
||||
//! `auto` picks, in order: with noise, `density_matrix` then `trajectories`;
|
||||
//! for counts, `stabilizer`, `exact`, `statevector`, `mps`; for expectations
|
||||
//! only, `stabilizer`, `exact`, `statevector`, `pauli`, `mps`; for amplitudes
|
||||
//! only, `exact`, `statevector`, `tensor_network`, `mps`. A method that cannot take the circuit is
|
||||
//! `auto` picks, in order: with noise, `density_matrix`, `stabilizer`, then `trajectories`;
|
||||
//! for counts, `stabilizer`, `exact`, `statevector`, `statevector_single`,
|
||||
//! `mps`; for expectations only, `stabilizer`, `exact`, `statevector`,
|
||||
//! `statevector_single`, `pauli`, `mps`; for amplitudes only, `exact`,
|
||||
//! `statevector`, `statevector_single`, `tensor_network`, `mps`. A method that cannot take the circuit is
|
||||
//! skipped, and the plan says why the one chosen was chosen.
|
||||
|
||||
use crate::qasm::{self, Program};
|
||||
|
|
@ -35,7 +37,10 @@ pub const MAX_OUTCOMES: usize = 65_536;
|
|||
pub const EXACT_MAX_QUBITS: u32 = 26;
|
||||
pub const STABILIZER_MAX_QUBITS: u32 = wai_quantum::quantum_tableau::MAX_QUBITS;
|
||||
pub const PAULI_MAX_QUBITS: u32 = 1024;
|
||||
pub const DENSITY_MAX_QUBITS: u32 = 12;
|
||||
pub const DENSITY_MAX_QUBITS: u32 = wai_quantum::quantum_noisy::MAX_EVOLVE_QUBITS;
|
||||
/// The most a density-matrix job may do: its gates times the `4ⁿ` entries
|
||||
/// of `ρ`, a few minutes of passes.
|
||||
pub const DENSITY_MAX_WORK: u64 = 1_000_000_000_000;
|
||||
pub const TRAJECTORY_MAX_QUBITS: u32 = 32;
|
||||
pub const MPS_MAX_QUBITS: u32 = wai_quantum::quantum_mpsim::MAX_QUBITS;
|
||||
/// The bond dimension an `mps` job gets unless it asks, and the most it may.
|
||||
|
|
@ -80,11 +85,12 @@ pub fn describe(job_bytes: u64) -> Value {
|
|||
json!([
|
||||
{ "method": "exact", "accuracy": "bit-exact", "circuits": "every angle a dyadic multiple of 2π (2πm/2^k, k ≤ 32)", "max_qubits": EXACT_MAX_QUBITS, "answers": ["counts", "expectations", "amplitudes", "state_hash"] },
|
||||
{ "method": "statevector", "accuracy": "exact", "circuits": "any", "max_qubits": statevector_max_qubits(job_bytes), "answers": ["counts", "expectations", "amplitudes"] },
|
||||
{ "method": "stabilizer", "accuracy": "exact", "circuits": "Clifford (rotations by whole quarter turns, controlled phases by whole half turns)", "max_qubits": STABILIZER_MAX_QUBITS, "answers": ["counts", "expectations"] },
|
||||
{ "method": "pauli", "accuracy": "bounded", "circuits": "any", "max_qubits": PAULI_MAX_QUBITS, "answers": ["expectations"], "options": { "threshold": "coefficients below it are dropped (default 1e-5); the bound covers what was dropped", "max_weight": "Pauli strings heavier than this are dropped", "ladder": "also report the value at up to three thresholds a decade apart above it, to show it settling (default true)" } },
|
||||
{ "method": "statevector_single", "accuracy": "single", "circuits": "any, gates on at most four qubits", "max_qubits": gpu_max_qubits().map(|m| m.min(single_max_qubits(job_bytes))).unwrap_or(0), "available": gpu_max_qubits().is_ok(), "answers": ["counts", "expectations", "amplitudes"], "options": { "fusion": "the widest fused gate (default 3, at most 3)" } },
|
||||
{ "method": "stabilizer", "accuracy": "exact; with Pauli noise, counts are noisy samples and expectations estimates", "circuits": "Clifford (rotations by whole quarter turns, controlled phases by whole half turns), optionally with depolarizing, dephasing and readout noise", "max_qubits": STABILIZER_MAX_QUBITS, "answers": ["counts", "expectations"], "options": { "trajectories": "with noise: the shots expectation values are estimated from (default: the shots, at least 10000)" } },
|
||||
{ "method": "pauli", "accuracy": "bounded", "circuits": "any", "max_qubits": PAULI_MAX_QUBITS, "answers": ["expectations"], "options": { "threshold": "coefficients below it are dropped (default 1e-5); the bound covers what was dropped", "max_weight": "Pauli strings heavier than this are dropped", "ladder": "also report the value at up to three thresholds a decade apart above it, to show it settling (default true)", "max_terms": "the most Pauli terms to hold (default: what the job's memory holds); past it the threshold rises, and the bound covers what drops" } },
|
||||
{ "method": "mps", "accuracy": "exact when nothing is truncated, otherwise heuristic", "circuits": "any; best on low-entanglement circuits", "max_qubits": MPS_MAX_QUBITS, "answers": ["counts", "expectations", "amplitudes"], "options": { "max_bond": "the largest bond dimension kept (default 64, at most 256)", "cutoff": "singular values whose squared share is below it are dropped (default 1e-12)" } },
|
||||
{ "method": "tensor_network", "accuracy": "exact", "circuits": "any", "max_qubits": qasm::MAX_QUBITS, "answers": ["amplitudes"] },
|
||||
{ "method": "density_matrix", "accuracy": "exact", "circuits": "any, with noise", "max_qubits": DENSITY_MAX_QUBITS, "answers": ["counts", "expectations"] },
|
||||
{ "method": "density_matrix", "accuracy": "exact", "circuits": "any, with noise", "max_qubits": density_max_qubits(job_bytes), "answers": ["counts", "expectations"] },
|
||||
{ "method": "trajectories", "accuracy": "estimate", "circuits": "any, with noise", "max_qubits": TRAJECTORY_MAX_QUBITS.min(statevector_max_qubits(job_bytes)), "answers": ["counts", "expectations"], "options": { "trajectories": "how many (default: the shots, or 1000 for expectations alone); shots that take the same errors share one state" } },
|
||||
])
|
||||
}
|
||||
|
|
@ -100,6 +106,23 @@ pub fn amplitude_bytes(n: u32) -> u64 {
|
|||
amps + amps / 32
|
||||
}
|
||||
|
||||
/// The most qubits the density matrix holds in `job_bytes`: `ρ` is a state
|
||||
/// vector on twice as many.
|
||||
pub fn density_max_qubits(job_bytes: u64) -> u32 {
|
||||
(1..=DENSITY_MAX_QUBITS).take_while(|&n| statevector_bytes(2 * n) <= job_bytes).last().unwrap_or(0)
|
||||
}
|
||||
|
||||
/// The most qubits the GPU's largest buffer holds, asked of the device once.
|
||||
pub fn gpu_max_qubits() -> Result<u32, String> {
|
||||
static MAX: std::sync::OnceLock<Result<u32, String>> = std::sync::OnceLock::new();
|
||||
MAX.get_or_init(|| wai_quantum::quantum_sv_metal::GpuState::max_qubits().map_err(|e| e.to_string())).clone()
|
||||
}
|
||||
|
||||
/// The most qubits the single-precision state vector holds in `job_bytes`.
|
||||
pub fn single_max_qubits(job_bytes: u64) -> u32 {
|
||||
(3..=40u32).take_while(|&n| (8u64 << n) + 256 * MIB <= job_bytes).last().unwrap_or(0)
|
||||
}
|
||||
|
||||
/// The most qubits the dense state vector holds in `job_bytes`.
|
||||
pub fn statevector_max_qubits(job_bytes: u64) -> u32 {
|
||||
(1..=wai_quantum::quantum_sv::MAX_QUBITS).rev().find(|&n| statevector_bytes(n) <= job_bytes).unwrap_or(0)
|
||||
|
|
@ -245,7 +268,7 @@ fn try_method(method: &str, r: &Request, dyadic: &Result<wai_quantum::quantum::C
|
|||
let counts = r.shots > 0;
|
||||
let wants_amplitudes = !r.amplitudes.is_empty();
|
||||
let wants_expectations = !r.observables.is_empty();
|
||||
if r.noise.is_some() && !matches!(method, "density_matrix" | "trajectories") {
|
||||
if r.noise.is_some() && !matches!(method, "density_matrix" | "trajectories" | "stabilizer") {
|
||||
return Err(format!("{method} does not model noise"));
|
||||
}
|
||||
match method {
|
||||
|
|
@ -266,6 +289,17 @@ fn try_method(method: &str, r: &Request, dyadic: &Result<wai_quantum::quantum::C
|
|||
}
|
||||
Ok(plan(bytes + r.shots * 24, gates << n, format!("dense double-precision state vector of {n} qubits")))
|
||||
}
|
||||
"statevector_single" => {
|
||||
let max = gpu_max_qubits().map_err(|e| format!("statevector_single needs a GPU: {e}"))?;
|
||||
let bytes = (8u64 << n.min(40)) + 256 * MIB;
|
||||
if n > max || bytes > cap {
|
||||
return Err(format!("statevector_single holds at most {} qubits on this GPU in a job's memory", max.min(single_max_qubits(cap))));
|
||||
}
|
||||
if r.program.ops.iter().any(|o| o.qubits.len() > wai_quantum::quantum_sv_metal::MAX_WIDTH) {
|
||||
return Err("statevector_single takes gates on at most four qubits".into());
|
||||
}
|
||||
Ok(plan(bytes + r.shots * 24, gates << n, format!("dense single-precision state vector of {n} qubits on the GPU")))
|
||||
}
|
||||
"stabilizer" => {
|
||||
if let Err(e) = clifford {
|
||||
return Err(format!("stabilizer needs a Clifford circuit: {e}"));
|
||||
|
|
@ -276,18 +310,35 @@ fn try_method(method: &str, r: &Request, dyadic: &Result<wai_quantum::quantum::C
|
|||
if wants_amplitudes {
|
||||
return Err("stabilizer answers counts and expectations".into());
|
||||
}
|
||||
if r.noise.is_some_and(|z| z.damping > 0.0) {
|
||||
return Err("stabilizer models Pauli noise (depolarizing, dephasing, readout), not amplitude damping".into());
|
||||
}
|
||||
let (nn, measured) = (n as u64, r.program.measured.len() as u64);
|
||||
let words = nn.div_ceil(64);
|
||||
// The gates; the elimination for counts; each shot; each observable.
|
||||
let work = gates * 2 * words
|
||||
let mut work = gates * 2 * words
|
||||
+ if counts { nn * nn * words + r.shots * nn * measured.div_ceil(64) } else { 0 }
|
||||
+ r.observables.len() as u64 * 2 * nn * words;
|
||||
if r.noise.is_some() {
|
||||
// Pauli frames: every step on a word of 64 shots, and each
|
||||
// observable's parity over the shots.
|
||||
let lowered = clifford.as_ref().map_or(0, |c| c.len() as u64);
|
||||
let steps = lowered + 3 * r.program.ops.len() as u64;
|
||||
let frame_shots = r.shots.max(if wants_expectations { frame_expectation_shots(r) } else { 0 });
|
||||
work += steps * frame_shots.div_ceil(64) + r.observables.len() as u64 * frame_shots.div_ceil(64) * 4;
|
||||
}
|
||||
if work > STABILIZER_MAX_WORK {
|
||||
return Err(format!("this Clifford job is more stabilizer work than one job may take ({work} word operations, at most {STABILIZER_MAX_WORK}); ask for fewer shots"));
|
||||
}
|
||||
// The tableau, its rows and the outcome subspace, about n² bytes
|
||||
// each; and the counts.
|
||||
let memory = 3 * nn * nn / 2 + 64 * MIB + r.shots.min(1 << 20) * (r.program.n_bits as u64 / 8 + 64);
|
||||
let mut memory = 3 * nn * nn / 2 + 64 * MIB + r.shots.min(1 << 20) * (r.program.n_bits as u64 / 8 + 64);
|
||||
if r.noise.is_some() {
|
||||
// Up to eight batches of frames at once, 8192 shots each, and
|
||||
// their outcomes.
|
||||
memory += 8 * (2 * nn * 128 * 8 + 8192 * (measured.div_ceil(64) * 8 + 64));
|
||||
return Ok(plan(memory, work, format!("noisy Clifford circuit: Pauli frames of {n} qubits against a stabilizer-tableau reference")));
|
||||
}
|
||||
Ok(plan(memory, work, format!("Clifford circuit: the stabilizer tableau of {n} qubits")))
|
||||
}
|
||||
"pauli" => {
|
||||
|
|
@ -321,10 +372,14 @@ fn try_method(method: &str, r: &Request, dyadic: &Result<wai_quantum::quantum::C
|
|||
if r.noise.is_none() {
|
||||
return Err("density_matrix is for noisy circuits (give \"noise\")".into());
|
||||
}
|
||||
if n > DENSITY_MAX_QUBITS || wants_amplitudes {
|
||||
return Err(format!("density_matrix holds at most {DENSITY_MAX_QUBITS} qubits and answers counts and expectations"));
|
||||
let most = density_max_qubits(cap);
|
||||
if n > most || wants_amplitudes {
|
||||
return Err(format!("density_matrix holds at most {most} qubits in a job's memory and answers counts and expectations"));
|
||||
}
|
||||
Ok(plan((32u64 << (2 * n)) + 128 * MIB, gates << (2 * n), format!("density matrix of {n} qubits under the noise model")))
|
||||
if (gates << (2 * n)) > DENSITY_MAX_WORK {
|
||||
return Err(format!("density_matrix on {n} qubits holds 4^{n} entries; {gates} gates of passes over them is more than one job may take"));
|
||||
}
|
||||
Ok(plan(statevector_bytes(2 * n), gates << (2 * n), format!("density matrix of {n} qubits under the noise model, each gate and its noise one superoperator pass")))
|
||||
}
|
||||
"trajectories" => {
|
||||
if r.noise.is_none() {
|
||||
|
|
@ -369,6 +424,13 @@ pub fn trajectory_states(n: u32, cap: u64) -> usize {
|
|||
((cap.saturating_sub(256 * MIB) / 2 / (amplitude_bytes(n) * trees).max(1)) as usize).clamp(1, 16)
|
||||
}
|
||||
|
||||
/// Shots a noisy stabilizer job estimates expectation values from: the
|
||||
/// `trajectories` option, or the job's shots, at least 10000 (frames are
|
||||
/// cheap, so the default buys a small standard error).
|
||||
pub fn frame_expectation_shots(r: &Request) -> u64 {
|
||||
r.options.get("trajectories").and_then(Value::as_u64).unwrap_or(r.shots.max(10_000)).clamp(1, 10_000_000)
|
||||
}
|
||||
|
||||
/// Trajectories: the option, or one a shot, or 1000 for expectations alone.
|
||||
pub fn trajectory_count(r: &Request) -> u64 {
|
||||
let default = if r.shots > 0 { r.shots } else { 1000 };
|
||||
|
|
@ -384,13 +446,13 @@ pub fn plan(v: &Value, job_bytes: u64) -> Result<Plan, String> {
|
|||
return try_method(&r.method, &r, &dyadic, &clifford, job_bytes);
|
||||
}
|
||||
let order: &[&str] = if r.noise.is_some() {
|
||||
&["density_matrix", "trajectories"]
|
||||
&["density_matrix", "stabilizer", "trajectories"]
|
||||
} else if r.shots > 0 {
|
||||
&["stabilizer", "exact", "statevector", "mps"]
|
||||
&["stabilizer", "exact", "statevector", "statevector_single", "mps"]
|
||||
} else if !r.observables.is_empty() && r.amplitudes.is_empty() {
|
||||
&["stabilizer", "exact", "statevector", "pauli", "mps"]
|
||||
&["stabilizer", "exact", "statevector", "statevector_single", "pauli", "mps"]
|
||||
} else {
|
||||
&["exact", "statevector", "tensor_network", "mps"]
|
||||
&["exact", "statevector", "statevector_single", "tensor_network", "mps"]
|
||||
};
|
||||
let mut refusals = Vec::new();
|
||||
for m in order {
|
||||
|
|
|
|||
222
wai/wai-quantum-service/src/discovery.rs
Normal file
222
wai/wai-quantum-service/src/discovery.rs
Normal file
|
|
@ -0,0 +1,222 @@
|
|||
//! What a newcomer, person or program, needs to start: worked examples, an
|
||||
//! OpenAPI description of the API, and an index at the root.
|
||||
|
||||
use serde_json::{Value, json};
|
||||
|
||||
/// A chain of `n` qubits: `body(i)` for each qubit, as OpenQASM lines.
|
||||
fn lines(n: usize, body: impl Fn(usize) -> String) -> String {
|
||||
(0..n).map(body).collect::<Vec<_>>().join("\n")
|
||||
}
|
||||
|
||||
/// Ready-to-run jobs, each with what it shows. Every one runs as given.
|
||||
pub fn examples() -> Value {
|
||||
let ghz = |n: usize| format!("OPENQASM 3.0;\nqubit[{n}] q;\nbit[{n}] c;\nh q[0];\n{}\nc = measure q;", lines(n - 1, |i| format!("cx q[{i}], q[{}];", i + 1)));
|
||||
let ring_qaoa = {
|
||||
let n = 6;
|
||||
let mut s = format!("OPENQASM 3.0;\nqubit[{n}] q;\nbit[{n}] c;\nh q;\n");
|
||||
s += &lines(n, |i| format!("rzz(0.8) q[{i}], q[{}];", (i + 1) % n));
|
||||
s += "\nrx(0.7) q;\nc = measure q;";
|
||||
s
|
||||
};
|
||||
let trotter = |n: usize, steps: usize| {
|
||||
let mut s = format!("OPENQASM 3.0;\nqubit[{n}] q;\n");
|
||||
for _ in 0..steps {
|
||||
s += &format!("rx(0.6) q;\n{}\n", lines(n - 1, |i| format!("rzz(-1.5707963267948966) q[{i}], q[{}];", i + 1)));
|
||||
}
|
||||
s
|
||||
};
|
||||
let brick = |n: usize, layers: usize| {
|
||||
let mut s = format!("OPENQASM 3.0;\nqubit[{n}] q;\nbit[{n}] c;\n");
|
||||
for l in 0..layers {
|
||||
s += &format!("ry(0.3) q;\nrz(0.2) q;\n{}\n", lines(n / 2 - usize::from(l % 2 == 1), |i| format!("cx q[{}], q[{}];", 2 * i + l % 2, 2 * i + l % 2 + 1)));
|
||||
}
|
||||
s += "c = measure q;";
|
||||
s
|
||||
};
|
||||
json!([
|
||||
{
|
||||
"name": "bell",
|
||||
"title": "A Bell pair",
|
||||
"shows": "Two qubits entangled: half the shots read 00, half 11, never 01 or 10. The service sees a Clifford circuit and uses the stabilizer tableau, exactly.",
|
||||
"job": { "qasm": ghz(2), "shots": 1000, "observables": ["Z0 Z1", "X0 X1"] }
|
||||
},
|
||||
{
|
||||
"name": "rotation",
|
||||
"title": "A rotation by any angle",
|
||||
"shows": "rx(0.3) on one qubit: ⟨Z⟩ = cos 0.3 = 0.9553, ⟨Y⟩ = −sin 0.3. A continuous angle takes the double-precision state vector.",
|
||||
"job": { "qasm": "OPENQASM 3.0;\nqubit[1] q;\nrx(0.3) q[0];", "observables": ["Z0", "Y0"], "amplitudes": ["0", "1"] }
|
||||
},
|
||||
{
|
||||
"name": "qaoa-ring",
|
||||
"title": "QAOA on a six-node ring",
|
||||
"shows": "One layer of the quantum approximate optimisation algorithm for max-cut on a ring: cost terms ZZ on each edge, then a mixer. The counts favour cuts; the ZZ expectations score them.",
|
||||
"job": { "qasm": ring_qaoa, "shots": 2000, "observables": ["Z0 Z1", "Z1 Z2", "Z2 Z3", "Z3 Z4", "Z4 Z5", "Z5 Z0"] }
|
||||
},
|
||||
{
|
||||
"name": "noisy-bell",
|
||||
"title": "A Bell pair on a noisy device",
|
||||
"shows": "The same pair with gate errors and readout flips: 01 and 10 appear, and ⟨Z0 Z1⟩ falls below 1. With noise on few qubits the exact density matrix runs.",
|
||||
"job": { "qasm": ghz(2), "shots": 4000, "observables": ["Z0 Z1"], "noise": { "depolarizing_1q": 0.01, "depolarizing_2q": 0.03, "readout": 0.02 } }
|
||||
},
|
||||
{
|
||||
"name": "ghz-1000",
|
||||
"title": "A thousand-qubit GHZ state",
|
||||
"shows": "Far past any state vector: a Clifford circuit on 1000 qubits, sampled exactly by the stabilizer tableau. Every shot is all zeros or all ones.",
|
||||
"job": { "qasm": ghz(1000), "shots": 100, "observables": ["Z0 Z999"] }
|
||||
},
|
||||
{
|
||||
"name": "noisy-ghz-500",
|
||||
"title": "A noisy 500-qubit GHZ state",
|
||||
"shows": "Pauli noise on a large Clifford circuit, sampled by Pauli frames: ⟨Z0 Z499⟩ decays with every gate between them.",
|
||||
"job": { "qasm": ghz(500), "shots": 1000, "observables": ["Z0 Z499"], "noise": { "depolarizing_2q": 0.0005, "readout": 0.001 } }
|
||||
},
|
||||
{
|
||||
"name": "chain-100-dynamics",
|
||||
"title": "Quench dynamics on 100 qubits",
|
||||
"shows": "Four Trotter steps of a transverse-field Ising chain. Expectation values only, so sparse Pauli dynamics runs: each value comes with a bound on its error, and a ladder of thresholds shows it settling.",
|
||||
"job": { "qasm": trotter(100, 4), "observables": ["Z50", "Z49 Z50"] }
|
||||
},
|
||||
{
|
||||
"name": "mps-200",
|
||||
"title": "200 qubits at low entanglement",
|
||||
"shows": "Shallow brickwork on 200 qubits with any angles: a matrix product state holds it, and says whether anything was truncated.",
|
||||
"job": { "qasm": brick(200, 4), "shots": 200, "observables": ["Z100"] }
|
||||
}
|
||||
])
|
||||
}
|
||||
|
||||
/// The API as an OpenAPI 3.1 document.
|
||||
pub fn openapi(server: &str) -> Value {
|
||||
let pauli = json!({ "oneOf": [
|
||||
{ "type": "string", "description": "A Pauli string: factors such as `Z0`, `X3`, `Y12`, separated by spaces; qubit indices count from 0.", "examples": ["Z0 Z1", "X3 Y4"] },
|
||||
{ "type": "object", "required": ["pauli"], "properties": { "pauli": { "type": "string" }, "coeff": { "type": "number", "default": 1.0 } } }
|
||||
]});
|
||||
let job_request = json!({
|
||||
"type": "object",
|
||||
"required": ["qasm"],
|
||||
"properties": {
|
||||
"qasm": { "type": "string", "description": "The program: OpenQASM 3, or OpenQASM 2. Registers, standard gates (h, x, y, z, s, sdg, t, tdg, sx, rx, ry, rz, p, u, cx, cy, cz, ch, swap, iswap, cp, crx, cry, crz, ccx, ccz, cswap, rxx, ryy, rzz), the control modifiers ctrl @, ctrl(n) @ and negctrl @ on x y z h s sdg t tdg p and the controlled names, user gate definitions and terminal measurement. Mid-circuit measurement, classical control, inv @ and pow @ are refused, with the line." },
|
||||
"method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "density_matrix", "trajectories"], "default": "auto" },
|
||||
"shots": { "type": "integer", "minimum": 0, "maximum": 1000000, "description": "Measurement shots (default 1024 when nothing else is asked for)." },
|
||||
"seed": { "type": "integer", "default": 1, "description": "The sampling seed: the same job and seed give the same counts." },
|
||||
"observables": { "type": "array", "items": pauli },
|
||||
"amplitudes": { "type": "array", "items": { "type": "string", "pattern": "^[01]+$" }, "description": "Bit strings, one character per qubit, qubit n−1 first." },
|
||||
"noise": { "type": "object", "properties": {
|
||||
"depolarizing_1q": { "type": "number", "minimum": 0, "maximum": 1 },
|
||||
"depolarizing_2q": { "type": "number", "minimum": 0, "maximum": 1 },
|
||||
"amplitude_damping": { "type": "number", "minimum": 0, "maximum": 1 },
|
||||
"dephasing": { "type": "number", "minimum": 0, "maximum": 1 },
|
||||
"readout": { "type": "number", "minimum": 0, "maximum": 1 }
|
||||
}, "description": "Each a probability, applied after every gate on the qubits it touched (readout: each measured bit)." },
|
||||
"options": { "type": "object", "description": "Per method: threshold, max_weight, ladder, max_terms (pauli); max_bond, cutoff (mps); fusion (statevector); trajectories; memory_gb." }
|
||||
}
|
||||
});
|
||||
let job = json!({
|
||||
"type": "object",
|
||||
"properties": {
|
||||
"id": { "type": "string" },
|
||||
"status": { "type": "string", "enum": ["queued", "running", "done", "failed", "cancelled"] },
|
||||
"method": { "type": "string" },
|
||||
"why": { "type": "string", "description": "Why this method: what it is, and what the circuit is." },
|
||||
"request": { "$ref": "#/components/schemas/JobRequest", "description": "The request as submitted: the receipt is checked against it." },
|
||||
"result": { "type": "object", "description": "counts (highest classical bit first), expectations, amplitudes, accuracy, elapsed_ms and a receipt." },
|
||||
"error": { "type": "string" }
|
||||
}
|
||||
});
|
||||
let err = json!({ "type": "object", "properties": { "error": { "type": "string" } } });
|
||||
json!({
|
||||
"openapi": "3.1.0",
|
||||
"info": {
|
||||
"title": "WAI quantum service",
|
||||
"version": env!("CARGO_PKG_VERSION"),
|
||||
"summary": "Run quantum programs: counts, expectation values and amplitudes, each with the method that produced it, its accuracy and a signed receipt. No account and no key.",
|
||||
"license": { "name": "Apache-2.0" }
|
||||
},
|
||||
"servers": [{ "url": server }],
|
||||
"components": { "schemas": { "JobRequest": job_request, "Job": job, "Error": err } },
|
||||
"paths": {
|
||||
"/v1": { "get": { "summary": "The service: methods, limits, load and receipt key", "responses": { "200": { "description": "The service" } } } },
|
||||
"/v1/methods": { "get": { "summary": "The methods and what each takes", "responses": { "200": { "description": "The methods" } } } },
|
||||
"/v1/examples": { "get": { "summary": "Ready-to-run example jobs", "responses": { "200": { "description": "The examples" } } } },
|
||||
"/v1/run": { "post": {
|
||||
"summary": "Submit a job and wait up to 25 s for it",
|
||||
"requestBody": { "required": true, "content": { "application/json": { "schema": { "$ref": "#/components/schemas/JobRequest" } } } },
|
||||
"responses": {
|
||||
"200": { "description": "The finished job", "content": { "application/json": { "schema": { "$ref": "#/components/schemas/Job" } } } },
|
||||
"202": { "description": "Still running: poll GET /v1/jobs/{id}", "content": { "application/json": { "schema": { "$ref": "#/components/schemas/Job" } } } },
|
||||
"400": { "description": "A program or job that does not read", "content": { "application/json": { "schema": { "$ref": "#/components/schemas/Error" } } } },
|
||||
"422": { "description": "Larger than one job may be" },
|
||||
"429": { "description": "Past the per-client limits" },
|
||||
"503": { "description": "The queue is full" }
|
||||
}
|
||||
} },
|
||||
"/v1/jobs": { "post": {
|
||||
"summary": "Submit a job",
|
||||
"requestBody": { "required": true, "content": { "application/json": { "schema": { "$ref": "#/components/schemas/JobRequest" } } } },
|
||||
"responses": { "202": { "description": "Accepted", "content": { "application/json": { "schema": { "$ref": "#/components/schemas/Job" } } } } }
|
||||
} },
|
||||
"/v1/jobs/{id}": {
|
||||
"parameters": [{ "name": "id", "in": "path", "required": true, "schema": { "type": "string" } }],
|
||||
"get": { "summary": "A job: its status, and its result when done", "responses": { "200": { "description": "The job", "content": { "application/json": { "schema": { "$ref": "#/components/schemas/Job" } } } }, "404": { "description": "No such job" } } },
|
||||
"delete": { "summary": "Cancel a queued or running job", "responses": { "200": { "description": "Cancelled" }, "409": { "description": "Already finished" } } }
|
||||
},
|
||||
"/v1/verify": { "post": {
|
||||
"summary": "Whether a result's receipt seals that request and that result",
|
||||
"requestBody": { "required": true, "content": { "application/json": { "schema": { "type": "object", "required": ["request", "result"], "properties": { "request": { "$ref": "#/components/schemas/JobRequest" }, "result": { "type": "object" } } } } } },
|
||||
"responses": { "200": { "description": "{\"valid\": bool}" } }
|
||||
} },
|
||||
"/mcp": { "post": { "summary": "The Model Context Protocol endpoint (streamable HTTP, JSON responses): tools to run circuits, read jobs, verify receipts and list examples", "responses": { "200": { "description": "A JSON-RPC response" } } } }
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
/// Where everything is, for a program that lands on the root.
|
||||
pub fn index(base: &str, console: &str) -> Value {
|
||||
json!({
|
||||
"service": "wai-quantum-service",
|
||||
"summary": "Run quantum programs (OpenQASM 3 or 2). No account and no key.",
|
||||
"console": console,
|
||||
"api": format!("{base}/v1"),
|
||||
"run": format!("POST {base}/v1/run"),
|
||||
"examples": format!("{base}/v1/examples"),
|
||||
"openapi": format!("{base}/v1/openapi.json"),
|
||||
"mcp": format!("{base}/mcp"),
|
||||
})
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
/// Every example plans, with the method its description promises.
|
||||
#[test]
|
||||
fn every_example_plans_as_described() {
|
||||
let want = [
|
||||
("bell", "stabilizer"),
|
||||
("rotation", "statevector"),
|
||||
("qaoa-ring", "statevector"),
|
||||
("noisy-bell", "density_matrix"),
|
||||
("ghz-1000", "stabilizer"),
|
||||
("noisy-ghz-500", "stabilizer"),
|
||||
("chain-100-dynamics", "pauli"),
|
||||
("mps-200", "mps"),
|
||||
];
|
||||
let ex = examples();
|
||||
let ex = ex.as_array().unwrap();
|
||||
assert_eq!(ex.len(), want.len());
|
||||
for (e, (name, method)) in ex.iter().zip(want) {
|
||||
assert_eq!(e["name"], name);
|
||||
let plan = crate::catalog::plan(&e["job"], 72 << 30).unwrap_or_else(|why| panic!("{name}: {why}"));
|
||||
assert_eq!(plan.method, method, "{name}: {}", plan.why);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_openapi_document_names_every_route() {
|
||||
let d = openapi("https://example.org/q");
|
||||
for p in ["/v1", "/v1/methods", "/v1/examples", "/v1/run", "/v1/jobs", "/v1/jobs/{id}", "/v1/verify", "/mcp"] {
|
||||
assert!(d["paths"].get(p).is_some(), "{p}");
|
||||
}
|
||||
assert_eq!(d["openapi"], "3.1.0");
|
||||
}
|
||||
}
|
||||
|
|
@ -10,6 +10,10 @@
|
|||
//! POST /v1/run submit and wait (up to 25 s): the finished job,
|
||||
//! or 202 and the job if it is still running
|
||||
//! POST /v1/verify a receipt checked against its request and result
|
||||
//! GET /v1/examples ready-to-run example jobs
|
||||
//! GET /v1/openapi.json this API as an OpenAPI 3.1 document
|
||||
//! POST /mcp the Model Context Protocol endpoint (see `mcp`)
|
||||
//! GET / an index (a browser is sent to the console page)
|
||||
//! ```
|
||||
//!
|
||||
//! A job is `{"qasm": "<OpenQASM 3 program>", "method": "auto", "shots": N,
|
||||
|
|
@ -40,6 +44,10 @@ pub struct App {
|
|||
pub signer: Arc<Signer>,
|
||||
/// Peers whose `X-Forwarded-For` names the client: the front proxy.
|
||||
pub trusted: Vec<IpAddr>,
|
||||
/// The service's public base URL, for links it hands out.
|
||||
pub base: String,
|
||||
/// The page people use the service from.
|
||||
pub console: String,
|
||||
}
|
||||
|
||||
type Shared = Arc<App>;
|
||||
|
|
@ -89,7 +97,9 @@ async fn describe(State(app): State<Shared>) -> Response {
|
|||
"load": { "queued": queued, "running": running, "memory_held_bytes": held << 20 },
|
||||
"uptime_seconds": app.service.started.elapsed().map(|d| d.as_secs()).unwrap_or(0),
|
||||
"receipt_key": app.signer.public_key(),
|
||||
"endpoints": ["GET /v1", "GET /v1/health", "GET /v1/methods", "POST /v1/jobs", "GET /v1/jobs/{id}", "DELETE /v1/jobs/{id}", "POST /v1/run", "POST /v1/verify"],
|
||||
"endpoints": ["GET /v1", "GET /v1/health", "GET /v1/methods", "GET /v1/examples", "GET /v1/openapi.json", "POST /v1/jobs", "GET /v1/jobs/{id}", "DELETE /v1/jobs/{id}", "POST /v1/run", "POST /v1/verify", "POST /mcp"],
|
||||
"mcp": format!("{}/mcp", app.base),
|
||||
"console": app.console,
|
||||
}),
|
||||
)
|
||||
}
|
||||
|
|
@ -104,7 +114,10 @@ async fn health() -> Response {
|
|||
|
||||
/// Plans and admits `request` for `who`: the job's id, or the refusal's
|
||||
/// status and reason.
|
||||
fn admit(app: &App, who: &str, request: Value) -> Result<String, (u16, String)> {
|
||||
pub(crate) fn admit(app: &App, who: &str, request: Value) -> Result<String, (u16, String)> {
|
||||
if let Some(at) = crate::receipt::inexact_number(&request) {
|
||||
return Err((400, format!("`{at}`: a whole number here must be within ±9007199254740991 (2^53 − 1), so that every JSON reader reads it alike and the receipt binds it")));
|
||||
}
|
||||
let plan = catalog::plan(&request, app.service.config.job_mib << 20).map_err(|e| (400, e))?;
|
||||
app.service.submit(who, request, plan).map_err(|Refused(code, why)| (code, why))
|
||||
}
|
||||
|
|
@ -169,12 +182,58 @@ async fn verify(body: Result<Json<Value>, axum::extract::rejection::JsonRejectio
|
|||
reply(200, json!({ "valid": receipt::verify(request, &bare, receipt), "key": receipt.get("key") }))
|
||||
}
|
||||
|
||||
async fn examples() -> Response {
|
||||
reply(200, json!({ "examples": crate::discovery::examples() }))
|
||||
}
|
||||
|
||||
async fn openapi(State(app): State<Shared>) -> Response {
|
||||
reply(200, crate::discovery::openapi(&app.base))
|
||||
}
|
||||
|
||||
/// The root: a browser goes to the console page; a program gets an index.
|
||||
async fn index(State(app): State<Shared>, headers: HeaderMap) -> Response {
|
||||
let html = headers.get("accept").and_then(|v| v.to_str().ok()).is_some_and(|a| a.contains("text/html"));
|
||||
if html && let Ok(to) = HeaderValue::from_str(&app.console) {
|
||||
let mut r = StatusCode::FOUND.into_response();
|
||||
r.headers_mut().insert("location", to);
|
||||
return r;
|
||||
}
|
||||
reply(200, crate::discovery::index(&app.base, &app.console))
|
||||
}
|
||||
|
||||
/// One JSON-RPC message (MCP streamable HTTP, answered as JSON): a request
|
||||
/// gets its response; a notification gets 202 and no body.
|
||||
async fn mcp(State(app): State<Shared>, ConnectInfo(peer): ConnectInfo<SocketAddr>, headers: HeaderMap, body: Result<Json<Value>, axum::extract::rejection::JsonRejection>) -> Response {
|
||||
let Ok(Json(msg)) = body else {
|
||||
return reply(400, json!({ "jsonrpc": "2.0", "id": null, "error": { "code": -32700, "message": "not JSON" } }));
|
||||
};
|
||||
if !msg.is_object() {
|
||||
return reply(400, json!({ "jsonrpc": "2.0", "id": null, "error": { "code": -32600, "message": "one JSON-RPC message a request" } }));
|
||||
}
|
||||
let who = client(&app, peer, &headers);
|
||||
match crate::mcp::handle(&app, &who, &msg, &app.base).await {
|
||||
Some(answer) => reply(200, answer),
|
||||
None => {
|
||||
let mut r = StatusCode::ACCEPTED.into_response();
|
||||
r.headers_mut().insert("access-control-allow-origin", HeaderValue::from_static("*"));
|
||||
r
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// No event stream is offered: `GET /mcp` is 405, as the protocol allows.
|
||||
async fn mcp_get() -> Response {
|
||||
let mut r = error(405, "this endpoint takes POST; it offers no event stream");
|
||||
r.headers_mut().insert("allow", HeaderValue::from_static("POST"));
|
||||
r
|
||||
}
|
||||
|
||||
async fn preflight() -> Response {
|
||||
let mut r = StatusCode::NO_CONTENT.into_response();
|
||||
let h = r.headers_mut();
|
||||
h.insert("access-control-allow-origin", HeaderValue::from_static("*"));
|
||||
h.insert("access-control-allow-methods", HeaderValue::from_static("GET, POST, DELETE, OPTIONS"));
|
||||
h.insert("access-control-allow-headers", HeaderValue::from_static("content-type"));
|
||||
h.insert("access-control-allow-headers", HeaderValue::from_static("content-type, accept, mcp-protocol-version, mcp-session-id"));
|
||||
h.insert("access-control-max-age", HeaderValue::from_static("86400"));
|
||||
r
|
||||
}
|
||||
|
|
@ -188,6 +247,10 @@ pub fn router(app: Shared) -> Router {
|
|||
.route("/v1/jobs/{id}", get(job).delete(cancel))
|
||||
.route("/v1/run", post(run))
|
||||
.route("/v1/verify", post(verify))
|
||||
.route("/v1/examples", get(examples))
|
||||
.route("/v1/openapi.json", get(openapi))
|
||||
.route("/mcp", post(mcp).get(mcp_get))
|
||||
.route("/", get(index))
|
||||
.route_layer(axum::middleware::from_fn(|req: axum::extract::Request, next: axum::middleware::Next| async move {
|
||||
if req.method() == Method::OPTIONS { preflight().await } else { next.run(req).await }
|
||||
}))
|
||||
|
|
|
|||
|
|
@ -138,6 +138,9 @@ impl Job {
|
|||
"method": self.plan.method,
|
||||
"why": self.plan.why,
|
||||
"estimate": { "memory_bytes": self.plan.memory_bytes, "work": self.plan.work },
|
||||
// The request as submitted, so whoever holds the id can check the
|
||||
// result's receipt against it.
|
||||
"request": self.request,
|
||||
"submitted_ms": unix_ms(self.submitted),
|
||||
});
|
||||
if let Some(t) = self.started {
|
||||
|
|
|
|||
|
|
@ -14,8 +14,10 @@
|
|||
//! ```
|
||||
|
||||
mod catalog;
|
||||
mod discovery;
|
||||
mod http;
|
||||
mod jobs;
|
||||
mod mcp;
|
||||
mod qasm;
|
||||
mod receipt;
|
||||
mod run;
|
||||
|
|
@ -99,7 +101,10 @@ fn serve() {
|
|||
let signer = Arc::new(signer);
|
||||
service.seal_with(Arc::clone(&signer));
|
||||
tokio::spawn(Arc::clone(&service).schedule());
|
||||
let app = Arc::new(http::App { service, signer, trusted });
|
||||
// Where the service is seen from outside, and its page for people.
|
||||
let base = std::env::var("WAI_QSVC_PUBLIC_BASE").unwrap_or_else(|_| "https://wai.transaction.science/quantum-jobs".into());
|
||||
let console = std::env::var("WAI_QSVC_CONSOLE").unwrap_or_else(|_| "https://wai.transaction.science/quantum-run/".into());
|
||||
let app = Arc::new(http::App { service, signer, trusted, base, console });
|
||||
let listener = tokio::net::TcpListener::bind(addr).await.unwrap_or_else(|e| {
|
||||
eprintln!("binding {addr}: {e}");
|
||||
std::process::exit(1)
|
||||
|
|
|
|||
294
wai/wai-quantum-service/src/mcp.rs
Normal file
294
wai/wai-quantum-service/src/mcp.rs
Normal file
|
|
@ -0,0 +1,294 @@
|
|||
//! The service as Model Context Protocol tools: `POST /mcp`, streamable HTTP
|
||||
//! answered with plain JSON, no session and no key. An agent adds the URL and
|
||||
//! gets five tools: run a circuit, read a job, verify a receipt, describe the
|
||||
//! service, and list worked examples.
|
||||
//!
|
||||
//! Each tool result carries a short text summary written for a model to read
|
||||
//! (method, why, accuracy, the leading counts, the expectation values) and
|
||||
//! the whole job as structured content.
|
||||
|
||||
use crate::catalog;
|
||||
use crate::discovery;
|
||||
use crate::http::App;
|
||||
use crate::receipt;
|
||||
use serde_json::{Value, json};
|
||||
use std::time::Duration;
|
||||
|
||||
/// Protocol versions this server speaks, newest first.
|
||||
pub const VERSIONS: [&str; 3] = ["2025-11-25", "2025-06-18", "2025-03-26"];
|
||||
|
||||
/// How long a `run_quantum_circuit` call waits before handing back the job id.
|
||||
pub const CALL_WAIT: Duration = Duration::from_secs(50);
|
||||
|
||||
const INSTRUCTIONS: &str = "Run quantum programs written in OpenQASM 3 (or 2). Call run_quantum_circuit with the program and what you want back: shots for measurement counts, observables (Pauli strings such as \"Z0 Z1\") for expectation values, amplitudes (bit strings, qubit n-1 first) for amplitudes. The service picks the strongest method the circuit allows and says which and why; every result names its accuracy and carries a signature verify_quantum_receipt can check. Counts are keyed by the classical bits, the highest bit first. If a call answers with status \"running\", call get_quantum_job with its id. quantum_circuit_examples gives ready-to-run programs.";
|
||||
|
||||
fn tools() -> Value {
|
||||
let job_props = json!({
|
||||
"qasm": { "type": "string", "description": "The program, OpenQASM 3 or 2. Supported: qubit/bit registers and broadcast over them; gates h x y z s sdg t tdg sx sxdg rx ry rz p u cx cy cz ch swap iswap cp crx cry crz cu ccx ccz cswap rxx ryy rzz; the control modifiers `ctrl @`, `ctrl(n) @`, `negctrl @` on x y z h s sdg t tdg p and the controlled names; user `gate` definitions; angle expressions with pi; `reset` before a qubit's first gate; terminal measurement (`c = measure q;`). Not supported: mid-circuit measurement or reset, classical control, `inv @`, `pow @`." },
|
||||
"shots": { "type": "integer", "minimum": 0, "maximum": 1000000, "description": "Measurement shots. Default 1024 when no observables or amplitudes are asked for." },
|
||||
"observables": { "type": "array", "items": { "type": "string" }, "description": "Pauli strings to take expectation values of, e.g. [\"Z0\", \"Z0 Z1\", \"X2 Y3\"]. Qubits count from 0." },
|
||||
"amplitudes": { "type": "array", "items": { "type": "string" }, "description": "Basis states to give amplitudes of: one character per qubit, qubit n-1 first, e.g. \"0011\"." },
|
||||
"noise": { "type": "object", "description": "Optional noise, each a probability: depolarizing_1q, depolarizing_2q, amplitude_damping, dephasing (after each gate on its qubits), readout (each measured bit).", "properties": {
|
||||
"depolarizing_1q": { "type": "number" }, "depolarizing_2q": { "type": "number" }, "amplitude_damping": { "type": "number" }, "dephasing": { "type": "number" }, "readout": { "type": "number" }
|
||||
} },
|
||||
"method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "density_matrix", "trajectories"], "description": "Leave as auto unless you need a particular method." },
|
||||
"seed": { "type": "integer", "description": "Sampling seed (default 1); the same job and seed give the same counts." },
|
||||
"options": { "type": "object", "description": "Advanced, per method: max_bond (mps), threshold and max_terms (pauli), fusion (statevector), trajectories." }
|
||||
});
|
||||
json!([
|
||||
{
|
||||
"name": "run_quantum_circuit",
|
||||
"title": "Run a quantum circuit",
|
||||
"description": "Run an OpenQASM program and return measurement counts, Pauli expectation values and/or amplitudes. The service chooses the method (bit-exact or double-precision state vector up to 32 qubits, single precision on a GPU to 33, stabilizer tableau for Clifford circuits up to 32768 qubits, matrix product states up to 4096, sparse Pauli dynamics for expectation values up to 1024 with an error bound, an exact density matrix or trajectories for noise), reports why and how accurate, and signs the result. Waits up to 50 s; a longer job answers with status running and an id for get_quantum_job.",
|
||||
"inputSchema": { "type": "object", "required": ["qasm"], "properties": job_props },
|
||||
"annotations": { "readOnlyHint": false, "idempotentHint": true, "openWorldHint": false }
|
||||
},
|
||||
{
|
||||
"name": "get_quantum_job",
|
||||
"title": "Read a quantum job",
|
||||
"description": "A job's status, and its result when it is done. Results are kept for a day.",
|
||||
"inputSchema": { "type": "object", "required": ["id"], "properties": { "id": { "type": "string" } } },
|
||||
"annotations": { "readOnlyHint": true }
|
||||
},
|
||||
{
|
||||
"name": "verify_quantum_receipt",
|
||||
"title": "Verify a result's receipt",
|
||||
"description": "Whether a result's receipt (an Ed25519 signature by the service) seals exactly that request and that result. Pass the request as submitted and the result as returned.",
|
||||
"inputSchema": { "type": "object", "required": ["request", "result"], "properties": { "request": { "type": "object" }, "result": { "type": "object" } } },
|
||||
"annotations": { "readOnlyHint": true }
|
||||
},
|
||||
{
|
||||
"name": "quantum_service_info",
|
||||
"title": "What the quantum service can do",
|
||||
"description": "The methods with their qubit limits and what each answers, the accuracy classes, the per-client limits, the current load, and the conventions (bit order, Pauli strings).",
|
||||
"inputSchema": { "type": "object", "properties": {} },
|
||||
"annotations": { "readOnlyHint": true }
|
||||
},
|
||||
{
|
||||
"name": "quantum_circuit_examples",
|
||||
"title": "Ready-to-run quantum examples",
|
||||
"description": "Worked examples, each a complete job (OpenQASM program and arguments) with what it shows: a Bell pair, a rotation, QAOA, noise, a 1000-qubit GHZ state, 100-qubit dynamics, 200 qubits as a matrix product state. Pass an example's job to run_quantum_circuit as it is, or adapt it.",
|
||||
"inputSchema": { "type": "object", "properties": {} },
|
||||
"annotations": { "readOnlyHint": true }
|
||||
}
|
||||
])
|
||||
}
|
||||
|
||||
/// A bit string short enough to read: as it is up to 48 bits; past that, as
|
||||
/// runs (`0×1000`, `0×3 1×997`), or its two ends if the runs are many.
|
||||
pub fn compact_bits(b: &str) -> String {
|
||||
if b.len() <= 48 {
|
||||
return b.to_string();
|
||||
}
|
||||
let mut runs: Vec<(char, usize)> = Vec::new();
|
||||
for c in b.chars() {
|
||||
match runs.last_mut() {
|
||||
Some((d, n)) if *d == c => *n += 1,
|
||||
_ => runs.push((c, 1)),
|
||||
}
|
||||
}
|
||||
if runs.len() <= 6 {
|
||||
return runs.iter().map(|(c, n)| format!("{c}×{n}")).collect::<Vec<_>>().join(" ");
|
||||
}
|
||||
format!("{}…{} ({} bits)", &b[..16], &b[b.len() - 16..], b.len())
|
||||
}
|
||||
|
||||
/// A job as a model should read it: a few lines, the leading outcomes first.
|
||||
pub fn summary(job: &Value) -> String {
|
||||
let mut s = Vec::new();
|
||||
let status = job["status"].as_str().unwrap_or("?");
|
||||
s.push(format!("status: {status} · job {}", job["id"].as_str().unwrap_or("?")));
|
||||
if let Some(m) = job["method"].as_str() {
|
||||
s.push(format!("method: {m} — {}", job["why"].as_str().unwrap_or("")));
|
||||
}
|
||||
if let Some(e) = job["error"].as_str() {
|
||||
s.push(format!("error: {e}"));
|
||||
}
|
||||
let r = &job["result"];
|
||||
if r.is_object() {
|
||||
let acc = &r["accuracy"];
|
||||
let mut a = format!("accuracy: {}", acc["class"].as_str().unwrap_or("?"));
|
||||
for k in ["fidelity_estimate", "discarded_weight", "threshold", "trajectories", "histories", "device"] {
|
||||
if let Some(v) = acc.get(k) {
|
||||
a += &format!(", {k} {v}");
|
||||
}
|
||||
}
|
||||
s.push(a);
|
||||
if let Some(w) = r["warning"].as_str() {
|
||||
s.push(format!("warning: {w}"));
|
||||
}
|
||||
if let Some(c) = r["counts"].as_object() {
|
||||
let total: u64 = c.values().filter_map(Value::as_u64).sum::<u64>() + r["omitted_shots"].as_u64().unwrap_or(0);
|
||||
let mut v: Vec<(&String, u64)> = c.iter().filter_map(|(k, n)| n.as_u64().map(|n| (k, n))).collect();
|
||||
v.sort_by(|x, y| y.1.cmp(&x.1).then(x.0.cmp(y.0)));
|
||||
let shown: Vec<String> = v.iter().take(16).map(|(k, n)| format!("{}: {n} ({:.1}%)", compact_bits(k), 100.0 * *n as f64 / total.max(1) as f64)).collect();
|
||||
let more = if v.len() > 16 { format!(" … and {} more outcomes", v.len() - 16) } else { String::new() };
|
||||
s.push(format!("counts ({total} shots, highest classical bit first): {}{more}", shown.join(", ")));
|
||||
}
|
||||
if let Some(es) = r["expectations"].as_array() {
|
||||
let parts: Vec<String> = es
|
||||
.iter()
|
||||
.map(|e| {
|
||||
let mut p = format!("⟨{}⟩ = {:.6}", e["observable"].as_str().unwrap_or("?"), e["value"].as_f64().unwrap_or(f64::NAN));
|
||||
if let Some(b) = e["bound"].as_f64() {
|
||||
p += &format!(" (error at most {b:.2e})");
|
||||
}
|
||||
if let Some(se) = e["stderr"].as_f64() {
|
||||
p += &format!(" ± {se:.2e}");
|
||||
}
|
||||
p
|
||||
})
|
||||
.collect();
|
||||
s.push(format!("expectations: {}", parts.join("; ")));
|
||||
}
|
||||
if let Some(am) = r["amplitudes"].as_array() {
|
||||
let parts: Vec<String> = am.iter().map(|a| format!("{}: {:.6}{:+.6}i", compact_bits(a["bits"].as_str().unwrap_or("?")), a["re"].as_f64().unwrap_or(f64::NAN), a["im"].as_f64().unwrap_or(f64::NAN))).collect();
|
||||
s.push(format!("amplitudes: {}", parts.join("; ")));
|
||||
}
|
||||
if let Some(ms) = r["elapsed_ms"].as_u64() {
|
||||
s.push(format!("ran in {ms} ms; signed receipt attached (verify_quantum_receipt checks it)"));
|
||||
}
|
||||
} else if status == "queued" || status == "running" {
|
||||
s.push("still running: call get_quantum_job with this id".into());
|
||||
}
|
||||
s.join("\n")
|
||||
}
|
||||
|
||||
fn ok(id: &Value, result: Value) -> Value {
|
||||
json!({ "jsonrpc": "2.0", "id": id, "result": result })
|
||||
}
|
||||
|
||||
fn fail(id: &Value, code: i64, message: &str) -> Value {
|
||||
json!({ "jsonrpc": "2.0", "id": id, "error": { "code": code, "message": message } })
|
||||
}
|
||||
|
||||
/// A tool's answer: text for the model, the whole object as structured content.
|
||||
fn tool_result(text: String, structured: Value, is_error: bool) -> Value {
|
||||
json!({ "content": [{ "type": "text", "text": text }], "structuredContent": structured, "isError": is_error })
|
||||
}
|
||||
|
||||
/// Handle one JSON-RPC message from `who`. `None` for a notification, which
|
||||
/// gets no answer.
|
||||
pub async fn handle(app: &App, who: &str, msg: &Value, base: &str) -> Option<Value> {
|
||||
let id = msg.get("id")?.clone();
|
||||
let Some(method) = msg["method"].as_str() else {
|
||||
return Some(fail(&id, -32600, "not a JSON-RPC request"));
|
||||
};
|
||||
let params = &msg["params"];
|
||||
Some(match method {
|
||||
"initialize" => {
|
||||
let asked = params["protocolVersion"].as_str().unwrap_or("");
|
||||
let version = VERSIONS.iter().find(|v| **v == asked).copied().unwrap_or(VERSIONS[0]);
|
||||
ok(
|
||||
&id,
|
||||
json!({
|
||||
"protocolVersion": version,
|
||||
"capabilities": { "tools": { "listChanged": false } },
|
||||
"serverInfo": { "name": "wai-quantum", "title": "WAI quantum service", "version": env!("CARGO_PKG_VERSION") },
|
||||
"instructions": INSTRUCTIONS,
|
||||
}),
|
||||
)
|
||||
}
|
||||
"ping" => ok(&id, json!({})),
|
||||
"tools/list" => ok(&id, json!({ "tools": tools() })),
|
||||
"tools/call" => {
|
||||
let args = params.get("arguments").cloned().unwrap_or(json!({}));
|
||||
ok(&id, call(app, who, params["name"].as_str().unwrap_or(""), args, base).await)
|
||||
}
|
||||
other => fail(&id, -32601, &format!("no method {other}; this server offers initialize, ping, tools/list and tools/call")),
|
||||
})
|
||||
}
|
||||
|
||||
async fn call(app: &App, who: &str, name: &str, args: Value, base: &str) -> Value {
|
||||
match name {
|
||||
"run_quantum_circuit" => {
|
||||
let id = match crate::http::admit(app, who, args) {
|
||||
Ok(id) => id,
|
||||
Err((_, why)) => return tool_result(format!("refused: {why}"), json!({ "error": why }), true),
|
||||
};
|
||||
app.service.wait(&id, CALL_WAIT).await;
|
||||
let job = app.service.get(&id).unwrap_or(json!({ "id": id }));
|
||||
let failed = job["status"] == "failed";
|
||||
tool_result(summary(&job), job, failed)
|
||||
}
|
||||
"get_quantum_job" => match args["id"].as_str().and_then(|i| app.service.get(i)) {
|
||||
Some(job) => tool_result(summary(&job), job, false),
|
||||
None => tool_result("no such job (results are kept for a day)".into(), json!({ "error": "no such job" }), true),
|
||||
},
|
||||
"verify_quantum_receipt" => {
|
||||
let (request, result) = (&args["request"], &args["result"]);
|
||||
let Some(r) = result.get("receipt") else {
|
||||
return tool_result("the result carries no receipt".into(), json!({ "valid": false }), true);
|
||||
};
|
||||
let mut bare = result.clone();
|
||||
if let Some(o) = bare.as_object_mut() {
|
||||
o.remove("receipt");
|
||||
}
|
||||
let valid = receipt::verify(request, &bare, r);
|
||||
tool_result(if valid { "valid: the receipt seals exactly this request and this result".into() } else { "not valid: the request, the result or the receipt was changed".into() }, json!({ "valid": valid, "key": r.get("key") }), false)
|
||||
}
|
||||
"quantum_service_info" => {
|
||||
let c = &app.service.config;
|
||||
let (queued, running, _) = app.service.load();
|
||||
let info = json!({
|
||||
"methods": catalog::describe(c.job_mib << 20),
|
||||
"accuracy_classes": catalog::accuracy_classes(),
|
||||
"conventions": {
|
||||
"counts": "keyed by the program's classical bits, the highest bit first; with no measurement, every qubit is measured into the bit of its own index",
|
||||
"amplitudes": "bit strings with one character per qubit, qubit n-1 first",
|
||||
"observables": "Pauli strings such as \"Z0 Z1\" or {\"pauli\": \"X3\", \"coeff\": 0.5}; qubits count from 0"
|
||||
},
|
||||
"limits": { "job_seconds": c.job_timeout.as_secs(), "active_jobs_per_client": c.per_client_active, "submissions_per_window": c.per_client_burst, "window_seconds": c.per_client_window.as_secs() },
|
||||
"load": { "queued": queued, "running": running },
|
||||
"api": discovery::index(base, &app.console),
|
||||
});
|
||||
let lines: Vec<String> = info["methods"]
|
||||
.as_array()
|
||||
.into_iter()
|
||||
.flatten()
|
||||
.map(|m| format!("{} — up to {} qubits, {}; {}", m["method"].as_str().unwrap_or("?"), m["max_qubits"], m["accuracy"].as_str().unwrap_or(""), m["circuits"].as_str().unwrap_or("")))
|
||||
.collect();
|
||||
tool_result(format!("methods (auto picks the first that fits):\n{}\ncounts: highest classical bit first; amplitudes: qubit n-1 first", lines.join("\n")), info, false)
|
||||
}
|
||||
"quantum_circuit_examples" => {
|
||||
let ex = discovery::examples();
|
||||
let lines: Vec<String> = ex.as_array().into_iter().flatten().map(|e| format!("{} — {}: {}", e["name"].as_str().unwrap_or(""), e["title"].as_str().unwrap_or(""), e["shows"].as_str().unwrap_or(""))).collect();
|
||||
tool_result(format!("{}\n(each example's `job` in the structured content runs as given)", lines.join("\n")), json!({ "examples": ex }), false)
|
||||
}
|
||||
other => tool_result(format!("no tool {other}"), json!({ "error": "no such tool" }), true),
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn the_summary_reads_a_finished_job() {
|
||||
let job = json!({
|
||||
"id": "abc", "status": "done", "method": "stabilizer", "why": "Clifford circuit",
|
||||
"result": { "accuracy": { "class": "exact" }, "counts": { "00": 503, "11": 497 }, "expectations": [{ "observable": "Z0 Z1", "value": 1.0 }], "elapsed_ms": 3 }
|
||||
});
|
||||
let s = summary(&job);
|
||||
assert!(s.contains("method: stabilizer"), "{s}");
|
||||
assert!(s.contains("00: 503 (50.3%)") && s.contains("11: 497"), "{s}");
|
||||
assert!(s.contains("⟨Z0 Z1⟩ = 1.000000"), "{s}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn wide_outcomes_are_written_short() {
|
||||
assert_eq!(compact_bits("0110"), "0110");
|
||||
assert_eq!(compact_bits(&"0".repeat(1000)), "0×1000");
|
||||
assert_eq!(compact_bits(&format!("{}{}", "1".repeat(3), "0".repeat(97))), "1×3 0×97");
|
||||
let alt: String = (0..100).map(|i| if i % 2 == 0 { '0' } else { '1' }).collect();
|
||||
assert_eq!(compact_bits(&alt), format!("{}…{} (100 bits)", &alt[..16], &alt[84..]));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn every_tool_has_a_schema_and_a_description() {
|
||||
for t in tools().as_array().unwrap() {
|
||||
assert!(t["description"].as_str().unwrap().len() > 40, "{t}");
|
||||
assert_eq!(t["inputSchema"]["type"], "object", "{t}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -11,8 +11,14 @@
|
|||
//! `c = measure q;`, `measure q -> c;`, `measure q;`). With no measurement,
|
||||
//! every qubit is measured into a bit of its own index.
|
||||
//!
|
||||
//! Refused, each with its line: classical control, gate modifiers, a gate on a
|
||||
//! qubit after its measurement, and anything else not listed.
|
||||
//! Control modifiers, `ctrl @`, `ctrl(n) @`, `negctrl @` and `negctrl(n) @`,
|
||||
//! on `x y z h s sdg t tdg p` and on the controlled names `cx cy cz ch cp ccx
|
||||
//! ccz`: a gate the standard names is read as that gate, and any other as
|
||||
//! standard gates, exactly (see [`Parser::phase_n`]).
|
||||
//!
|
||||
//! Refused, each with its line: classical control, the `inv` and `pow`
|
||||
//! modifiers, a gate on a qubit after its measurement, and anything else not
|
||||
//! listed.
|
||||
//!
|
||||
//! The standard gates: `id x y z h s sdg t tdg sx sxdg rx ry rz p phase u1 u2
|
||||
//! u u3 U cx CX cnot cy cz ch cp cphase cu1 crx cry crz cu cu3 swap iswap rxx
|
||||
|
|
@ -34,6 +40,9 @@ pub const MAX_QUBITS: u32 = 65_536;
|
|||
pub const MAX_OPS: usize = 2_000_000;
|
||||
/// How deep user gate definitions may nest.
|
||||
const MAX_DEPTH: usize = 32;
|
||||
/// The most controls modifiers may add: a phase on m qubits is written as
|
||||
/// 2^m − 1 `cx` and as many `p`.
|
||||
const MAX_CONTROLS: u32 = 16;
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
|
||||
pub enum Kind {
|
||||
|
|
@ -146,6 +155,9 @@ fn normalise(kind: Kind, given: usize, args: Vec<f64>) -> Vec<f64> {
|
|||
}
|
||||
}
|
||||
|
||||
/// One op's Clifford gates and the qubits the op acts on.
|
||||
pub type CliffordOp = (Vec<Clifford>, Vec<u32>);
|
||||
|
||||
/// One gate on named qubits.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct Op {
|
||||
|
|
@ -498,7 +510,10 @@ impl Parser {
|
|||
let regs = if bits { &self.cregs } else { &self.qregs };
|
||||
match regs.get(&name) {
|
||||
Some(r) => *r,
|
||||
None => return self.err(format!("no {} register named {name}", if bits { "bit" } else { "qubit" })),
|
||||
None => {
|
||||
let (kind, decl) = if bits { ("bit", "bit") } else { ("qubit", "qubit") };
|
||||
return self.err(format!("no {kind} register named {name}; declare it first, e.g. `{decl}[2] {name};`"));
|
||||
}
|
||||
}
|
||||
};
|
||||
if self.eat("[") {
|
||||
|
|
@ -592,8 +607,22 @@ impl Parser {
|
|||
self.expect(";")?;
|
||||
self.measure(&qs, &bs, line)?;
|
||||
}
|
||||
"if" | "while" | "for" | "def" | "defcal" | "cal" | "box" | "let" | "const" | "input" | "output" | "ctrl" | "negctrl" | "inv" | "pow" | "delay" | "gphase" => {
|
||||
return Err(format!("line {line}: `{word}` is not supported by this service"));
|
||||
"ctrl" | "negctrl" => {
|
||||
let polarity = self.modifiers(word == "ctrl")?;
|
||||
let name = self.ident()?;
|
||||
if matches!(name.as_str(), "inv" | "pow") {
|
||||
return Err(format!("line {line}: `{name}` is not supported by this service; write the inverse or powered gate out, or define it with `gate`"));
|
||||
}
|
||||
self.call_with(name, line, &polarity)?;
|
||||
}
|
||||
"if" | "while" | "for" | "def" | "defcal" | "cal" | "box" | "let" | "const" | "input" | "output" | "inv" | "pow" | "delay" | "gphase" => {
|
||||
let hint = match word.as_str() {
|
||||
"if" | "while" => "; classical control is not simulated here: measure every qubit at the end instead",
|
||||
"for" => "; unroll the loop into gates",
|
||||
"inv" | "pow" => "; write the inverse or powered gate out, or define it with `gate`",
|
||||
_ => "",
|
||||
};
|
||||
return Err(format!("line {line}: `{word}` is not supported by this service{hint}"));
|
||||
}
|
||||
_ => {
|
||||
// `c = measure q;` / `c[0] = measure q[0];`, or a gate call.
|
||||
|
|
@ -690,6 +719,41 @@ impl Parser {
|
|||
}
|
||||
|
||||
fn call(&mut self, name: String, line: usize) -> Result<(), String> {
|
||||
self.call_with(name, line, &[])
|
||||
}
|
||||
|
||||
/// The controls a run of `ctrl` and `negctrl` modifiers adds, in order
|
||||
/// (`true`: a control on |1⟩), the first modifier's word already read.
|
||||
fn modifiers(&mut self, mut on: bool) -> Result<Vec<bool>, String> {
|
||||
let mut polarity = Vec::new();
|
||||
loop {
|
||||
let count = if self.eat("(") {
|
||||
let n = self.integer()?;
|
||||
self.expect(")")?;
|
||||
n
|
||||
} else {
|
||||
1
|
||||
};
|
||||
if count == 0 {
|
||||
return self.err("a control modifier adds at least one control");
|
||||
}
|
||||
if polarity.len() as u32 + count > MAX_CONTROLS {
|
||||
return self.err(format!("modifiers add at most {MAX_CONTROLS} controls here"));
|
||||
}
|
||||
self.expect("@")?;
|
||||
polarity.extend(std::iter::repeat_n(on, count as usize));
|
||||
match self.peek() {
|
||||
Some(Tok::Ident(w)) if w == "ctrl" || w == "negctrl" => {
|
||||
on = w == "ctrl";
|
||||
self.pos += 1;
|
||||
}
|
||||
_ => return Ok(polarity),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A gate call, under the controls `polarity` adds (none for a plain call).
|
||||
fn call_with(&mut self, name: String, line: usize, polarity: &[bool]) -> Result<(), String> {
|
||||
let mut params = Vec::new();
|
||||
if self.eat("(") && !self.eat(")") {
|
||||
loop {
|
||||
|
|
@ -716,7 +780,142 @@ impl Parser {
|
|||
}
|
||||
for k in 0..width {
|
||||
let qs: Vec<u32> = operands.iter().map(|o| if o.len() == 1 { o[0] } else { o[k] }).collect();
|
||||
self.apply(&name, ¶ms, &qs, line, 0)?;
|
||||
if polarity.is_empty() {
|
||||
self.apply(&name, ¶ms, &qs, line, 0)?;
|
||||
} else {
|
||||
self.controlled(&name, ¶ms, &qs, polarity, line)?;
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// A gate under control modifiers, its added controls the first qubits
|
||||
/// named; a control on |0⟩ is a control between two `x`. The result is the
|
||||
/// gate the standard names where there is one, and otherwise the base
|
||||
/// gate's target turned into a phase: `x` = `h·z·h`, `y` = `s·x·sdg`, `h` =
|
||||
/// `ry(π/4)·z·ry(−π/4)`, each exact.
|
||||
fn controlled(&mut self, name: &str, params: &[f64], qubits: &[u32], polarity: &[bool], line: usize) -> Result<(), String> {
|
||||
use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, PI};
|
||||
enum Base {
|
||||
X,
|
||||
Y,
|
||||
H,
|
||||
Phase(f64),
|
||||
}
|
||||
// (the base gate, the controls its name already has)
|
||||
let (base, inner) = match (name, params) {
|
||||
("x", []) => (Base::X, 0),
|
||||
("y", []) => (Base::Y, 0),
|
||||
("h", []) => (Base::H, 0),
|
||||
("z", []) => (Base::Phase(PI), 0),
|
||||
("s", []) => (Base::Phase(FRAC_PI_2), 0),
|
||||
("sdg", []) => (Base::Phase(-FRAC_PI_2), 0),
|
||||
("t", []) => (Base::Phase(FRAC_PI_4), 0),
|
||||
("tdg", []) => (Base::Phase(-FRAC_PI_4), 0),
|
||||
("p" | "phase" | "u1", [t]) => (Base::Phase(*t), 0),
|
||||
("cx" | "CX" | "cnot", []) => (Base::X, 1),
|
||||
("cy", []) => (Base::Y, 1),
|
||||
("ch", []) => (Base::H, 1),
|
||||
("cz", []) => (Base::Phase(PI), 1),
|
||||
("cp" | "cphase" | "cu1", [t]) => (Base::Phase(*t), 1),
|
||||
("ccx" | "toffoli", []) => (Base::X, 2),
|
||||
("ccz", []) => (Base::Phase(PI), 2),
|
||||
_ => {
|
||||
return Err(format!(
|
||||
"line {line}: a control modifier on `{name}` is not supported (under one: x y z h s sdg t tdg p, and cx cy cz ch cp ccx ccz); define the controlled gate with `gate`"
|
||||
));
|
||||
}
|
||||
};
|
||||
if let Base::Phase(t) = base
|
||||
&& !t.is_finite()
|
||||
{
|
||||
return Err(format!("line {line}: {name} has a parameter that is not finite"));
|
||||
}
|
||||
let want = polarity.len() + inner + 1;
|
||||
if qubits.len() != want {
|
||||
return Err(format!("line {line}: {name} under {} added control(s) takes {want} qubits, given {}", polarity.len(), qubits.len()));
|
||||
}
|
||||
for (i, q) in qubits.iter().enumerate() {
|
||||
if qubits[..i].contains(q) {
|
||||
return Err(format!("line {line}: qubit {q} named twice in one gate"));
|
||||
}
|
||||
}
|
||||
let flip: Vec<u32> = qubits.iter().zip(polarity).filter(|(_, on)| !**on).map(|(q, _)| *q).collect();
|
||||
for &q in &flip {
|
||||
self.apply("x", &[], &[q], line, 0)?;
|
||||
}
|
||||
let t = qubits[want - 1];
|
||||
match base {
|
||||
Base::X => self.x_n(qubits, line)?,
|
||||
Base::Y if qubits.len() <= 2 => self.apply(if qubits.len() == 1 { "y" } else { "cy" }, &[], qubits, line, 0)?,
|
||||
Base::Y => {
|
||||
self.apply("sdg", &[], &[t], line, 0)?;
|
||||
self.x_n(qubits, line)?;
|
||||
self.apply("s", &[], &[t], line, 0)?;
|
||||
}
|
||||
Base::H => match qubits.len() {
|
||||
1 => self.apply("h", &[], qubits, line, 0)?,
|
||||
2 => self.apply("ch", &[], qubits, line, 0)?,
|
||||
_ => {
|
||||
self.apply("ry", &[-FRAC_PI_4], &[t], line, 0)?;
|
||||
self.phase_n(PI, qubits, line)?;
|
||||
self.apply("ry", &[FRAC_PI_4], &[t], line, 0)?;
|
||||
}
|
||||
},
|
||||
Base::Phase(theta) => self.phase_n(theta, qubits, line)?,
|
||||
}
|
||||
for &q in &flip {
|
||||
self.apply("x", &[], &[q], line, 0)?;
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// `x` on the last of `qs`, controlled by the others.
|
||||
fn x_n(&mut self, qs: &[u32], line: usize) -> Result<(), String> {
|
||||
let t = qs[qs.len() - 1];
|
||||
match qs.len() {
|
||||
1 => self.apply("x", &[], qs, line, 0),
|
||||
2 => self.apply("cx", &[], qs, line, 0),
|
||||
3 => self.apply("ccx", &[], qs, line, 0),
|
||||
_ => {
|
||||
self.apply("h", &[], &[t], line, 0)?;
|
||||
self.phase_n(std::f64::consts::PI, qs, line)?;
|
||||
self.apply("h", &[], &[t], line, 0)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The phase e^{iθ} on the state in which all of `qs` are 1, and on no
|
||||
/// other: `p`, `cp`, `cz` or `ccz` where one is that gate, and otherwise
|
||||
/// from x₁⋯x_m = 2^{1−m} Σ_S (−1)^{|S|+1} ⊕_{i∈S} x_i over the nonempty
|
||||
/// subsets S of the m qubits. Each subset's parity is gathered on its
|
||||
/// highest qubit by `cx` in Gray-code order and turned by `p(±θ/2^{m−1})`:
|
||||
/// 2^m − 1 of each, exact, with no global phase, and dyadic when θ is.
|
||||
fn phase_n(&mut self, theta: f64, qs: &[u32], line: usize) -> Result<(), String> {
|
||||
let z = theta == std::f64::consts::PI;
|
||||
match qs.len() {
|
||||
1 if z => return self.apply("z", &[], qs, line, 0),
|
||||
1 => return self.apply("p", &[theta], qs, line, 0),
|
||||
2 if z => return self.apply("cz", &[], qs, line, 0),
|
||||
2 => return self.apply("cp", &[theta], qs, line, 0),
|
||||
3 if z => return self.apply("ccz", &[], qs, line, 0),
|
||||
_ => {}
|
||||
}
|
||||
let a = theta / (1u64 << (qs.len() - 1)) as f64;
|
||||
for j in 0..qs.len() {
|
||||
let lead = qs[j];
|
||||
self.apply("p", &[a], &[lead], line, 0)?;
|
||||
for i in 1u64..1 << j {
|
||||
// The Gray code's step from i − 1 to i flips this bit.
|
||||
let bit = i.trailing_zeros() as usize;
|
||||
self.apply("cx", &[], &[qs[bit], lead], line, 0)?;
|
||||
let size = (i ^ (i >> 1)).count_ones() + 1;
|
||||
self.apply("p", &[if size % 2 == 1 { a } else { -a }], &[lead], line, 0)?;
|
||||
}
|
||||
// The code ends on the single bit j − 1: release it.
|
||||
if j > 0 {
|
||||
self.apply("cx", &[], &[qs[j - 1], lead], line, 0)?;
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
|
@ -1357,8 +1556,31 @@ impl Program {
|
|||
/// phase; refused when a gate is not Clifford at its angles (rotations
|
||||
/// must be whole quarter turns, controlled phases whole half turns).
|
||||
pub fn clifford(&self) -> Result<Vec<Clifford>, String> {
|
||||
self.clifford_bounds().map(|(gates, _)| gates)
|
||||
}
|
||||
|
||||
/// The Clifford gates of each op, with the qubits the op acts on: what
|
||||
/// noise after each gate attaches to, as the other noisy methods apply it.
|
||||
pub fn clifford_ops(&self) -> Result<Vec<CliffordOp>, String> {
|
||||
let (gates, ends) = self.clifford_bounds()?;
|
||||
let mut start = 0;
|
||||
Ok(self
|
||||
.ops
|
||||
.iter()
|
||||
.zip(ends)
|
||||
.map(|(op, end)| {
|
||||
let g = gates[start..end].to_vec();
|
||||
start = end;
|
||||
(g, op.qubits.clone())
|
||||
})
|
||||
.collect())
|
||||
}
|
||||
|
||||
/// The Clifford lowering, and where each op's gates end in it.
|
||||
fn clifford_bounds(&self) -> Result<(Vec<Clifford>, Vec<usize>), String> {
|
||||
use Kind::*;
|
||||
let mut out = Vec::new();
|
||||
let mut ends = Vec::with_capacity(self.ops.len());
|
||||
for op in &self.ops {
|
||||
let q = &op.qubits;
|
||||
let p = &op.params;
|
||||
|
|
@ -1445,8 +1667,9 @@ impl Program {
|
|||
out.extend(post);
|
||||
}
|
||||
}
|
||||
ends.push(out.len());
|
||||
}
|
||||
Ok(out)
|
||||
Ok((out, ends))
|
||||
}
|
||||
|
||||
/// The classical outcome a measurement of the qubits gives, as the bit
|
||||
|
|
@ -1503,6 +1726,13 @@ mod tests {
|
|||
for (src, needle) in [
|
||||
("qubit[1] q;\nbit c;\nc = measure q;\nx q;", "line 4"),
|
||||
("qubit[1] q;\nif (c) x q;", "not supported"),
|
||||
("qubit[1] q;\nif (c) x q;", "measure every qubit at the end"),
|
||||
("qubit[1] q;\nfor uint i in [0:3] { x q; }", "unroll the loop"),
|
||||
("qubit[1] q;\ninv @ x q[0];", "write the inverse"),
|
||||
("qubit[2] q;\nctrl @ rz(0.1) q[0], q[1];", "a control modifier on `rz` is not supported"),
|
||||
("qubit[3] q;\nctrl @ x q[0], q[1], q[2];", "takes 2 qubits"),
|
||||
("qubit[2] q;\nctrl(0) @ x q[0], q[1];", "at least one control"),
|
||||
("qubit[18] q;\nctrl(17) @ x q;", "at most 16 controls"),
|
||||
("qubit[1] q;\nfoo q;", "unknown gate foo"),
|
||||
("qubit[1] q;\nrx q;", "rx takes 1 parameters"),
|
||||
("qubit[2] q;\ncx q[0], q[0];", "named twice"),
|
||||
|
|
@ -1562,6 +1792,72 @@ mod tests {
|
|||
sv
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_controlled_gate_is_its_gate_where_the_controls_hold_and_nothing_elsewhere() {
|
||||
use std::f64::consts::FRAC_1_SQRT_2 as R;
|
||||
// A state with no special structure on 6 qubits; then each modified
|
||||
// gate, against the 2×2 gate applied to the amplitudes directly where
|
||||
// the controls hold. Exactly, phase and all: a controlled gate exposes
|
||||
// any global phase its expansion has.
|
||||
let prep = "qubit[6] q;\nry(0.3) q[0]; rx(1.1) q[1]; ry(2.2) q[2]; rz(0.7) q[2]; rx(0.4) q[3]; ry(1.7) q[4]; ry(0.6) q[5]; cx q[0], q[3]; rz(0.9) q[4]; crx(0.8) q[5], q[1];\n";
|
||||
let (c4, s4) = (std::f64::consts::FRAC_PI_4.cos(), std::f64::consts::FRAC_PI_4.sin());
|
||||
let (cp, sp) = (0.37f64.cos(), 0.37f64.sin());
|
||||
let gates: [(&str, [(f64, f64); 4]); 9] = [
|
||||
("x", [(0.0, 0.0), (1.0, 0.0), (1.0, 0.0), (0.0, 0.0)]),
|
||||
("y", [(0.0, 0.0), (0.0, -1.0), (0.0, 1.0), (0.0, 0.0)]),
|
||||
("z", [(1.0, 0.0), (0.0, 0.0), (0.0, 0.0), (-1.0, 0.0)]),
|
||||
("h", [(R, 0.0), (R, 0.0), (R, 0.0), (-R, 0.0)]),
|
||||
("s", [(1.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 1.0)]),
|
||||
("sdg", [(1.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, -1.0)]),
|
||||
("t", [(1.0, 0.0), (0.0, 0.0), (0.0, 0.0), (c4, s4)]),
|
||||
("tdg", [(1.0, 0.0), (0.0, 0.0), (0.0, 0.0), (c4, -s4)]),
|
||||
("p(0.37)", [(1.0, 0.0), (0.0, 0.0), (0.0, 0.0), (cp, sp)]),
|
||||
];
|
||||
let base = dense(&parse(prep).unwrap());
|
||||
let n = 6u32;
|
||||
let bit = |q: u32| 1usize << (n - 1 - q);
|
||||
let order = [3u32, 0, 5, 4, 1];
|
||||
for (g, u) in gates {
|
||||
for polarity in [vec![true], vec![false], vec![true, true], vec![false, true], vec![true, true, true], vec![true, false, true, false], vec![true; 5]] {
|
||||
let mods: String = polarity.iter().map(|&on| if on { "ctrl @ " } else { "negctrl @ " }).collect();
|
||||
let controls = &order[..polarity.len()];
|
||||
let target = 2u32;
|
||||
let operands: Vec<String> = controls.iter().chain([&target]).map(|q| format!("q[{q}]")).collect();
|
||||
let got = dense(&parse(&format!("{prep}{mods}{g} {};", operands.join(", "))).unwrap());
|
||||
let (mut re, mut im) = (base.re.clone(), base.im.clone());
|
||||
for i0 in 0..1usize << n {
|
||||
if i0 & bit(target) != 0 || !controls.iter().zip(&polarity).all(|(&c, &on)| (i0 & bit(c) != 0) == on) {
|
||||
continue;
|
||||
}
|
||||
let i1 = i0 | bit(target);
|
||||
let (a0, a1) = ((re[i0], im[i0]), (re[i1], im[i1]));
|
||||
let mul = |m: (f64, f64), a: (f64, f64)| (m.0 * a.0 - m.1 * a.1, m.0 * a.1 + m.1 * a.0);
|
||||
let (x0, x1) = (mul(u[0], a0), mul(u[1], a1));
|
||||
let (y0, y1) = (mul(u[2], a0), mul(u[3], a1));
|
||||
(re[i0], im[i0]) = (x0.0 + x1.0, x0.1 + x1.1);
|
||||
(re[i1], im[i1]) = (y0.0 + y1.0, y0.1 + y1.1);
|
||||
}
|
||||
for i in 0..1usize << n {
|
||||
assert!((got.re[i] - re[i]).abs() < 1e-12 && (got.im[i] - im[i]).abs() < 1e-12, "{mods}{g}: amplitude {i}");
|
||||
}
|
||||
}
|
||||
}
|
||||
// A controlled name under a modifier is the gate with both controls.
|
||||
for (a, b) in [("ctrl @ cx q[3], q[0], q[2];", "ctrl(2) @ x q[3], q[0], q[2];"), ("negctrl @ cp(0.37) q[4], q[1], q[2];", "negctrl @ ctrl @ p(0.37) q[4], q[1], q[2];"), ("ctrl(2) @ ccz q[0], q[1], q[3], q[4], q[2];", "ctrl(4) @ z q[0], q[1], q[3], q[4], q[2];")] {
|
||||
let (x, y) = (dense(&parse(&format!("{prep}{a}")).unwrap()), dense(&parse(&format!("{prep}{b}")).unwrap()));
|
||||
assert!(x.re.iter().zip(&y.re).chain(x.im.iter().zip(&y.im)).all(|(p, q)| (p - q).abs() < 1e-12), "{a} / {b}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_controlled_gate_the_standard_names_is_read_as_it() {
|
||||
let p = parse("qubit[3] q;\nctrl @ x q[0], q[1];\nctrl(2) @ z q[0], q[1], q[2];\nnegctrl @ z q[0], q[1];\nctrl @ ctrl @ x q[2], q[1], q[0];\nctrl @ h q[1], q[0];\nctrl @ y q[1], q[2];").unwrap();
|
||||
let kinds: Vec<Kind> = p.ops.iter().map(|o| o.kind).collect();
|
||||
assert_eq!(kinds, vec![Kind::Cx, Kind::Ccz, Kind::X, Kind::Cz, Kind::X, Kind::Ccx, Kind::Ch, Kind::Cy]);
|
||||
// So a Clifford program under modifiers keeps the stabilizer method.
|
||||
assert!(parse("qubit[2] q;\nctrl @ x q[0], q[1];\nnegctrl @ z q[0], q[1];\nctrl @ y q[1], q[0];").unwrap().clifford().is_ok());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_matrices_are_unitary() {
|
||||
for op in every_kind(4, 3).ops {
|
||||
|
|
|
|||
|
|
@ -1,12 +1,24 @@
|
|||
//! Signed receipts: every finished job's result is sealed with the service's
|
||||
//! Ed25519 key, binding the request as submitted and the result as returned,
|
||||
//! so a client can show anyone what was run and what came back.
|
||||
//!
|
||||
//! What is signed is the JSON Canonicalization Scheme (RFC 8785) text of each:
|
||||
//! keys sorted, no whitespace, and every number written as ECMAScript writes
|
||||
//! the double it reads. A result therefore verifies after any JSON reader and
|
||||
//! writer has carried it: a browser's `JSON.parse` and `JSON.stringify` turn
|
||||
//! `1.0` into `1`, and other readers reorder keys.
|
||||
|
||||
use serde_json::{Value, json};
|
||||
use std::path::{Path, PathBuf};
|
||||
use wai_quantum::ed25519_dalek::{Signer as _, SigningKey, Verifier as _, VerifyingKey};
|
||||
|
||||
const DOMAIN: &[u8] = b"wai:quantum-service-receipt\x01";
|
||||
/// The first receipts signed serde_json's own text, which only a reader that
|
||||
/// keeps each number's spelling and each object's key order can reproduce.
|
||||
/// They still verify; new ones are canonical.
|
||||
const DOMAIN_V1: &[u8] = b"wai:quantum-service-receipt\x01";
|
||||
const DOMAIN: &[u8] = b"wai:quantum-service-receipt\x02";
|
||||
/// The receipt's name for its canonical form.
|
||||
const CANON: &str = "rfc8785";
|
||||
|
||||
/// The service's signing key.
|
||||
pub struct Signer {
|
||||
|
|
@ -29,15 +41,229 @@ fn put(o: &mut Vec<u8>, bytes: &[u8]) {
|
|||
}
|
||||
|
||||
/// The digest a receipt signs: BLAKE3(domain ‖ put(request) ‖ put(result)),
|
||||
/// each the compact JSON text serde_json writes (keys in insertion order, as
|
||||
/// the service builds them and returns them).
|
||||
/// each in its RFC 8785 canonical text.
|
||||
pub fn digest(request: &Value, result: &Value) -> [u8; 32] {
|
||||
let mut o = DOMAIN.to_vec();
|
||||
put(&mut o, &canonical(request));
|
||||
put(&mut o, &canonical(result));
|
||||
*blake3::hash(&o).as_bytes()
|
||||
}
|
||||
|
||||
/// The digest the first receipts signed (see [`DOMAIN_V1`]).
|
||||
fn digest_v1(request: &Value, result: &Value) -> [u8; 32] {
|
||||
let mut o = DOMAIN_V1.to_vec();
|
||||
put(&mut o, request.to_string().as_bytes());
|
||||
put(&mut o, result.to_string().as_bytes());
|
||||
*blake3::hash(&o).as_bytes()
|
||||
}
|
||||
|
||||
/// The RFC 8785 canonical text of `v`.
|
||||
pub fn canonical(v: &Value) -> Vec<u8> {
|
||||
let mut o = Vec::new();
|
||||
write(v, &mut o);
|
||||
o
|
||||
}
|
||||
|
||||
fn write(v: &Value, o: &mut Vec<u8>) {
|
||||
match v {
|
||||
Value::Null => o.extend_from_slice(b"null"),
|
||||
Value::Bool(b) => o.extend_from_slice(if *b { b"true" } else { b"false" }),
|
||||
// Every number as the double a JSON reader in ECMAScript holds.
|
||||
Value::Number(n) => o.extend_from_slice(es_number(n.as_f64().unwrap_or(f64::NAN)).as_bytes()),
|
||||
Value::String(s) => write_str(s, o),
|
||||
Value::Array(a) => {
|
||||
o.push(b'[');
|
||||
for (i, x) in a.iter().enumerate() {
|
||||
if i > 0 {
|
||||
o.push(b',');
|
||||
}
|
||||
write(x, o);
|
||||
}
|
||||
o.push(b']');
|
||||
}
|
||||
Value::Object(m) => {
|
||||
// Keys in the order of their UTF-16 code units.
|
||||
let mut keys: Vec<&String> = m.keys().collect();
|
||||
keys.sort_by(|a, b| a.encode_utf16().cmp(b.encode_utf16()));
|
||||
o.push(b'{');
|
||||
for (i, k) in keys.into_iter().enumerate() {
|
||||
if i > 0 {
|
||||
o.push(b',');
|
||||
}
|
||||
write_str(k, o);
|
||||
o.push(b':');
|
||||
write(&m[k], o);
|
||||
}
|
||||
o.push(b'}');
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A string as RFC 8785 writes it: `"` and `\\` escaped, the five short
|
||||
/// control escapes, every other control character as `\u00xx`, and the rest
|
||||
/// as it is.
|
||||
fn write_str(s: &str, o: &mut Vec<u8>) {
|
||||
o.push(b'"');
|
||||
for c in s.chars() {
|
||||
match c {
|
||||
'"' => o.extend_from_slice(b"\\\""),
|
||||
'\\' => o.extend_from_slice(b"\\\\"),
|
||||
'\u{8}' => o.extend_from_slice(b"\\b"),
|
||||
'\t' => o.extend_from_slice(b"\\t"),
|
||||
'\n' => o.extend_from_slice(b"\\n"),
|
||||
'\u{c}' => o.extend_from_slice(b"\\f"),
|
||||
'\r' => o.extend_from_slice(b"\\r"),
|
||||
c if (c as u32) < 0x20 => o.extend_from_slice(format!("\\u{:04x}", c as u32).as_bytes()),
|
||||
c => {
|
||||
let mut b = [0u8; 4];
|
||||
o.extend_from_slice(c.encode_utf8(&mut b).as_bytes());
|
||||
}
|
||||
}
|
||||
}
|
||||
o.push(b'"');
|
||||
}
|
||||
|
||||
/// A double as ECMAScript's Number::toString writes it (ECMA-262, 6.1.6.1.20):
|
||||
/// the shortest digits that read back as `x`, laid out by where the decimal
|
||||
/// point falls. JSON holds no infinity or NaN; they write as `null`.
|
||||
pub fn es_number(x: f64) -> String {
|
||||
if !x.is_finite() {
|
||||
return "null".into();
|
||||
}
|
||||
if x == 0.0 {
|
||||
return "0".into();
|
||||
}
|
||||
let (s, n) = shortest(x.abs());
|
||||
let digits: String = s.iter().map(|d| char::from(b'0' + d)).collect();
|
||||
let k = s.len() as i32;
|
||||
let mut out = String::new();
|
||||
if x < 0.0 {
|
||||
out.push('-');
|
||||
}
|
||||
if k <= n && n <= 21 {
|
||||
out.push_str(&digits);
|
||||
out.extend(std::iter::repeat_n('0', (n - k) as usize));
|
||||
} else if 0 < n && n <= 21 {
|
||||
out.push_str(&digits[..n as usize]);
|
||||
out.push('.');
|
||||
out.push_str(&digits[n as usize..]);
|
||||
} else if -6 < n && n <= 0 {
|
||||
out.push_str("0.");
|
||||
out.extend(std::iter::repeat_n('0', (-n) as usize));
|
||||
out.push_str(&digits);
|
||||
} else {
|
||||
out.push_str(&digits[..1]);
|
||||
if k > 1 {
|
||||
out.push('.');
|
||||
out.push_str(&digits[1..]);
|
||||
}
|
||||
out.push('e');
|
||||
out.push(if n - 1 < 0 { '-' } else { '+' });
|
||||
out.push_str(&(n - 1).abs().to_string());
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// The digits s and the exponent n, x = 0.s × 10^n, that ECMAScript chooses:
|
||||
/// the fewest that read back as `x` (> 0), and of those the nearest to `x`,
|
||||
/// a tie going to the even.
|
||||
fn shortest(x: f64) -> (Vec<u8>, i32) {
|
||||
let e = format!("{x:e}");
|
||||
let (mant, exp) = e.split_once('e').unwrap_or((&e, "0"));
|
||||
let mut s: Vec<u8> = mant.bytes().filter(u8::is_ascii_digit).map(|b| b - b'0').collect();
|
||||
while s.len() > 1 && s.last() == Some(&0) {
|
||||
s.pop();
|
||||
}
|
||||
let n = exp.parse::<i32>().unwrap_or(0) + 1;
|
||||
// Rust's digits are the fewest, but a tie between two equally near it
|
||||
// rounds up. Only when another of the same length also reads back as x
|
||||
// (near 17 digits) is there a choice: then take x's exact expansion (a
|
||||
// double's has at most 767 significant digits) rounded to that length.
|
||||
if ![false, true].into_iter().any(|up| step(&s, up).is_some_and(|t| reads_as(&t, n, x))) {
|
||||
return (s, n);
|
||||
}
|
||||
let k = s.len();
|
||||
let exact = format!("{x:.800e}");
|
||||
let (em, ee) = exact.split_once('e').unwrap_or((&exact, "0"));
|
||||
let all: Vec<u8> = em.bytes().filter(u8::is_ascii_digit).map(|b| b - b'0').collect();
|
||||
let mut cn = ee.parse::<i32>().unwrap_or(0) + 1;
|
||||
let mut c = all[..k].to_vec();
|
||||
let rest = &all[k..];
|
||||
let up = match rest.first().copied().unwrap_or(0) {
|
||||
d if d > 5 => true,
|
||||
5 => rest[1..].iter().any(|&d| d != 0) || c[k - 1] % 2 == 1,
|
||||
_ => false,
|
||||
};
|
||||
if up {
|
||||
match step(&c, true) {
|
||||
Some(t) => c = t,
|
||||
None => {
|
||||
// 99…9 rounds to 10…0.
|
||||
c = std::iter::once(1).chain(std::iter::repeat_n(0, k - 1)).collect();
|
||||
cn += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
if reads_as(&c, cn, x) { (c, cn) } else { (s, n) }
|
||||
}
|
||||
|
||||
/// The k-digit decimal one unit above or below `s` in its last digit, if it
|
||||
/// still has k digits.
|
||||
fn step(s: &[u8], up: bool) -> Option<Vec<u8>> {
|
||||
let mut t = s.to_vec();
|
||||
for i in (0..t.len()).rev() {
|
||||
match (up, t[i]) {
|
||||
(true, 9) => t[i] = 0,
|
||||
(false, 0) => t[i] = 9,
|
||||
(true, d) | (false, d) => {
|
||||
t[i] = if up { d + 1 } else { d - 1 };
|
||||
return (t[0] != 0).then_some(t);
|
||||
}
|
||||
}
|
||||
}
|
||||
None
|
||||
}
|
||||
|
||||
/// Whether 0.s × 10^n reads as `x` (std's reader is correctly rounded).
|
||||
fn reads_as(s: &[u8], n: i32, x: f64) -> bool {
|
||||
let digits: String = s.iter().map(|d| char::from(b'0' + d)).collect();
|
||||
format!("0.{digits}e{n}").parse::<f64>().is_ok_and(|y| y == x)
|
||||
}
|
||||
|
||||
/// The first number in `v` that a double cannot hold exactly (an integer past
|
||||
/// ±2^53), by its path. Such a number would be read differently by readers that
|
||||
/// hold every number as a double, so the receipt could not bind it.
|
||||
pub fn inexact_number(v: &Value) -> Option<String> {
|
||||
const SAFE: u64 = (1 << 53) - 1;
|
||||
fn walk(v: &Value, path: &mut String) -> bool {
|
||||
match v {
|
||||
Value::Number(n) => n.as_u64().is_some_and(|u| u > SAFE) || n.as_i64().is_some_and(|i| i.unsigned_abs() > SAFE),
|
||||
Value::Array(a) => a.iter().enumerate().any(|(i, x)| {
|
||||
let at = path.len();
|
||||
path.push_str(&format!("[{i}]"));
|
||||
let hit = walk(x, path);
|
||||
if !hit {
|
||||
path.truncate(at);
|
||||
}
|
||||
hit
|
||||
}),
|
||||
Value::Object(m) => m.iter().any(|(k, x)| {
|
||||
let at = path.len();
|
||||
path.push('.');
|
||||
path.push_str(k);
|
||||
let hit = walk(x, path);
|
||||
if !hit {
|
||||
path.truncate(at);
|
||||
}
|
||||
hit
|
||||
}),
|
||||
_ => false,
|
||||
}
|
||||
}
|
||||
let mut path = String::new();
|
||||
walk(v, &mut path).then(|| path.trim_start_matches('.').to_string())
|
||||
}
|
||||
|
||||
impl Signer {
|
||||
/// The key at `path`, created (32 random bytes, readable by the owner
|
||||
/// only) when there is none.
|
||||
|
|
@ -70,7 +296,7 @@ impl Signer {
|
|||
pub fn seal(&self, request: &Value, result: &Value) -> Value {
|
||||
let d = digest(request, result);
|
||||
let sig = self.key.sign(&d);
|
||||
json!({ "alg": "ed25519", "key": self.public_key(), "digest": hex(&d), "signature": hex(&sig.to_bytes()) })
|
||||
json!({ "alg": "ed25519", "canon": CANON, "key": self.public_key(), "digest": hex(&d), "signature": hex(&sig.to_bytes()) })
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -81,7 +307,11 @@ pub fn verify(request: &Value, result: &Value, receipt: &Value) -> bool {
|
|||
let (Some(key), Some(sig)) = (field("key"), field("signature")) else { return false };
|
||||
let (Ok(key), Ok(sig)) = (<[u8; 32]>::try_from(key.as_slice()), <[u8; 64]>::try_from(sig.as_slice())) else { return false };
|
||||
let Ok(vk) = VerifyingKey::from_bytes(&key) else { return false };
|
||||
let d = digest(request, result);
|
||||
let d = match receipt.get("canon").and_then(Value::as_str) {
|
||||
Some(CANON) => digest(request, result),
|
||||
None => digest_v1(request, result),
|
||||
Some(_) => return false,
|
||||
};
|
||||
field("digest").is_some_and(|g| g == d) && vk.verify(&d, &wai_quantum::ed25519_dalek::Signature::from_bytes(&sig)).is_ok()
|
||||
}
|
||||
|
||||
|
|
@ -125,6 +355,86 @@ mod tests {
|
|||
assert!(!verify(&req, &res, &forged));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn numbers_are_written_as_ecmascript_writes_them() {
|
||||
// Bit patterns and ECMAScript's Number::toString of each, written by a
|
||||
// JavaScript engine: RFC 8785's appendix B, the edges of each layout,
|
||||
// and patterns spread over every exponent.
|
||||
let mut n = 0;
|
||||
for line in include_str!("../testdata/es-number-tostring.txt").lines() {
|
||||
let (bits, want) = line.split_once(' ').unwrap();
|
||||
let x = f64::from_bits(u64::from_str_radix(bits, 16).unwrap());
|
||||
assert_eq!(es_number(x), want, "{bits}");
|
||||
// And the text reads back as the same double, as it must for a
|
||||
// result to verify after a round trip (serde_json's
|
||||
// float_roundtrip: its default reader rounds twice).
|
||||
if x != 0.0 {
|
||||
assert_eq!(serde_json::from_str::<f64>(want).unwrap().to_bits(), x.to_bits(), "{bits} {want}");
|
||||
}
|
||||
n += 1;
|
||||
}
|
||||
assert!(n > 3000);
|
||||
}
|
||||
|
||||
/// The same against a larger set, in the same format, written by a
|
||||
/// JavaScript engine: `WAI_ES_VECTORS=path cargo test -- --ignored`.
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn numbers_match_a_javascript_engine_on_a_large_set() {
|
||||
let path = std::env::var("WAI_ES_VECTORS").expect("WAI_ES_VECTORS: a file of `hexbits text` lines");
|
||||
let text = std::fs::read_to_string(path).unwrap();
|
||||
let mut n = 0u64;
|
||||
for line in text.lines() {
|
||||
let (bits, want) = line.split_once(' ').unwrap();
|
||||
assert_eq!(es_number(f64::from_bits(u64::from_str_radix(bits, 16).unwrap())), want, "{bits}");
|
||||
n += 1;
|
||||
}
|
||||
eprintln!("{n} numbers agree");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_canonical_text_sorts_keys_and_writes_numbers_and_strings_one_way() {
|
||||
let v = json!({ "b": 1.0, "a": [true, null, "x\n\u{1}\"\\é"], "\u{e000}": 0.000001, "\u{1f600}": -0.0, "c": 1e21 });
|
||||
// U+1F600 is a surrogate pair, so it sorts before U+E000.
|
||||
assert_eq!(String::from_utf8(canonical(&v)).unwrap(), "{\"a\":[true,null,\"x\\n\\u0001\\\"\\\\é\"],\"b\":1,\"c\":1e+21,\"\u{1f600}\":0,\"\u{e000}\":0.000001}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_receipt_verifies_after_a_round_trip_through_another_json_writer() {
|
||||
let s = Signer::from_seed([7; 32]);
|
||||
let req = json!({ "qasm": "qubit[2] q; h q[0]; cx q[0], q[1];", "shots": 1000, "observables": ["Z0 Z1"] });
|
||||
let res = json!({ "method": "stabilizer", "expectations": [{ "observable": "Z0 Z1", "coeff": 1.0, "value": 1.0 }], "threshold": 0.00001, "elapsed_ms": 0 });
|
||||
let r = s.seal(&req, &res);
|
||||
// What a browser sends back: integral doubles without a fraction, keys
|
||||
// in another order, the threshold in plain decimals.
|
||||
let req2: Value = serde_json::from_str(r#"{"observables":["Z0 Z1"],"shots":1000,"qasm":"qubit[2] q; h q[0]; cx q[0], q[1];"}"#).unwrap();
|
||||
let res2: Value = serde_json::from_str(r#"{"elapsed_ms":0,"threshold":0.00001,"method":"stabilizer","expectations":[{"value":1,"coeff":1,"observable":"Z0 Z1"}]}"#).unwrap();
|
||||
assert!(verify(&req2, &res2, &r));
|
||||
let res3: Value = serde_json::from_str(r#"{"elapsed_ms":0,"threshold":0.00001,"method":"stabilizer","expectations":[{"value":0.999,"coeff":1,"observable":"Z0 Z1"}]}"#).unwrap();
|
||||
assert!(!verify(&req2, &res3, &r));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_first_receipts_still_verify_and_an_unknown_form_does_not() {
|
||||
let s = Signer::from_seed([7; 32]);
|
||||
let req = json!({ "qasm": "OPENQASM 3.0;", "shots": 10 });
|
||||
let res = json!({ "method": "statevector", "counts": { "0": 10 } });
|
||||
let d = digest_v1(&req, &res);
|
||||
let old = json!({ "alg": "ed25519", "key": s.public_key(), "digest": hex(&d), "signature": hex(&s.key.sign(&d).to_bytes()) });
|
||||
assert!(verify(&req, &res, &old));
|
||||
let mut other = s.seal(&req, &res);
|
||||
other["canon"] = json!("something-else");
|
||||
assert!(!verify(&req, &res, &other));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn integers_a_double_cannot_hold_are_found_by_path() {
|
||||
assert_eq!(inexact_number(&json!({ "seed": 9007199254740991u64, "shots": 10 })), None);
|
||||
assert_eq!(inexact_number(&json!({ "seed": 9007199254740992u64 })).as_deref(), Some("seed"));
|
||||
assert_eq!(inexact_number(&json!({ "options": { "x": [1, -9007199254740993i64] } })).as_deref(), Some("options.x[1]"));
|
||||
assert_eq!(inexact_number(&json!({ "noise": { "readout": 1e300 } })), None);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_key_is_created_once_private_and_then_reused() {
|
||||
let dir = std::env::temp_dir().join(format!("wai-qsvc-key-{}", std::process::id()));
|
||||
|
|
|
|||
|
|
@ -127,6 +127,36 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
|
|||
);
|
||||
}
|
||||
}
|
||||
"statevector_single" => {
|
||||
use wai_quantum::quantum_sv_metal::GpuState;
|
||||
let gates = r.program.matrices();
|
||||
// Wider fused gates run slower on the GPU than the passes they save.
|
||||
let width = r.options.get("fusion").and_then(Value::as_u64).unwrap_or(3).clamp(1, 3) as usize;
|
||||
let mut sv = GpuState::zero(n).map_err(|e| e.to_string())?;
|
||||
sv.run(&gates, width).map_err(|e| e.to_string())?;
|
||||
out["accuracy"] = json!({ "class": "single", "fusion": width, "device": sv.device(), "note": "single-precision amplitudes on the GPU, about 1e-7 relative rounding per gate; sums in double precision" });
|
||||
if r.shots > 0 {
|
||||
let mut map: HashMap<String, u64> = HashMap::new();
|
||||
for x in sv.sample(r.shots as usize, r.seed) {
|
||||
*map.entry(r.program.bitstring(|q| x >> (n - 1 - q) & 1 == 1)).or_insert(0) += 1;
|
||||
}
|
||||
merge(&mut out, counts(map));
|
||||
}
|
||||
if !r.observables.is_empty() {
|
||||
out["expectations"] = expectation_entries(&r, |s| sv.expectation(&[(1.0, s.to_vec())]));
|
||||
}
|
||||
if !r.amplitudes.is_empty() {
|
||||
out["amplitudes"] = Value::Array(
|
||||
r.amplitudes
|
||||
.iter()
|
||||
.map(|b| {
|
||||
let a = sv.amplitude(index_of(b, n, |q| n - 1 - q));
|
||||
json!({ "bits": b, "re": a.re, "im": a.im })
|
||||
})
|
||||
.collect(),
|
||||
);
|
||||
}
|
||||
}
|
||||
"statevector" => {
|
||||
let gates = r.program.matrices();
|
||||
let width = r.options.get("fusion").and_then(Value::as_u64).unwrap_or(3).clamp(1, 5) as usize;
|
||||
|
|
@ -155,6 +185,60 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
|
|||
);
|
||||
}
|
||||
}
|
||||
"stabilizer" if r.noise.is_some() => {
|
||||
use wai_quantum::quantum_tableau::{FrameStep, frame_counts, frame_expectations};
|
||||
let noise = r.noise.ok_or("no noise model")?;
|
||||
// Each op's gates, then its noise on the qubits it touched, as the
|
||||
// other noisy methods apply it.
|
||||
let mut steps = Vec::new();
|
||||
for (gates, qs) in r.program.clifford_ops()? {
|
||||
steps.extend(gates.into_iter().map(FrameStep::Gate));
|
||||
if qs.len() == 2 && noise.depol2 > 0.0 {
|
||||
steps.push(FrameStep::Depolarize2(qs[0], qs[1], noise.depol2));
|
||||
}
|
||||
for &q in &qs {
|
||||
if qs.len() != 2 && noise.depol1 > 0.0 {
|
||||
steps.push(FrameStep::Depolarize1(q, noise.depol1));
|
||||
}
|
||||
if noise.dephasing > 0.0 {
|
||||
steps.push(FrameStep::Dephase(q, noise.dephasing));
|
||||
}
|
||||
}
|
||||
}
|
||||
let threads = std::thread::available_parallelism().map_or(1, |t| t.get()).min(8);
|
||||
let frame_shots = catalog::frame_expectation_shots(&r);
|
||||
out["accuracy"] = json!({
|
||||
"class": "estimate",
|
||||
"shots": frame_shots,
|
||||
"note": "Pauli frames against a noiseless stabilizer-tableau reference: counts are samples with the law of independent noisy shots; expectations are the noiseless values times the mean sign of the frames, with standard errors"
|
||||
});
|
||||
if r.shots > 0 {
|
||||
let measured: Vec<u32> = r.program.measured.iter().map(|&(q, _)| q).collect();
|
||||
let nb = r.program.n_bits as usize;
|
||||
let mut map: HashMap<String, u64> = HashMap::new();
|
||||
for (bits, k) in frame_counts(n, &steps, &measured, noise.readout, r.shots, r.seed, threads).map_err(|e| e.to_string())? {
|
||||
let mut s = vec![b'0'; nb];
|
||||
for (j, &(_, b)) in r.program.measured.iter().enumerate() {
|
||||
if bits[j / 64] >> (j % 64) & 1 == 1 {
|
||||
s[nb - 1 - b as usize] = b'1';
|
||||
}
|
||||
}
|
||||
*map.entry(String::from_utf8(s).expect("ASCII")).or_insert(0) += k;
|
||||
}
|
||||
merge(&mut out, counts(map));
|
||||
}
|
||||
if !r.observables.is_empty() {
|
||||
let strings: Vec<Vec<(u32, char)>> = r.observables.iter().map(|(_, _, s)| s.clone()).collect();
|
||||
let est = frame_expectations(n, &steps, &strings, frame_shots, r.seed, threads).map_err(|e| e.to_string())?;
|
||||
out["expectations"] = Value::Array(
|
||||
r.observables
|
||||
.iter()
|
||||
.zip(est)
|
||||
.map(|((text, coeff, _), (v, se))| json!({ "observable": text, "coeff": coeff, "value": coeff * v, "stderr": coeff.abs() * se }))
|
||||
.collect(),
|
||||
);
|
||||
}
|
||||
}
|
||||
"stabilizer" => {
|
||||
use wai_quantum::quantum_tableau::Tableau;
|
||||
let t = Tableau::run(n, &r.program.clifford()?).map_err(|e| e.to_string())?;
|
||||
|
|
@ -187,8 +271,16 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
|
|||
}
|
||||
}
|
||||
"pauli" => {
|
||||
use wai_quantum::quantum_spd::{Pauli, SpdConfig, expectation};
|
||||
use wai_quantum::quantum_spd::{Pauli, SpdConfig, expectation_capped};
|
||||
let rc = r.program.rotations();
|
||||
// As many terms as the job's memory holds: a term's string (two
|
||||
// words of bits a 64 qubits, padded to a power of two), its
|
||||
// coefficient and index slots, and room for a rotation's new terms.
|
||||
let words = (n as u64).div_ceil(64).next_power_of_two();
|
||||
let per_term = 2 * (16 * words + 8) + 16;
|
||||
let fits = ((plan.memory_bytes.saturating_sub(256 << 20)) / per_term).max(1 << 16) as usize;
|
||||
let max_terms = r.options.get("max_terms").and_then(Value::as_u64).map_or(fits, |m| (m as usize).clamp(1, fits));
|
||||
let threads = std::thread::available_parallelism().map_or(1, |t| t.get());
|
||||
let threshold = r.options.get("threshold").and_then(Value::as_f64).unwrap_or(1e-5).clamp(0.0, 1.0);
|
||||
let max_weight = r.options.get("max_weight").and_then(Value::as_u64).map(|w| w as u32);
|
||||
// The ladder: the same value at up to three thresholds a decade
|
||||
|
|
@ -202,12 +294,22 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
|
|||
let mut steps = Vec::new();
|
||||
let mut last = None;
|
||||
for &t in &rungs {
|
||||
let res = expectation(&rc, &[(*coeff, obs.clone())], &[], &SpdConfig { threshold: t, max_weight }).map_err(|e| e.to_string())?;
|
||||
steps.push(json!({ "threshold": t, "value": res.value, "bound": res.truncation_bound, "peak_terms": res.peak_terms }));
|
||||
last = Some(res);
|
||||
let (res, used) = expectation_capped(&rc, &[(*coeff, obs.clone())], &[], &SpdConfig { threshold: t, max_weight }, threads, max_terms).map_err(|e| e.to_string())?;
|
||||
let mut step = json!({ "threshold": t, "value": res.value, "bound": res.truncation_bound, "peak_terms": res.peak_terms });
|
||||
if used > t {
|
||||
step["threshold_reached"] = json!(used);
|
||||
}
|
||||
steps.push(step);
|
||||
last = Some((res, used));
|
||||
}
|
||||
let res = last.expect("one rung at least");
|
||||
let (res, used) = last.expect("one rung at least");
|
||||
let mut entry = json!({ "observable": text, "coeff": coeff, "value": res.value, "bound": res.truncation_bound, "peak_terms": res.peak_terms });
|
||||
if used > threshold {
|
||||
// The sum outgrew the job's memory: the threshold rose to
|
||||
// fit, and the bound covers what that dropped.
|
||||
entry["threshold_reached"] = json!(used);
|
||||
entry["max_terms"] = json!(max_terms);
|
||||
}
|
||||
if steps.len() > 1 {
|
||||
entry["ladder"] = Value::Array(steps);
|
||||
}
|
||||
|
|
@ -277,7 +379,7 @@ pub fn run(request: &Value, plan: &Plan) -> Result<Value, String> {
|
|||
"density_matrix" => {
|
||||
use wai_quantum::quantum_noisy::DensityMatrix;
|
||||
let noise = r.noise.ok_or("no noise model")?;
|
||||
let rho = DensityMatrix::run(&r.program.matrices(), n, &noise).map_err(|e| format!("{e:?}"))?;
|
||||
let rho = DensityMatrix::evolve(&r.program.matrices(), n, &noise).map_err(|e| format!("{e:?}"))?;
|
||||
out["accuracy"] = json!({ "class": "exact", "note": "the noisy state itself, under the model; counts are seeded samples of its distribution" });
|
||||
if r.shots > 0 {
|
||||
let p = rho.probabilities(noise.readout);
|
||||
|
|
@ -378,6 +480,22 @@ mod tests {
|
|||
run(&job, &p).unwrap()
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn grover_as_the_ide_exports_it_finds_the_marked_item() {
|
||||
// Its oracle and its diffusion are each a Z under four controls,
|
||||
// written `ctrl(4) @ z`. Four rounds on five qubits find the marked
|
||||
// item with probability 0.9992.
|
||||
let marked = 0b10110u32;
|
||||
let qasm = wai_quantum::quantum_toolchain::to_qasm(&wai_quantum::quantum_toolchain::grover(5, marked));
|
||||
assert!(qasm.contains("ctrl(4) @ z"), "{qasm}");
|
||||
let key: String = (0..5).rev().map(|q| if marked >> q & 1 == 1 { '1' } else { '0' }).collect();
|
||||
for method in ["exact", "statevector"] {
|
||||
let r = go(json!({ "qasm": qasm, "method": method, "shots": 1000, "seed": 3 }));
|
||||
let hits = r["counts"][&key].as_u64().unwrap_or(0);
|
||||
assert!(hits > 980, "{method}: {r}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_bell_pair_is_the_same_on_every_method_that_takes_it() {
|
||||
let qasm = "OPENQASM 3.0; qubit[2] q; bit[2] c; h q[0]; cx q[0], q[1]; c = measure q;";
|
||||
|
|
@ -455,6 +573,80 @@ mod tests {
|
|||
assert!(histories < shots / 4, "{histories} histories for {shots} shots");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noisy_clifford_circuits_run_on_frames_and_agree_with_the_density_matrix() {
|
||||
let qasm = "qubit[4] q; bit[4] c; h q[0]; cx q[0], q[1]; s q[1]; cx q[1], q[2]; h q[3]; cz q[3], q[2]; rx(pi/2) q[0]; c = measure q;";
|
||||
let noise = json!({ "depolarizing_1q": 0.01, "depolarizing_2q": 0.03, "dephasing": 0.01, "readout": 0.02 });
|
||||
let p = plan(&json!({ "qasm": qasm, "noise": noise, "shots": 100 }), CAP).unwrap();
|
||||
// Small enough for the exact density matrix, which auto prefers.
|
||||
assert_eq!(p.method, "density_matrix");
|
||||
let shots = 100_000u64;
|
||||
let dm = go(json!({ "qasm": qasm, "method": "density_matrix", "noise": noise, "shots": 1_000_000, "observables": ["Z0 Z1", "X3 Z2", "Y1"] }));
|
||||
let st = go(json!({ "qasm": qasm, "method": "stabilizer", "noise": noise, "shots": shots, "seed": 5, "observables": ["Z0 Z1", "X3 Z2", "Y1"] }));
|
||||
assert_eq!(st["method"], "stabilizer", "{st}");
|
||||
assert_eq!(st["accuracy"]["class"], "estimate");
|
||||
for (bits, c) in st["counts"].as_object().unwrap() {
|
||||
let p = dm["counts"][bits].as_u64().unwrap_or(0) as f64 / 1e6;
|
||||
let f = c.as_u64().unwrap() as f64 / shots as f64;
|
||||
let se = (p * (1.0 - p) / shots as f64).sqrt().max(1e-4);
|
||||
assert!((f - p).abs() < 5.0 * se + 2e-3, "{bits}: {f} vs {p}");
|
||||
}
|
||||
for k in 0..3 {
|
||||
let (v, se) = (st["expectations"][k]["value"].as_f64().unwrap(), st["expectations"][k]["stderr"].as_f64().unwrap());
|
||||
let want = dm["expectations"][k]["value"].as_f64().unwrap();
|
||||
assert!((v - want).abs() < 5.0 * se + 1e-9, "observable {k}: {v} ± {se} vs {want}");
|
||||
}
|
||||
// Wide: auto takes the frames for a noisy Clifford circuit past the density matrix.
|
||||
let wide = format!("qubit[300] q; h q[0]; {} ", (0..299).map(|i| format!("cx q[{i}], q[{}];", i + 1)).collect::<String>());
|
||||
let w = go(json!({ "qasm": wide, "noise": { "depolarizing_2q": 0.001 }, "shots": 1000, "observables": ["Z0 Z299"] }));
|
||||
assert_eq!(w["method"], "stabilizer", "{w}");
|
||||
let v = w["expectations"][0]["value"].as_f64().unwrap();
|
||||
assert!(v > 0.5 && v < 1.0, "{v}");
|
||||
// Amplitude damping is not a Pauli channel.
|
||||
let e = plan(&json!({ "qasm": wide, "method": "stabilizer", "noise": { "amplitude_damping": 0.01 } }), CAP).unwrap_err();
|
||||
assert!(e.contains("amplitude damping"), "{e}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_gpu_state_vector_agrees_where_there_is_one_and_is_refused_where_not() {
|
||||
let qasm = "qubit[14] q; bit[14] c; h q; rx(0.3) q[0]; cx q[0], q[13]; crz(0.7) q[3], q[9]; ry(1.1) q[5]; c = measure q;";
|
||||
let job = json!({ "qasm": qasm, "method": "statevector_single", "shots": 2000, "observables": ["Z0 Z13", "X5"], "amplitudes": ["00000000000000"] });
|
||||
match crate::catalog::gpu_max_qubits() {
|
||||
Ok(_) => {
|
||||
let g = go(job);
|
||||
let d = go(json!({ "qasm": qasm, "method": "statevector", "observables": ["Z0 Z13", "X5"], "amplitudes": ["00000000000000"] }));
|
||||
assert_eq!(g["accuracy"]["class"], "single");
|
||||
for k in 0..2 {
|
||||
let (a, b) = (g["expectations"][k]["value"].as_f64().unwrap(), d["expectations"][k]["value"].as_f64().unwrap());
|
||||
assert!((a - b).abs() < 1e-5, "{a} vs {b}");
|
||||
}
|
||||
let (a, b) = (g["amplitudes"][0]["re"].as_f64().unwrap(), d["amplitudes"][0]["re"].as_f64().unwrap());
|
||||
assert!((a - b).abs() < 1e-6);
|
||||
assert_eq!(g["counts"].as_object().unwrap().values().map(|v| v.as_u64().unwrap()).sum::<u64>(), 2000);
|
||||
}
|
||||
Err(_) => {
|
||||
let e = plan(&job, CAP).unwrap_err();
|
||||
assert!(e.contains("needs a GPU"), "{e}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_pauli_sum_held_to_a_term_cap_keeps_its_bound() {
|
||||
let qasm = "qubit[8] q; rx(0.7) q; rzz(0.4) q[0], q[1]; rzz(0.4) q[2], q[3]; rzz(0.4) q[4], q[5]; rzz(0.4) q[6], q[7]; ry(0.9) q; rzz(0.3) q[1], q[2]; rzz(0.3) q[3], q[4]; rzz(0.3) q[5], q[6]; rx(1.1) q; rzz(0.5) q[0], q[7]; ry(0.4) q;";
|
||||
let exact = go(json!({ "qasm": qasm, "method": "statevector", "observables": ["Z0 Z5"] }));
|
||||
let want = exact["expectations"][0]["value"].as_f64().unwrap();
|
||||
let capped = go(json!({ "qasm": qasm, "method": "pauli", "observables": ["Z0 Z5"], "options": { "threshold": 0.0, "max_terms": 40 } }));
|
||||
let e = &capped["expectations"][0];
|
||||
assert!(e["threshold_reached"].as_f64().unwrap() > 0.0, "{e}");
|
||||
let (v, b) = (e["value"].as_f64().unwrap(), e["bound"].as_f64().unwrap());
|
||||
assert!((v - want).abs() <= b + 1e-9, "{v} vs {want}, bound {b}");
|
||||
// Uncapped at threshold 0 it is exact.
|
||||
let free = go(json!({ "qasm": qasm, "method": "pauli", "observables": ["Z0 Z5"], "options": { "threshold": 0.0 } }));
|
||||
assert!(free["expectations"][0].get("threshold_reached").is_none());
|
||||
assert!((free["expectations"][0]["value"].as_f64().unwrap() - want).abs() < 1e-9);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn pauli_answers_carry_a_threshold_ladder() {
|
||||
let qasm = "qubit[6] q; rx(0.3) q[0]; rzz(0.4) q[0], q[1]; ry(0.9) q[1]; rzz(0.2) q[1], q[2]; rx(1.1) q[2]; cx q[2], q[3]; ry(0.5) q[3];";
|
||||
|
|
|
|||
3385
wai/wai-quantum-service/testdata/es-number-tostring.txt
vendored
Normal file
3385
wai/wai-quantum-service/testdata/es-number-tostring.txt
vendored
Normal file
File diff suppressed because it is too large
Load diff
2
wai/wai-quantum/Cargo.lock
generated
2
wai/wai-quantum/Cargo.lock
generated
|
|
@ -1079,7 +1079,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"block2",
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
[package]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
edition = "2024"
|
||||
rust-version = "1.97"
|
||||
license = "Apache-2.0"
|
||||
|
|
@ -178,6 +178,8 @@ quantum_fabric_fpga = ["quantum_fabric"]
|
|||
# The neural-processor member: exact fp16 limb products on the platform's neural processor.
|
||||
quantum_fabric_npu = ["quantum_fabric", "dep:objc2", "dep:block2", "dep:objc2-foundation", "dep:objc2-core-ml"]
|
||||
quantum_fabric_metal = ["quantum_fabric", "dep:objc2", "dep:objc2-foundation", "dep:objc2-metal"]
|
||||
# The dense state vector in single precision on the platform's GPU (macOS).
|
||||
quantum_sv_metal = ["quantum_sv", "dep:objc2", "dep:objc2-foundation", "dep:objc2-metal"]
|
||||
# The npu-cores member: the gate's arithmetic on a neural processor's compute core, programmed directly (Linux): the program built by the tile array's open toolchain, the host side compiled against the device runtime, both at run time.
|
||||
quantum_fabric_npu_cores = ["quantum_fabric"]
|
||||
# The tensor-engine member: the gate's products as exact bf16 limb convolutions on a tensor processor's matrix engine, through its kernel toolchain's Python run as a helper.
|
||||
|
|
|
|||
|
|
@ -62,7 +62,8 @@ pub mod gates;
|
|||
feature = "quantum_arch",
|
||||
feature = "quantum_qio",
|
||||
feature = "quantum_dmrg",
|
||||
feature = "quantum_tdvp"
|
||||
feature = "quantum_tdvp",
|
||||
feature = "quantum_tableau"
|
||||
))]
|
||||
mod repro;
|
||||
|
||||
|
|
@ -250,6 +251,9 @@ pub mod quantum_receipt;
|
|||
pub mod quantum_stabilizer;
|
||||
#[cfg(feature = "quantum_tableau")]
|
||||
pub mod quantum_tableau;
|
||||
|
||||
#[cfg(feature = "quantum_sv_metal")]
|
||||
pub mod quantum_sv_metal;
|
||||
#[cfg(feature = "quantum_mpsim")]
|
||||
pub mod quantum_mpsim;
|
||||
|
||||
|
|
|
|||
|
|
@ -731,9 +731,13 @@ pub fn linear_xeb(ideal: &StateVector, samples: &[u64]) -> Estimate {
|
|||
// The exact density matrix
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
/// Largest register the exact density matrix accepts.
|
||||
/// Largest register the exact density matrix accepts in [`DensityMatrix::run`].
|
||||
pub const MAX_DENSITY_QUBITS: u32 = 12;
|
||||
|
||||
/// Largest register [`DensityMatrix::evolve`] accepts: `ρ` is then a state on
|
||||
/// 32 qubits, 64 GiB.
|
||||
pub const MAX_EVOLVE_QUBITS: u32 = 16;
|
||||
|
||||
/// `ρ` on `n` qubits, stored as a state on `2n` qubits: entry `(r, c)` at
|
||||
/// index `r·2ⁿ + c`.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
|
|
@ -821,6 +825,50 @@ impl DensityMatrix {
|
|||
Ok(rho)
|
||||
}
|
||||
|
||||
/// Run `gates` under `noise` from `|0…0⟩⟨0…0|`, exactly, with each gate
|
||||
/// and its noise as one superoperator on `ρ`'s vector: the gate on the
|
||||
/// rows, its conjugate on the columns, and each channel's
|
||||
/// `Σ_i w_i K_i ⊗ K_i*` on its qubits' row and column bits, in the order
|
||||
/// [`run`](Self::run) applies them. They are fused (to at most four of the
|
||||
/// vector's qubits for gates on one or two) and applied in the state
|
||||
/// vector's tiled passes, so a gate costs one pass and no copy of `ρ`.
|
||||
/// The result is `run`'s to rounding, on up to [`MAX_EVOLVE_QUBITS`].
|
||||
pub fn evolve(gates: &[Gate], n: u32, noise: &NoiseModel) -> Result<DensityMatrix, SvError> {
|
||||
if n > MAX_EVOLVE_QUBITS {
|
||||
return Err(SvError::TooLarge);
|
||||
}
|
||||
let mut list = Vec::with_capacity(gates.len() * 4);
|
||||
for (i, g) in gates.iter().enumerate() {
|
||||
if g.qubits.iter().any(|&q| q >= n) {
|
||||
return Err(SvError::BadGate(i));
|
||||
}
|
||||
list.push(g.clone());
|
||||
list.push(conj(g, n));
|
||||
let k = g.qubits.len();
|
||||
if k == 2 && noise.depol2 > 0.0 {
|
||||
let (a, b) = (g.qubits[0], g.qubits[1]);
|
||||
let terms: Vec<(f64, Gate)> = (0..16u8).map(|j| (if j == 0 { 1.0 - noise.depol2 } else { noise.depol2 / 15.0 }, two_qubit(&pauli(a, j >> 2), &pauli(b, j & 3)))).collect();
|
||||
list.push(superoperator(n, &terms));
|
||||
}
|
||||
for &q in &g.qubits {
|
||||
if k != 2 && noise.depol1 > 0.0 {
|
||||
let terms: Vec<(f64, Gate)> = (0..4u8).map(|j| (if j == 0 { 1.0 - noise.depol1 } else { noise.depol1 / 3.0 }, pauli(q, j))).collect();
|
||||
list.push(superoperator(n, &terms));
|
||||
}
|
||||
if noise.dephasing > 0.0 {
|
||||
list.push(superoperator(n, &[(1.0 - noise.dephasing, pauli(q, 0)), (noise.dephasing, pauli(q, 3))]));
|
||||
}
|
||||
if noise.damping > 0.0 {
|
||||
let [k0, k1] = damping_kraus(q, noise.damping);
|
||||
list.push(superoperator(n, &[(1.0, k0), (1.0, k1)]));
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut vec = StateVector::zero(2 * n)?;
|
||||
vec.run(&list, 4)?;
|
||||
Ok(DensityMatrix { n, vec })
|
||||
}
|
||||
|
||||
/// Entry `ρ_{r,c}`.
|
||||
pub fn entry(&self, r: usize, c: usize) -> C {
|
||||
self.vec.amplitude((r << self.n) | c)
|
||||
|
|
@ -888,6 +936,41 @@ impl DensityMatrix {
|
|||
}
|
||||
}
|
||||
|
||||
/// A channel `ρ → Σ_i w_i K_i ρ K_i†` on `k` qubits as a gate on the vector
|
||||
/// of `ρ`: `Σ_i w_i K_i ⊗ K_i*` on the qubits' row bits then their column
|
||||
/// bits (`q` and `q + n`).
|
||||
fn superoperator(n: u32, terms: &[(f64, Gate)]) -> Gate {
|
||||
let qs = &terms[0].1.qubits;
|
||||
let k = qs.len();
|
||||
let d = 1usize << k;
|
||||
let big = d * d;
|
||||
let mut m = vec![complex(0.0, 0.0); big * big];
|
||||
for (w, kr) in terms {
|
||||
debug_assert_eq!(&kr.qubits, qs);
|
||||
for r in 0..d {
|
||||
for c in 0..d {
|
||||
for r2 in 0..d {
|
||||
let a = kr.matrix[r * d + r2];
|
||||
if a.re == 0.0 && a.im == 0.0 {
|
||||
continue;
|
||||
}
|
||||
for c2 in 0..d {
|
||||
let b = kr.matrix[c * d + c2];
|
||||
// a · conj(b) · w
|
||||
let (pr, pi) = (a.re * b.re + a.im * b.im, a.im * b.re - a.re * b.im);
|
||||
let e = &mut m[((r << k) | c) * big + ((r2 << k) | c2)];
|
||||
e.re += w * pr;
|
||||
e.im += w * pi;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut qubits = qs.clone();
|
||||
qubits.extend(qs.iter().map(|&q| q + n));
|
||||
Gate { qubits, matrix: m }
|
||||
}
|
||||
|
||||
/// The tensor product `a ⊗ b` of two one-qubit gates on different qubits.
|
||||
fn two_qubit(a: &Gate, b: &Gate) -> Gate {
|
||||
let mut m = vec![complex(0.0, 0.0); 16];
|
||||
|
|
@ -1017,6 +1100,68 @@ mod tests {
|
|||
assert_eq!(a, b);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn evolving_by_superoperators_is_the_density_matrix() {
|
||||
for (n, seed) in [(3u32, 1u64), (4, 2), (6, 3)] {
|
||||
let mut g = random_circuit(n, 4, seed);
|
||||
g.push(Gate::cx(0, n - 1));
|
||||
g.push(Gate { qubits: vec![1, 0, n - 1], matrix: {
|
||||
// A Toffoli on three qubits: dense, with its channels per qubit.
|
||||
let mut m = vec![complex(0.0, 0.0); 64];
|
||||
for i in 0..8 {
|
||||
m[i * 8 + if i >= 6 { i ^ 1 } else { i }] = complex(1.0, 0.0);
|
||||
}
|
||||
m
|
||||
} });
|
||||
let noise = NoiseModel { depol1: 0.02, depol2: 0.04, damping: 0.03, dephasing: 0.01, readout: 0.0 };
|
||||
let want = DensityMatrix::run(&g, n, &noise).unwrap();
|
||||
let got = DensityMatrix::evolve(&g, n, &noise).unwrap();
|
||||
let worst = (0..1usize << n).flat_map(|r| (0..1usize << n).map(move |c| (r, c))).map(|(r, c)| {
|
||||
let (a, b) = (got.entry(r, c), want.entry(r, c));
|
||||
(a.re - b.re).abs().max((a.im - b.im).abs())
|
||||
}).fold(0.0, f64::max);
|
||||
assert!(worst < 1e-12, "n {n}: {worst}");
|
||||
assert!((got.trace() - 1.0).abs() < 1e-12);
|
||||
}
|
||||
}
|
||||
|
||||
/// Times; run with `--ignored --nocapture`.
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn density_timings() {
|
||||
let noise = NoiseModel { depol1: 0.01, depol2: 0.03, damping: 0.02, dephasing: 0.01, readout: 0.0 };
|
||||
for n in [10u32, 12, 14, 15] {
|
||||
let g = random_circuit(n, 6, 2);
|
||||
let t0 = std::time::Instant::now();
|
||||
let rho = DensityMatrix::evolve(&g, n, &noise).unwrap();
|
||||
let te = t0.elapsed();
|
||||
let old = if n <= MAX_DENSITY_QUBITS {
|
||||
let t1 = std::time::Instant::now();
|
||||
DensityMatrix::run(&g, n, &noise).unwrap();
|
||||
format!("{:?}", t1.elapsed())
|
||||
} else {
|
||||
"-".into()
|
||||
};
|
||||
println!("n={n}: {} gates; evolve {te:?}, run {old}; trace {:.12}", g.len(), rho.trace());
|
||||
}
|
||||
}
|
||||
|
||||
/// Past the referee's size, against branching trajectories.
|
||||
#[test]
|
||||
fn a_fourteen_qubit_density_matrix_agrees_with_trajectories() {
|
||||
let n = 14;
|
||||
let g = random_circuit(n, 3, 21);
|
||||
let noise = NoiseModel { depol1: 0.01, depol2: 0.03, damping: 0.02, dephasing: 0.01, readout: 0.0 };
|
||||
let rho = DensityMatrix::evolve(&g, n, &noise).unwrap();
|
||||
assert!((rho.trace() - 1.0).abs() < 1e-10);
|
||||
let obs = observables(n);
|
||||
let tr = branching_trajectories(&g, n, &noise, &obs, 4000, 3, 8, 8).unwrap();
|
||||
for (k, o) in obs.iter().enumerate() {
|
||||
let (exact, e) = (rho.expectation(o), tr.expectations[k]);
|
||||
assert!((e.mean - exact).abs() < 4.0 * e.stderr + 1e-9, "observable {k}: {} ± {} vs {exact}", e.mean, e.stderr);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn the_binomial_sampler_has_the_binomial_moments() {
|
||||
for (c, p, draws) in [(1000u64, 0.003f64, 20_000usize), (50, 0.7, 20_000), (6, 0.2, 20_000), (1_000_000, 1e-4, 4_000)] {
|
||||
|
|
|
|||
|
|
@ -822,6 +822,22 @@ pub fn expectation_threaded(
|
|||
cfg: &SpdConfig,
|
||||
threads: usize,
|
||||
) -> Result<SpdResult, SpdError> {
|
||||
expectation_capped(circuit, observable, ones, cfg, threads, usize::MAX).map(|(r, _)| r)
|
||||
}
|
||||
|
||||
/// [`expectation_threaded`] holding at most `max_terms` terms: whenever a
|
||||
/// rotation leaves more, the threshold is doubled and the sum truncated again
|
||||
/// until it fits, and every coefficient dropped that way is added to the
|
||||
/// bound like any other. A run that never reaches the cap is the uncapped
|
||||
/// run, bit for bit. Returns the result and the threshold the run ended at.
|
||||
pub fn expectation_capped(
|
||||
circuit: &RotCircuit,
|
||||
observable: &[(f64, PauliString)],
|
||||
ones: &[u32],
|
||||
cfg: &SpdConfig,
|
||||
threads: usize,
|
||||
max_terms: usize,
|
||||
) -> Result<(SpdResult, f64), SpdError> {
|
||||
circuit.validate()?;
|
||||
for (_, s) in observable {
|
||||
check_string(s, circuit.n)?;
|
||||
|
|
@ -832,11 +848,11 @@ pub fn expectation_threaded(
|
|||
}
|
||||
}
|
||||
Ok(match circuit.n.div_ceil(64) {
|
||||
1 => run::<1>(circuit, observable, ones, cfg, threads),
|
||||
2 => run::<2>(circuit, observable, ones, cfg, threads),
|
||||
3 | 4 => run::<4>(circuit, observable, ones, cfg, threads),
|
||||
5..=8 => run::<8>(circuit, observable, ones, cfg, threads),
|
||||
_ => run::<16>(circuit, observable, ones, cfg, threads),
|
||||
1 => run::<1>(circuit, observable, ones, cfg, threads, max_terms),
|
||||
2 => run::<2>(circuit, observable, ones, cfg, threads, max_terms),
|
||||
3 | 4 => run::<4>(circuit, observable, ones, cfg, threads, max_terms),
|
||||
5..=8 => run::<8>(circuit, observable, ones, cfg, threads, max_terms),
|
||||
_ => run::<16>(circuit, observable, ones, cfg, threads, max_terms),
|
||||
})
|
||||
}
|
||||
|
||||
|
|
@ -851,11 +867,12 @@ fn run<const W: usize>(
|
|||
ones: &[u32],
|
||||
cfg: &SpdConfig,
|
||||
threads: usize,
|
||||
) -> SpdResult {
|
||||
crate::pool::scoped(threads, |pool| propagate::<W>(circuit, observable, ones, cfg, pool))
|
||||
max_terms: usize,
|
||||
) -> (SpdResult, f64) {
|
||||
crate::pool::scoped(threads, |pool| propagate::<W>(circuit, observable, ones, cfg, pool, max_terms))
|
||||
}
|
||||
|
||||
fn propagate<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliString)], ones: &[u32], cfg: &SpdConfig, pool: &Pool) -> SpdResult {
|
||||
fn propagate<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliString)], ones: &[u32], cfg: &SpdConfig, pool: &Pool, max_terms: usize) -> (SpdResult, f64) {
|
||||
let mut t = Terms::<W>::new();
|
||||
for (c, s) in observable {
|
||||
t.add(Key::from_string(s), *c, pool);
|
||||
|
|
@ -868,6 +885,7 @@ fn propagate<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliStri
|
|||
// before the next rotation and at the end. A Clifford gate keeps the
|
||||
// number of terms, so the work count and the peak are as gate by gate.
|
||||
let mut batch: Vec<Clifford<'_, W>> = Vec::new();
|
||||
let mut threshold = cfg.threshold;
|
||||
for op in circuit.ops.iter().rev() {
|
||||
work += t.live() as u64;
|
||||
match op {
|
||||
|
|
@ -878,7 +896,12 @@ fn propagate<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliStri
|
|||
}
|
||||
let (s, c) = sin_cos(*theta);
|
||||
rotate(&mut t, &Key::from_string(axis), s, c, &mut pending, pool);
|
||||
bound += t.truncate(cfg.threshold, cfg.max_weight, pool);
|
||||
bound += t.truncate(threshold, cfg.max_weight, pool);
|
||||
// Over the cap: a doubled threshold until the sum fits.
|
||||
while t.live() > max_terms {
|
||||
threshold = if threshold > 0.0 { threshold * 2.0 } else { f64::MIN_POSITIVE };
|
||||
bound += t.truncate(threshold, cfg.max_weight, pool);
|
||||
}
|
||||
}
|
||||
Op::QuarterRot { axis, k } => batch.push(Clifford::Quarter(Key::from_string(axis), *k)),
|
||||
other => batch.push(Clifford::Gate(other)),
|
||||
|
|
@ -897,7 +920,7 @@ fn propagate<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliStri
|
|||
value += if flips & 1 == 1 { -t.coef[i] } else { t.coef[i] };
|
||||
}
|
||||
}
|
||||
SpdResult { value, truncation_bound: bound, peak_terms: peak, final_terms: t.live(), term_updates: work }
|
||||
(SpdResult { value, truncation_bound: bound, peak_terms: peak, final_terms: t.live(), term_updates: work }, threshold)
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
|
|
@ -1458,6 +1481,30 @@ mod tests {
|
|||
assert_eq!((one.value.to_bits(), one.truncation_bound.to_bits(), one.peak_terms), (many.value.to_bits(), many.truncation_bound.to_bits(), many.peak_terms));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_term_cap_raises_the_threshold_and_the_bound_still_holds() {
|
||||
let ki = kicked_ising(127, &heavy_hex_127(), 0.7, 5);
|
||||
let obs = [(1.0, vec![(62, Pauli::Z)])];
|
||||
let cfg = SpdConfig { threshold: 1e-5, max_weight: None };
|
||||
let (free, t0) = expectation_capped(&ki, &obs, &[], &cfg, 4, usize::MAX).unwrap();
|
||||
assert_eq!(t0, 1e-5);
|
||||
// A cap above the peak changes nothing.
|
||||
let (same, t1) = expectation_capped(&ki, &obs, &[], &cfg, 4, free.peak_terms + 1).unwrap();
|
||||
assert_eq!((same.value.to_bits(), same.truncation_bound.to_bits(), t1), (free.value.to_bits(), free.truncation_bound.to_bits(), 1e-5));
|
||||
// A cap a tenth of the peak raises the threshold; the capped value
|
||||
// stays within its own bound of the far more exact one.
|
||||
let (capped, t2) = expectation_capped(&ki, &obs, &[], &cfg, 4, free.peak_terms / 10).unwrap();
|
||||
assert!(t2 > 1e-5 && capped.peak_terms <= free.peak_terms);
|
||||
assert!((capped.value - free.value).abs() <= capped.truncation_bound + free.truncation_bound + 1e-12);
|
||||
// A small register, exact reference: the bound covers the true error.
|
||||
let mut rng = Rng(5);
|
||||
let c = random_circuit(&mut rng, 6, 60);
|
||||
let o = vec![(1.0, vec![(0u32, Pauli::Z), (3, Pauli::X)])];
|
||||
let exact = dense_expect(&c, &o[0].1, &[]);
|
||||
let (r, _) = expectation_capped(&c, &o, &[], &SpdConfig { threshold: 0.0, max_weight: None }, 2, 40).unwrap();
|
||||
assert!((r.value - exact).abs() <= r.truncation_bound + 1e-12, "{} vs {exact}, bound {}", r.value, r.truncation_bound);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn word_width_does_not_change_the_answer() {
|
||||
// The same circuit, placed in a 30-qubit and a 900-qubit register, is
|
||||
|
|
|
|||
|
|
@ -134,7 +134,7 @@ impl Schedule {
|
|||
/// One pass of a run: a gate applied to the whole state, or a group of gates
|
||||
/// applied tile by tile, with the local bits their dense gates act on.
|
||||
#[derive(Clone, Debug, PartialEq, Eq)]
|
||||
enum Stage {
|
||||
pub(crate) enum Stage {
|
||||
Direct(usize),
|
||||
Tiled { gates: Vec<usize>, local: Vec<u32> },
|
||||
}
|
||||
|
|
@ -144,7 +144,7 @@ enum Stage {
|
|||
/// `tile − low`. Diagonal gates join any pass. A gate on more bits than a tile
|
||||
/// holds gets a pass of its own. The order of the gates never changes, so
|
||||
/// neither does any amplitude's arithmetic.
|
||||
fn plan(list: &[Gate], n: u32, tile: u32, low: u32) -> Vec<Stage> {
|
||||
pub(crate) fn plan(list: &[Gate], n: u32, tile: u32, low: u32) -> Vec<Stage> {
|
||||
let room = (tile - low) as usize;
|
||||
let mut out = Vec::new();
|
||||
let (mut gates, mut local): (Vec<usize>, Vec<u32>) = (Vec::new(), Vec::new());
|
||||
|
|
@ -290,6 +290,7 @@ thread_local! {
|
|||
/// Run `f` with this thread's state-vector kernels on this thread alone. For
|
||||
/// callers that run many small states in parallel; the result is the same
|
||||
/// bits either way.
|
||||
#[cfg_attr(not(feature = "quantum_noisy"), allow(dead_code))]
|
||||
pub(crate) fn on_one_thread<R>(f: impl FnOnce() -> R) -> R {
|
||||
let before = ONE_THREAD.with(|o| o.replace(true));
|
||||
let r = f();
|
||||
|
|
@ -369,9 +370,81 @@ fn kernel_lanes(re: &mut [f64], im: &mut [f64], m: &[C], offsets: &[usize], sort
|
|||
}
|
||||
}
|
||||
|
||||
/// A gate on `log₂ D` qubits, `L` groups side by side: each group's members
|
||||
/// are loaded, from one run of `L` neighbours where the lowest target stride
|
||||
/// allows, gathered otherwise, and every output is the [`row_dot`] sum, term
|
||||
/// for term, over the matrix split into its real and imaginary parts. The
|
||||
/// sizes are compile-time constants, so the loops unroll and the members stay
|
||||
/// in registers; the result is the scalar path's, bit for bit.
|
||||
fn kernel_fixed<const D: usize, const L: usize>(re: &mut [f64], im: &mut [f64], m: &[C], offsets: &[usize], sorted: &[u32]) {
|
||||
let mut mr = [[0.0f64; D]; D];
|
||||
let mut mi = [[0.0f64; D]; D];
|
||||
for r in 0..D {
|
||||
for c in 0..D {
|
||||
mr[r][c] = m[r * D + c].re;
|
||||
mi[r][c] = m[r * D + c].im;
|
||||
}
|
||||
}
|
||||
let mut off = [0usize; D];
|
||||
off.copy_from_slice(offsets);
|
||||
let top = off.iter().copied().max().unwrap_or(0);
|
||||
let len = re.len();
|
||||
assert!(im.len() == len);
|
||||
let k = D.trailing_zeros() as usize;
|
||||
let groups = len >> k;
|
||||
let run = 1usize << sorted[0];
|
||||
let mut t = 0;
|
||||
while t < groups {
|
||||
let lanes = L.min(groups - t);
|
||||
let mut bases = [0usize; L];
|
||||
if run >= L && lanes == L {
|
||||
let b0 = spread(t, sorted);
|
||||
for (l, b) in bases.iter_mut().enumerate() {
|
||||
*b = b0 + l;
|
||||
}
|
||||
} else {
|
||||
for (l, b) in bases.iter_mut().enumerate().take(lanes) {
|
||||
*b = spread(t + l, sorted);
|
||||
}
|
||||
}
|
||||
// Every index read or written below is at most this.
|
||||
for &b in &bases[..lanes] {
|
||||
assert!(b + top < len);
|
||||
}
|
||||
let (mut vr, mut vi) = ([[0.0f64; L]; D], [[0.0f64; L]; D]);
|
||||
for j in 0..D {
|
||||
for l in 0..lanes {
|
||||
// SAFETY: `bases[l] + off[j] ≤ bases[l] + top < len`, checked above.
|
||||
unsafe {
|
||||
vr[j][l] = *re.get_unchecked(bases[l] + off[j]);
|
||||
vi[j][l] = *im.get_unchecked(bases[l] + off[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
for r in 0..D {
|
||||
let (mut sr, mut si) = ([0.0f64; L], [0.0f64; L]);
|
||||
for c in 0..D {
|
||||
let (a, b) = (mr[r][c], mi[r][c]);
|
||||
for l in 0..L {
|
||||
sr[l] += a * vr[c][l] - b * vi[c][l];
|
||||
si[l] += a * vi[c][l] + b * vr[c][l];
|
||||
}
|
||||
}
|
||||
for l in 0..lanes {
|
||||
// SAFETY: as for the loads.
|
||||
unsafe {
|
||||
*re.get_unchecked_mut(bases[l] + off[r]) = sr[l];
|
||||
*im.get_unchecked_mut(bases[l] + off[r]) = si[l];
|
||||
}
|
||||
}
|
||||
}
|
||||
t += lanes;
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply a gate to a slice holding whole groups (block-aligned). One-qubit
|
||||
/// gates run over the two halves of each block side by side; wider gates take
|
||||
/// [`kernel_lanes`] when their groups come in runs, unrolled loops otherwise.
|
||||
/// gates run over the two halves of each block side by side; gates on two to
|
||||
/// five qubits take [`kernel_fixed`]; wider ones the general loop.
|
||||
fn kernel_block(re: &mut [f64], im: &mut [f64], m: &[C], offsets: &[usize], sorted: &[u32]) {
|
||||
let k = sorted.len();
|
||||
if k == 1 {
|
||||
|
|
@ -391,6 +464,13 @@ fn kernel_block(re: &mut [f64], im: &mut [f64], m: &[C], offsets: &[usize], sort
|
|||
}
|
||||
return;
|
||||
}
|
||||
match k {
|
||||
2 => return kernel_fixed::<4, 4>(re, im, m, offsets, sorted),
|
||||
3 => return kernel_fixed::<8, 4>(re, im, m, offsets, sorted),
|
||||
4 => return kernel_fixed::<16, 2>(re, im, m, offsets, sorted),
|
||||
5 => return kernel_fixed::<32, 2>(re, im, m, offsets, sorted),
|
||||
_ => {}
|
||||
}
|
||||
if 1usize << sorted[0] >= LANES {
|
||||
kernel_lanes(re, im, m, offsets, sorted);
|
||||
return;
|
||||
|
|
|
|||
1030
wai/wai-quantum/src/quantum_sv_metal.rs
Normal file
1030
wai/wai-quantum/src/quantum_sv_metal.rs
Normal file
File diff suppressed because it is too large
Load diff
|
|
@ -24,6 +24,23 @@
|
|||
//!
|
||||
//! Signs follow Aaronson and Gottesman's `g` (the power of `i` when two Pauli
|
||||
//! letters multiply), evaluated a word at a time.
|
||||
//!
|
||||
//! **Noisy circuits** ([`frame_counts`], [`frame_expectations`]). Under Pauli
|
||||
//! noise (depolarizing, dephasing, readout flips) a Clifford circuit is
|
||||
//! sampled by Pauli frames, 64 shots to a machine word:
|
||||
//! - One noiseless outcome from the tableau is the reference, and each shot
|
||||
//! carries only the difference its noise makes.
|
||||
//! - Frames start with random Z bits. Those leave `|0…0⟩` as it is, and
|
||||
//! carried through the circuit they spread each shot uniformly over the
|
||||
//! noiseless outcomes.
|
||||
//! - Errors are placed by geometric gaps, so a channel costs its error rate,
|
||||
//! not a draw a shot.
|
||||
//!
|
||||
//! An expectation value is the noiseless one times the frames' mean sign
|
||||
//! against the string. Both readouts are the same bits at any thread count,
|
||||
//! and agree with the exact density matrix on small circuits. On a 16-core
|
||||
//! desktop processor, 100,000 full-width shots of a 10,000-qubit circuit with
|
||||
//! 100,000 noisy gates take 1.6 s; at 1,000 qubits and 10,000 gates, 44 ms.
|
||||
|
||||
/// The most qubits a tableau holds: the tableau is `4n²` bits, so 32768
|
||||
/// qubits take 512 MiB, and the outcome subspace about as much again.
|
||||
|
|
@ -58,6 +75,8 @@ pub enum TableauError {
|
|||
RepeatedQubit(u32),
|
||||
/// A Pauli factor that is not `X`, `Y`, `Z` or `I`.
|
||||
BadPauli(char),
|
||||
/// A noise probability outside `[0, 1]` (as written).
|
||||
BadProbability(String),
|
||||
}
|
||||
|
||||
impl std::fmt::Display for TableauError {
|
||||
|
|
@ -67,6 +86,7 @@ impl std::fmt::Display for TableauError {
|
|||
TableauError::QubitOutOfRange { qubit, n } => write!(f, "qubit {qubit} of a {n}-qubit tableau"),
|
||||
TableauError::RepeatedQubit(q) => write!(f, "qubit {q} named twice"),
|
||||
TableauError::BadPauli(c) => write!(f, "{c:?} is not a Pauli factor"),
|
||||
TableauError::BadProbability(p) => write!(f, "{p} is not a probability"),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -567,6 +587,373 @@ impl Sampler {
|
|||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Noisy shots by Pauli frames
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
/// One step of a noisy Clifford circuit: a gate, or a Pauli channel.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub enum FrameStep {
|
||||
Gate(Clifford),
|
||||
/// X, Y or Z on the qubit, each with probability `p/3`.
|
||||
Depolarize1(u32, f64),
|
||||
/// One of the fifteen non-identity Paulis on the pair, each `p/15`.
|
||||
Depolarize2(u32, u32, f64),
|
||||
/// Z on the qubit with probability `p`.
|
||||
Dephase(u32, f64),
|
||||
}
|
||||
|
||||
/// Words of shots a qubit holds in one batch (64 shots a word).
|
||||
const FRAME_WORDS: usize = 128;
|
||||
|
||||
/// The gates of `steps`, checked, without the noise.
|
||||
fn frame_gates(n: u32, steps: &[FrameStep]) -> Result<Vec<Clifford>, TableauError> {
|
||||
let qubit = |q: u32| if q < n { Ok(()) } else { Err(TableauError::QubitOutOfRange { qubit: q, n }) };
|
||||
let prob = |p: f64| if (0.0..=1.0).contains(&p) { Ok(()) } else { Err(TableauError::BadProbability(p.to_string())) };
|
||||
let mut gates = Vec::new();
|
||||
for s in steps {
|
||||
match *s {
|
||||
FrameStep::Gate(g) => gates.push(g),
|
||||
FrameStep::Depolarize1(q, p) | FrameStep::Dephase(q, p) => {
|
||||
qubit(q)?;
|
||||
prob(p)?;
|
||||
}
|
||||
FrameStep::Depolarize2(a, b, p) => {
|
||||
qubit(a)?;
|
||||
qubit(b)?;
|
||||
if a == b {
|
||||
return Err(TableauError::RepeatedQubit(a));
|
||||
}
|
||||
prob(p)?;
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(gates)
|
||||
}
|
||||
|
||||
/// The next success, at least one trial after the last, in trials of
|
||||
/// probability `p`; the logarithm is the crate's own, so the draws are the
|
||||
/// same everywhere.
|
||||
fn gap(p: f64, rng: &mut u64) -> u64 {
|
||||
if p >= 1.0 {
|
||||
return 1;
|
||||
}
|
||||
let u = 1.0 - (splitmix(rng) >> 11) as f64 * (1.0 / 9_007_199_254_740_992.0);
|
||||
let g = (crate::repro::ln(u) / crate::repro::ln(1.0 - p)).floor();
|
||||
if g >= 1.8e19 { u64::MAX } else { g as u64 + 1 }
|
||||
}
|
||||
|
||||
/// The Pauli frames of one batch of shots: for qubit `q`, its X bits are
|
||||
/// `x[q·w..(q+1)·w]` and its Z bits `z[…]`, one bit a shot.
|
||||
struct Frames {
|
||||
w: usize,
|
||||
x: Vec<u64>,
|
||||
z: Vec<u64>,
|
||||
}
|
||||
|
||||
impl Frames {
|
||||
/// Frames for `n` qubits with every Z bit random: a Z on `|0⟩` changes
|
||||
/// nothing, and carried through the circuit it makes each shot's
|
||||
/// noiseless outcome a uniform draw from the outcome subspace.
|
||||
fn gauge(n: u32, w: usize, rng: &mut u64) -> Frames {
|
||||
let len = n as usize * w;
|
||||
Frames { w, x: vec![0; len], z: (0..len).map(|_| splitmix(rng)).collect() }
|
||||
}
|
||||
|
||||
fn run(&self, q: u32) -> core::ops::Range<usize> {
|
||||
let a = q as usize * self.w;
|
||||
a..a + self.w
|
||||
}
|
||||
|
||||
/// Carry the frames through a gate: `P ↦ G·P·G†`, signs dropped.
|
||||
fn gate(&mut self, g: Clifford) {
|
||||
let w = self.w;
|
||||
let (x, z) = (&mut self.x, &mut self.z);
|
||||
match g {
|
||||
Clifford::X(_) | Clifford::Y(_) | Clifford::Z(_) => {}
|
||||
Clifford::H(q) => {
|
||||
let r = q as usize * w..(q as usize + 1) * w;
|
||||
x[r.clone()].swap_with_slice(&mut z[r]);
|
||||
}
|
||||
Clifford::S(q) | Clifford::Sdg(q) => {
|
||||
let a = q as usize * w;
|
||||
for i in a..a + w {
|
||||
z[i] ^= x[i];
|
||||
}
|
||||
}
|
||||
Clifford::SX(q) | Clifford::SXdg(q) => {
|
||||
let a = q as usize * w;
|
||||
for i in a..a + w {
|
||||
x[i] ^= z[i];
|
||||
}
|
||||
}
|
||||
Clifford::CX(c, t) => {
|
||||
let (c, t) = (c as usize * w, t as usize * w);
|
||||
for i in 0..w {
|
||||
x[t + i] ^= x[c + i];
|
||||
z[c + i] ^= z[t + i];
|
||||
}
|
||||
}
|
||||
Clifford::CY(c, t) => {
|
||||
let (c, t) = (c as usize * w, t as usize * w);
|
||||
for i in 0..w {
|
||||
let (xt, zt) = (x[t + i], z[t + i]);
|
||||
z[c + i] ^= xt ^ zt;
|
||||
x[t + i] ^= x[c + i];
|
||||
z[t + i] ^= x[c + i];
|
||||
}
|
||||
}
|
||||
Clifford::CZ(a, b) => {
|
||||
let (a, b) = (a as usize * w, b as usize * w);
|
||||
for i in 0..w {
|
||||
z[a + i] ^= x[b + i];
|
||||
z[b + i] ^= x[a + i];
|
||||
}
|
||||
}
|
||||
Clifford::Swap(a, b) => {
|
||||
let (a, b) = (a as usize * w, b as usize * w);
|
||||
for i in 0..w {
|
||||
x.swap(a + i, b + i);
|
||||
z.swap(a + i, b + i);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Multiply shot `s`'s frame by Pauli `p` (1 X, 2 Y, 3 Z) on `q`.
|
||||
fn flip(&mut self, q: u32, p: u64, s: usize) {
|
||||
let (i, b) = (q as usize * self.w + s / 64, 1u64 << (s % 64));
|
||||
if p != 3 {
|
||||
self.x[i] ^= b;
|
||||
}
|
||||
if p >= 2 {
|
||||
self.z[i] ^= b;
|
||||
}
|
||||
}
|
||||
|
||||
/// Draw a channel's errors over the batch's shots and apply them.
|
||||
fn channel(&mut self, step: FrameStep, rng: &mut u64) {
|
||||
let shots = (64 * self.w) as u64;
|
||||
let p = match step {
|
||||
FrameStep::Gate(_) => return,
|
||||
FrameStep::Depolarize1(_, p) | FrameStep::Depolarize2(_, _, p) | FrameStep::Dephase(_, p) => p,
|
||||
};
|
||||
if p <= 0.0 {
|
||||
return;
|
||||
}
|
||||
let mut at = 0u64;
|
||||
loop {
|
||||
at = at.saturating_add(gap(p, rng));
|
||||
if at > shots {
|
||||
return;
|
||||
}
|
||||
let s = (at - 1) as usize;
|
||||
match step {
|
||||
FrameStep::Depolarize1(q, _) => self.flip(q, 1 + splitmix(rng) % 3, s),
|
||||
FrameStep::Depolarize2(a, b, _) => {
|
||||
let k = 1 + splitmix(rng) % 15;
|
||||
if k >> 2 != 0 {
|
||||
self.flip(a, k >> 2, s);
|
||||
}
|
||||
if k & 3 != 0 {
|
||||
self.flip(b, k & 3, s);
|
||||
}
|
||||
}
|
||||
FrameStep::Dephase(q, _) => self.flip(q, 3, s),
|
||||
FrameStep::Gate(_) => unreachable!(),
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Run `f(batch)` for every batch of `shots`, on up to `threads` threads, and
|
||||
/// return the results in batch order. A batch's stream depends only on
|
||||
/// `(seed, batch)`.
|
||||
fn frame_batches<R: Send>(shots: u64, w: usize, threads: usize, f: &(dyn Fn(u64, usize) -> R + Sync)) -> Vec<R> {
|
||||
let per = 64 * w as u64;
|
||||
let batches = shots.div_ceil(per);
|
||||
let size = |b: u64| ((shots - b * per).min(per)) as usize;
|
||||
let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.clamp(1, batches.max(1) as usize) };
|
||||
if threads == 1 {
|
||||
return (0..batches).map(|b| f(b, size(b))).collect();
|
||||
}
|
||||
// A word-sized counter: some targets have no 64-bit atomics.
|
||||
let next = std::sync::atomic::AtomicUsize::new(0);
|
||||
let slots: Vec<std::sync::Mutex<Option<R>>> = (0..batches).map(|_| std::sync::Mutex::new(None)).collect();
|
||||
std::thread::scope(|sc| {
|
||||
for _ in 0..threads {
|
||||
sc.spawn(|| loop {
|
||||
let i = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
|
||||
if i as u64 >= batches {
|
||||
break;
|
||||
}
|
||||
let b = i as u64;
|
||||
*slots[i].lock().unwrap() = Some(f(b, size(b)));
|
||||
});
|
||||
}
|
||||
});
|
||||
slots.into_iter().map(|m| m.into_inner().unwrap().expect("every batch ran")).collect()
|
||||
}
|
||||
|
||||
fn batch_rng(seed: u64, batch: u64) -> u64 {
|
||||
let mut s = seed ^ batch.wrapping_mul(0xD1B5_4A32_D192_ED03) ^ 0xF4A3_E5A5_0001;
|
||||
splitmix(&mut s);
|
||||
s
|
||||
}
|
||||
|
||||
/// Words a qubit holds per batch for `shots` shots.
|
||||
fn frame_words(shots: u64) -> usize {
|
||||
(shots.div_ceil(64) as usize).clamp(1, FRAME_WORDS)
|
||||
}
|
||||
|
||||
/// `shots` measurement outcomes of a noisy Clifford circuit from `|0…0⟩`, the
|
||||
/// `measured` qubits read in the computational basis with each bit flipped
|
||||
/// with probability `readout`, counted as [`Sampler::counts`] counts them.
|
||||
///
|
||||
/// Pauli-frame sampling: one noiseless outcome from the tableau is the
|
||||
/// reference, and each shot carries a Pauli frame (the difference its noise
|
||||
/// makes), 64 shots to a machine word. Frames start with random Z bits,
|
||||
/// which leave `|0…0⟩` unchanged and, carried through the circuit, spread each
|
||||
/// shot uniformly over the noiseless outcomes. Noise is drawn by geometric
|
||||
/// gaps between errors, so a channel costs its error rate, not a draw a shot.
|
||||
/// A shot's outcome is the reference plus its frame's X bits on the measured
|
||||
/// qubits. The counts have the law of independent noisy shots, and are the
|
||||
/// same bits at any thread count.
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
pub fn frame_counts(n: u32, steps: &[FrameStep], measured: &[u32], readout: f64, shots: u64, seed: u64, threads: usize) -> Result<Vec<(Vec<u64>, u64)>, TableauError> {
|
||||
let gates = frame_gates(n, steps)?;
|
||||
if !(0.0..=1.0).contains(&readout) {
|
||||
return Err(TableauError::BadProbability(readout.to_string()));
|
||||
}
|
||||
let t = Tableau::run(n, &gates)?;
|
||||
for &q in measured {
|
||||
if q >= n {
|
||||
return Err(TableauError::QubitOutOfRange { qubit: q, n });
|
||||
}
|
||||
}
|
||||
let reference = t.support().sampler(measured).x0;
|
||||
let m = measured.len();
|
||||
let mw = m.div_ceil(64);
|
||||
let w = frame_words(shots);
|
||||
let maps = frame_batches(shots, w, threads, &|b, size| {
|
||||
let mut rng = batch_rng(seed, b);
|
||||
let mut f = Frames::gauge(n, w, &mut rng);
|
||||
for &s in steps {
|
||||
match s {
|
||||
FrameStep::Gate(g) => f.gate(g),
|
||||
other => f.channel(other, &mut rng),
|
||||
}
|
||||
}
|
||||
// Each shot's outcome bits, measured qubit j at bit j.
|
||||
let mut out = vec![0u64; size * mw];
|
||||
for (j, &q) in measured.iter().enumerate() {
|
||||
let r = if reference[j / 64] >> (j % 64) & 1 == 1 { !0u64 } else { 0 };
|
||||
let xs = &f.x[f.run(q)];
|
||||
for (k, &word) in xs.iter().enumerate() {
|
||||
let mut bits = word ^ r;
|
||||
while bits != 0 {
|
||||
let s = k * 64 + bits.trailing_zeros() as usize;
|
||||
bits &= bits - 1;
|
||||
if s < size {
|
||||
out[s * mw + j / 64] |= 1 << (j % 64);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if readout > 0.0 && m > 0 {
|
||||
let trials = (size * m) as u64;
|
||||
let mut at = 0u64;
|
||||
loop {
|
||||
at = at.saturating_add(gap(readout, &mut rng));
|
||||
if at > trials {
|
||||
break;
|
||||
}
|
||||
let k = (at - 1) as usize;
|
||||
let (s, j) = (k / m, k % m);
|
||||
out[s * mw + j / 64] ^= 1 << (j % 64);
|
||||
}
|
||||
}
|
||||
let mut map: std::collections::HashMap<Vec<u64>, u64> = std::collections::HashMap::new();
|
||||
for s in 0..size {
|
||||
*map.entry(out[s * mw..(s + 1) * mw].to_vec()).or_insert(0) += 1;
|
||||
}
|
||||
map
|
||||
});
|
||||
let mut total: std::collections::HashMap<Vec<u64>, u64> = std::collections::HashMap::new();
|
||||
for map in maps {
|
||||
for (k, c) in map {
|
||||
*total.entry(k).or_insert(0) += c;
|
||||
}
|
||||
}
|
||||
let mut v: Vec<(Vec<u64>, u64)> = total.into_iter().collect();
|
||||
v.sort_by(|a, b| b.1.cmp(&a.1).then(a.0.cmp(&b.0)));
|
||||
Ok(v)
|
||||
}
|
||||
|
||||
/// Each Pauli string's expectation value in the noisy state, estimated from
|
||||
/// `shots` frames, with its standard error.
|
||||
///
|
||||
/// Under Pauli noise the state is a mixture of the noiseless state with Pauli
|
||||
/// errors applied, so `⟨P⟩` is the noiseless `⟨P⟩` (`±1` or `0`, from the
|
||||
/// tableau) times the average sign of `P` against the shots' frames. A string
|
||||
/// whose noiseless value is `0` stays `0` exactly.
|
||||
pub fn frame_expectations(n: u32, steps: &[FrameStep], strings: &[Vec<(u32, char)>], shots: u64, seed: u64, threads: usize) -> Result<Vec<(f64, f64)>, TableauError> {
|
||||
let gates = frame_gates(n, steps)?;
|
||||
let t = Tableau::run(n, &gates)?;
|
||||
let ideal: Vec<i8> = strings.iter().map(|s| t.expectation(s)).collect::<Result<_, _>>()?;
|
||||
let shots = shots.max(1);
|
||||
let w = frame_words(shots);
|
||||
// Per batch, each string's count of shots whose frame anticommutes with it.
|
||||
let flips = frame_batches(shots, w, threads, &|b, size| {
|
||||
let mut rng = batch_rng(seed, b);
|
||||
let mut f = Frames::gauge(n, w, &mut rng);
|
||||
for &s in steps {
|
||||
match s {
|
||||
FrameStep::Gate(g) => f.gate(g),
|
||||
other => f.channel(other, &mut rng),
|
||||
}
|
||||
}
|
||||
strings
|
||||
.iter()
|
||||
.zip(&ideal)
|
||||
.map(|(s, &e)| {
|
||||
if e == 0 {
|
||||
return 0u64;
|
||||
}
|
||||
let mut count = 0u64;
|
||||
for k in 0..w {
|
||||
let mut par = 0u64;
|
||||
for &(q, c) in s {
|
||||
let i = q as usize * w + k;
|
||||
par ^= match c.to_ascii_uppercase() {
|
||||
'X' => f.z[i],
|
||||
'Z' => f.x[i],
|
||||
'Y' => f.x[i] ^ f.z[i],
|
||||
_ => 0,
|
||||
};
|
||||
}
|
||||
let valid = if (k + 1) * 64 <= size { !0u64 } else if k * 64 >= size { 0 } else { (1u64 << (size - k * 64)) - 1 };
|
||||
count += u64::from((par & valid).count_ones());
|
||||
}
|
||||
count
|
||||
})
|
||||
.collect::<Vec<u64>>()
|
||||
});
|
||||
Ok(ideal
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(i, &e)| {
|
||||
if e == 0 {
|
||||
return (0.0, 0.0);
|
||||
}
|
||||
let flipped: u64 = flips.iter().map(|v| v[i]).sum();
|
||||
let mean = 1.0 - 2.0 * flipped as f64 / shots as f64;
|
||||
let var = if shots > 1 { (1.0 - mean * mean) * shots as f64 / (shots - 1) as f64 } else { 0.0 };
|
||||
(f64::from(e) * mean, (var / shots as f64).sqrt())
|
||||
})
|
||||
.collect())
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
|
@ -749,3 +1136,183 @@ mod tests {
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(all(test, feature = "quantum_noisy"))]
|
||||
mod frame_tests {
|
||||
use super::*;
|
||||
use crate::gates::{complex, Gate};
|
||||
use crate::quantum_noisy::{DensityMatrix, NoiseModel};
|
||||
|
||||
fn next(s: &mut u64) -> u64 {
|
||||
splitmix(s)
|
||||
}
|
||||
|
||||
fn gate(c: Clifford) -> Gate {
|
||||
let (z, o, i) = (complex(0.0, 0.0), complex(1.0, 0.0), complex(0.0, 1.0));
|
||||
match c {
|
||||
Clifford::H(q) => Gate::h(q),
|
||||
Clifford::S(q) => Gate::s(q),
|
||||
Clifford::Sdg(q) => Gate::sdg(q),
|
||||
Clifford::X(q) => Gate::x(q),
|
||||
Clifford::Y(q) => Gate::y(q),
|
||||
Clifford::Z(q) => Gate::z(q),
|
||||
Clifford::SX(q) => Gate::sqrt_x(q),
|
||||
Clifford::SXdg(q) => {
|
||||
let (a, b) = (complex(0.5, -0.5), complex(0.5, 0.5));
|
||||
Gate { qubits: vec![q], matrix: vec![a, b, b, a] }
|
||||
}
|
||||
Clifford::CX(c, t) => Gate::cx(c, t),
|
||||
Clifford::CY(c, t) => Gate { qubits: vec![c, t], matrix: vec![o, z, z, z, z, o, z, z, z, z, z, complex(0.0, -1.0), z, z, i, z] },
|
||||
Clifford::CZ(a, b) => Gate::cz(a, b),
|
||||
Clifford::Swap(a, b) => Gate::swap(a, b),
|
||||
}
|
||||
}
|
||||
|
||||
fn random_clifford(n: u32, len: usize, seed: u64) -> Vec<Clifford> {
|
||||
let mut s = seed;
|
||||
(0..len)
|
||||
.map(|_| {
|
||||
let a = (next(&mut s) % u64::from(n)) as u32;
|
||||
let b = (a + 1 + (next(&mut s) % u64::from(n - 1)) as u32) % n;
|
||||
match next(&mut s) % 12 {
|
||||
0 => Clifford::H(a),
|
||||
1 => Clifford::S(a),
|
||||
2 => Clifford::Sdg(a),
|
||||
3 => Clifford::SX(a),
|
||||
4 => Clifford::SXdg(a),
|
||||
5 => Clifford::X(a),
|
||||
6 => Clifford::Y(a),
|
||||
7 => Clifford::CX(a, b),
|
||||
8 => Clifford::CY(a, b),
|
||||
9 => Clifford::CZ(a, b),
|
||||
10 => Clifford::Swap(a, b),
|
||||
_ => Clifford::H(b),
|
||||
}
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// The steps the other noisy methods apply: each gate, then its
|
||||
/// depolarizing channel, then dephasing on each qubit it touched.
|
||||
fn noisy(gates: &[Clifford], noise: &NoiseModel) -> Vec<FrameStep> {
|
||||
let mut out = Vec::new();
|
||||
for &g in gates {
|
||||
out.push(FrameStep::Gate(g));
|
||||
let qs = gate(g).qubits;
|
||||
if qs.len() == 2 {
|
||||
out.push(FrameStep::Depolarize2(qs[0], qs[1], noise.depol2));
|
||||
}
|
||||
for &q in &qs {
|
||||
if qs.len() != 2 {
|
||||
out.push(FrameStep::Depolarize1(q, noise.depol1));
|
||||
}
|
||||
out.push(FrameStep::Dephase(q, noise.dephasing));
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noisy_counts_follow_the_density_matrix() {
|
||||
for (n, seed) in [(3u32, 1u64), (4, 2), (5, 3)] {
|
||||
let gates = random_clifford(n, 30, seed);
|
||||
let noise = NoiseModel { depol1: 0.02, depol2: 0.05, dephasing: 0.01, damping: 0.0, readout: 0.03 };
|
||||
let matrices: Vec<Gate> = gates.iter().map(|&g| gate(g)).collect();
|
||||
let p = DensityMatrix::run(&matrices, n, &noise).unwrap().probabilities(noise.readout);
|
||||
let measured: Vec<u32> = (0..n).collect();
|
||||
let shots = 200_000u64;
|
||||
let counts = frame_counts(n, &noisy(&gates, &noise), &measured, noise.readout, shots, seed, 4).unwrap();
|
||||
assert_eq!(counts.iter().map(|c| c.1).sum::<u64>(), shots);
|
||||
let mut got = vec![0u64; 1 << n];
|
||||
for (bits, c) in &counts {
|
||||
// Bit j is qubit j; the density matrix puts qubit q at bit n − 1 − q.
|
||||
let x = (0..n).fold(0usize, |x, q| x | (((bits[0] >> q) & 1) as usize) << (n - 1 - q));
|
||||
got[x] = *c;
|
||||
}
|
||||
for x in 0..1usize << n {
|
||||
let want = p[x] * shots as f64;
|
||||
let sd = (want * (1.0 - p[x])).sqrt().max(1.0);
|
||||
assert!((got[x] as f64 - want).abs() < 5.0 * sd, "n {n} x {x}: {} vs {want:.0}", got[x]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn noisy_expectations_follow_the_density_matrix() {
|
||||
for (n, seed) in [(4u32, 7u64), (5, 8)] {
|
||||
let gates = random_clifford(n, 40, seed);
|
||||
let noise = NoiseModel { depol1: 0.01, depol2: 0.04, dephasing: 0.02, damping: 0.0, readout: 0.0 };
|
||||
let matrices: Vec<Gate> = gates.iter().map(|&g| gate(g)).collect();
|
||||
let rho = DensityMatrix::run(&matrices, n, &noise).unwrap();
|
||||
// The noiseless stabilizers, so that the values are not all zero.
|
||||
let t = Tableau::run(n, &gates).unwrap();
|
||||
let mut s = seed;
|
||||
let mut strings: Vec<Vec<(u32, char)>> = Vec::new();
|
||||
while strings.len() < 12 {
|
||||
let st: Vec<(u32, char)> = (0..n).filter_map(|q| match next(&mut s) % 4 { 0 => None, 1 => Some((q, 'X')), 2 => Some((q, 'Y')), _ => Some((q, 'Z')) }).collect();
|
||||
if !st.is_empty() && (t.expectation(&st).unwrap() != 0 || strings.len().is_multiple_of(3)) {
|
||||
strings.push(st);
|
||||
}
|
||||
}
|
||||
let est = frame_expectations(n, &noisy(&gates, &noise), &strings, 100_000, seed, 3).unwrap();
|
||||
for (st, &(v, se)) in strings.iter().zip(&est) {
|
||||
let exact = rho.expectation(&[(1.0, st.clone())]);
|
||||
assert!((v - exact).abs() < 5.0 * se + 1e-9, "n {n} {st:?}: {v} ± {se} vs {exact}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn frames_are_the_same_bits_at_any_thread_count() {
|
||||
let n = 30;
|
||||
let gates = random_clifford(n, 300, 5);
|
||||
let noise = NoiseModel { depol1: 0.01, depol2: 0.02, dephasing: 0.01, damping: 0.0, readout: 0.02 };
|
||||
let steps = noisy(&gates, &noise);
|
||||
let measured: Vec<u32> = (0..n).step_by(3).collect();
|
||||
// A shot count that leaves a partial word and a partial batch.
|
||||
let one = frame_counts(n, &steps, &measured, 0.02, 20_011, 9, 1).unwrap();
|
||||
assert_eq!(one.iter().map(|c| c.1).sum::<u64>(), 20_011);
|
||||
assert_eq!(frame_counts(n, &steps, &measured, 0.02, 20_011, 9, 5).unwrap(), one);
|
||||
let strings = vec![vec![(0, 'Z'), (3, 'X')], vec![(7, 'Y')]];
|
||||
assert_eq!(frame_expectations(n, &steps, &strings, 9_999, 2, 1).unwrap(), frame_expectations(n, &steps, &strings, 9_999, 2, 4).unwrap());
|
||||
}
|
||||
|
||||
/// Times at scale; run with `--ignored --nocapture`.
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn frame_timings_at_scale() {
|
||||
let noise = NoiseModel { depol1: 1e-4, depol2: 1e-3, dephasing: 1e-4, damping: 0.0, readout: 1e-3 };
|
||||
for (n, len, shots) in [(1000u32, 10_000usize, 100_000u64), (10_000, 100_000, 100_000)] {
|
||||
let gates = random_clifford(n, len, 11);
|
||||
let steps = noisy(&gates, &noise);
|
||||
let measured: Vec<u32> = (0..n).collect();
|
||||
let t0 = std::time::Instant::now();
|
||||
let counts = frame_counts(n, &steps, &measured, noise.readout, shots, 1, 8).unwrap();
|
||||
let tc = t0.elapsed();
|
||||
let t1 = std::time::Instant::now();
|
||||
let e = frame_expectations(n, &steps, &[vec![(0, 'Z')], vec![(1, 'X'), (2, 'Z')]], 1_000_000, 2, 8).unwrap();
|
||||
println!("n={n} gates={len} steps={}: {shots} full-width shots in {tc:?} ({} distinct); 2 expectations from 1e6 shots in {:?}: {e:?}", steps.len(), counts.len(), t1.elapsed());
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_wide_noisy_ghz_state_mostly_agrees() {
|
||||
let n = 2000;
|
||||
let mut gates = vec![Clifford::H(0)];
|
||||
gates.extend((0..n - 1).map(|q| Clifford::CX(q, q + 1)));
|
||||
let noise = NoiseModel { depol1: 0.0, depol2: 1e-4, dephasing: 0.0, damping: 0.0, readout: 0.0 };
|
||||
let steps = noisy(&gates, &noise);
|
||||
let measured: Vec<u32> = (0..n).collect();
|
||||
let counts = frame_counts(n, &steps, &measured, 0.0, 2000, 1, 4).unwrap();
|
||||
let words = n.div_ceil(64) as usize;
|
||||
let all = |one: bool| -> Vec<u64> { (0..words).map(|w| if !one { 0 } else if (w + 1) * 64 <= n as usize { !0 } else { (1u64 << (n as usize % 64)) - 1 }).collect() };
|
||||
let clean: u64 = counts.iter().filter(|(b, _)| *b == all(false) || *b == all(true)).map(|c| c.1).sum();
|
||||
// About 1999 pair channels at 1e-4, of which 12/15 flip some bit:
|
||||
// e^(−0.16) of the shots stay clean.
|
||||
let f = clean as f64 / 2000.0;
|
||||
assert!((f - (-0.16f64).exp()).abs() < 0.05, "{f}");
|
||||
// ⟨Z₀ Z_{n−1}⟩: flipped by an X or Y on the chain between them.
|
||||
let e = frame_expectations(n, &steps, &[vec![(0, 'Z'), (n - 1, 'Z')]], 20_000, 3, 4).unwrap();
|
||||
assert!(e[0].0 > 0.8 && e[0].0 < 1.0, "{e:?}");
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -72,11 +72,14 @@
|
|||
//! essentially free: the same cost at a cap 16 times smaller than the
|
||||
//! unsliced intermediate. Far below the natural size it is not: at `m = 14`
|
||||
//! with a `2³⁰` cap, a plan costs about twelve times its unsliced cost.
|
||||
//! - **The executor is scalar.** It runs about a billion complex multiply-adds
|
||||
//! per second per core, threaded over slices as memory allows, with no SIMD
|
||||
//! kernels or GPU.
|
||||
//! Plans worth `10¹³` operations or more are for planning and verification
|
||||
//! arithmetic here, not for contraction.
|
||||
//! - **The executor is a blocked, threaded loop, without vector kernels or a
|
||||
//! GPU.** Each pairwise contraction is a complex matrix product blocked for
|
||||
//! the cache, its rows shared between threads (slices run in parallel
|
||||
//! instead when there are several). Each output adds its terms in the same
|
||||
//! order however the work is split, so the result is the same bits. On a
|
||||
//! 16-core desktop processor it runs about 25 billion complex multiply-adds a
|
||||
//! second on large products. Plans worth `10¹⁴` operations or more are for
|
||||
//! planning and verification arithmetic here, not for contraction.
|
||||
//! - **Networks are simple graphs.** Every index joins at most two tensors.
|
||||
//! Diagonal-gate hyperedges are not exploited.
|
||||
|
||||
|
|
@ -276,25 +279,73 @@ pub fn contract(a: &Tensor, b: &Tensor, dims: &[usize], open: &[u32]) -> Tensor
|
|||
let ad = permute(a, &ao, dims);
|
||||
let bd = permute(b, &bo, dims);
|
||||
let mut out = vec![C::ZERO; m * n];
|
||||
for i in 0..m {
|
||||
let row = &mut out[i * n..(i + 1) * n];
|
||||
for kk in 0..k {
|
||||
let x = ad[i * k + kk];
|
||||
if x.re == 0.0 && x.im == 0.0 {
|
||||
continue;
|
||||
}
|
||||
let brow = &bd[kk * n..(kk + 1) * n];
|
||||
for (o, y) in row.iter_mut().zip(brow) {
|
||||
o.re += x.re * y.re - x.im * y.im;
|
||||
o.im += x.re * y.im + x.im * y.re;
|
||||
}
|
||||
}
|
||||
}
|
||||
matmul_into(&mut out, &ad, &bd, m, k, n);
|
||||
let mut inds = a_rem;
|
||||
inds.extend(b_rem);
|
||||
Tensor { inds, data: out }
|
||||
}
|
||||
|
||||
/// Columns and inner terms a block of [`matmul_into`] holds: a block of `b`
|
||||
/// (`KB × JB` complex values, 64 KiB) stays in cache while every row of the
|
||||
/// block uses it.
|
||||
const JB: usize = 256;
|
||||
const KB: usize = 16;
|
||||
|
||||
#[cfg(not(target_arch = "wasm32"))]
|
||||
std::thread_local! {
|
||||
/// Set inside callers that already run contractions in parallel.
|
||||
static ONE_THREAD: core::cell::Cell<bool> = const { core::cell::Cell::new(false) };
|
||||
}
|
||||
|
||||
/// `out (m × n) += a (m × k) · b (k × n)`, row-major. Each output element
|
||||
/// adds its terms one at a time in `k` order (a zero factor of `a` skipped),
|
||||
/// so blocking for the cache, and splitting rows between threads, give the
|
||||
/// unblocked loop's bits.
|
||||
fn matmul_into(out: &mut [C], a: &[C], b: &[C], m: usize, k: usize, n: usize) {
|
||||
let rows = |out: &mut [C], i0: usize| {
|
||||
let mrows = out.len() / n.max(1);
|
||||
for j0 in (0..n).step_by(JB) {
|
||||
let j1 = (j0 + JB).min(n);
|
||||
for k0 in (0..k).step_by(KB) {
|
||||
let k1 = (k0 + KB).min(k);
|
||||
for i in 0..mrows {
|
||||
let row = &mut out[i * n + j0..i * n + j1];
|
||||
for kk in k0..k1 {
|
||||
let x = a[(i0 + i) * k + kk];
|
||||
if x.re == 0.0 && x.im == 0.0 {
|
||||
continue;
|
||||
}
|
||||
let brow = &b[kk * n + j0..kk * n + j1];
|
||||
for (o, y) in row.iter_mut().zip(brow) {
|
||||
o.re += x.re * y.re - x.im * y.im;
|
||||
o.im += x.re * y.im + x.im * y.re;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
#[cfg(not(target_arch = "wasm32"))]
|
||||
{
|
||||
let threads = if ONE_THREAD.with(|o| o.get()) { 1 } else { std::thread::available_parallelism().map_or(1, |t| t.get()) };
|
||||
// Worth threading past a few million multiply-adds, with rows to share.
|
||||
let threads = threads.min(m / 4).max(1);
|
||||
if threads > 1 && (m as u64) * (k as u64) * (n as u64) >= 1 << 22 && n > 0 {
|
||||
let per = m.div_ceil(threads);
|
||||
std::thread::scope(|sc| {
|
||||
for (c, chunk) in out.chunks_mut(per * n).enumerate() {
|
||||
let rows = &rows;
|
||||
sc.spawn(move || rows(chunk, c * per));
|
||||
}
|
||||
});
|
||||
return;
|
||||
}
|
||||
}
|
||||
if n > 0 {
|
||||
rows(out, 0);
|
||||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Contraction paths
|
||||
// ---------------------------------------------------------------------------
|
||||
|
|
@ -1288,7 +1339,12 @@ pub fn contract_plan(net: &Network, plan: &Plan) -> Result<Tensor, TnError> {
|
|||
for (c, chunk) in chunks {
|
||||
let run = &run;
|
||||
let base = c * count.div_ceil(threads);
|
||||
let inner_one = threads > 1;
|
||||
sc.spawn(move || {
|
||||
// Slices in parallel: each slice's contractions on its own thread.
|
||||
if inner_one {
|
||||
ONE_THREAD.with(|o| o.set(true));
|
||||
}
|
||||
for (k, slot) in chunk.iter_mut().enumerate() {
|
||||
*slot = Some(run(base + k));
|
||||
}
|
||||
|
|
@ -1320,6 +1376,62 @@ pub fn contract_plan(net: &Network, plan: &Plan) -> Result<Tensor, TnError> {
|
|||
pub(crate) mod tests {
|
||||
use super::*;
|
||||
|
||||
/// The plain i-k-j loop the executor began with.
|
||||
fn matmul_plain(a: &[C], b: &[C], m: usize, k: usize, n: usize) -> Vec<C> {
|
||||
let mut out = vec![C::ZERO; m * n];
|
||||
for i in 0..m {
|
||||
let row = &mut out[i * n..(i + 1) * n];
|
||||
for kk in 0..k {
|
||||
let x = a[i * k + kk];
|
||||
if x.re == 0.0 && x.im == 0.0 {
|
||||
continue;
|
||||
}
|
||||
for (o, y) in row.iter_mut().zip(&b[kk * n..(kk + 1) * n]) {
|
||||
o.re += x.re * y.re - x.im * y.im;
|
||||
o.im += x.re * y.im + x.im * y.re;
|
||||
}
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// Times; run with `--ignored --nocapture`.
|
||||
#[test]
|
||||
#[ignore]
|
||||
fn product_timings() {
|
||||
let mut rng = Rng(4);
|
||||
for (m, k, n) in [(1024usize, 1024usize, 1024usize), (4096, 256, 4096), (64, 4096, 4096)] {
|
||||
let a: Vec<C> = (0..m * k).map(|_| C::new(rng.unit() - 0.5, rng.unit() - 0.5)).collect();
|
||||
let b: Vec<C> = (0..k * n).map(|_| C::new(rng.unit() - 0.5, rng.unit() - 0.5)).collect();
|
||||
let t0 = std::time::Instant::now();
|
||||
let want = matmul_plain(&a, &b, m, k, n);
|
||||
let tp = t0.elapsed().as_secs_f64();
|
||||
let mut got = vec![C::ZERO; m * n];
|
||||
let t1 = std::time::Instant::now();
|
||||
matmul_into(&mut got, &a, &b, m, k, n);
|
||||
let tb = t1.elapsed().as_secs_f64();
|
||||
let macs = (m * k * n) as f64;
|
||||
println!("{m}x{k}x{n}: plain {tp:.3} s ({:.2} G/s), blocked+threaded {tb:.3} s ({:.2} G/s); same {}", macs / tp / 1e9, macs / tb / 1e9, got == want);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn blocked_threaded_products_are_the_plain_loops_bits() {
|
||||
let mut rng = Rng(3);
|
||||
for (m, k, n) in [(1usize, 1usize, 1usize), (3, 17, 5), (64, 64, 64), (300, 40, 513), (1024, 33, 700), (9, 1000, 2)] {
|
||||
let mut a: Vec<C> = (0..m * k).map(|_| C::new(rng.unit() - 0.5, rng.unit() - 0.5)).collect();
|
||||
// Some exact zeros, which both skip.
|
||||
for x in a.iter_mut().step_by(7) {
|
||||
*x = C::ZERO;
|
||||
}
|
||||
let b: Vec<C> = (0..k * n).map(|_| C::new(rng.unit() - 0.5, rng.unit() - 0.5)).collect();
|
||||
let want = matmul_plain(&a, &b, m, k, n);
|
||||
let mut got = vec![C::ZERO; m * n];
|
||||
matmul_into(&mut got, &a, &b, m, k, n);
|
||||
assert!(got.iter().zip(&want).all(|(x, y)| x.re.to_bits() == y.re.to_bits() && x.im.to_bits() == y.im.to_bits()), "{m}×{k}×{n}");
|
||||
}
|
||||
}
|
||||
|
||||
/// A random unitary on `k` qubits from the QR of a random matrix
|
||||
/// (Gram–Schmidt on columns): enough to exercise the contractions.
|
||||
#[allow(clippy::needless_range_loop)]
|
||||
|
|
|
|||
|
|
@ -620,9 +620,22 @@ pub fn to_qasm(c: &Circuit) -> String {
|
|||
(BaseGate::Y, [c0]) => format!("cy q[{c0}],q[{t}];"),
|
||||
(BaseGate::Z, [c0]) => format!("cz q[{c0}],q[{t}];"),
|
||||
(BaseGate::P, [c0]) => format!("cp(2*pi/{}) q[{c0}],q[{t}];", 1u64 << g.param),
|
||||
(BaseGate::H, [c0]) => format!("ch q[{c0}],q[{t}];"),
|
||||
(BaseGate::X, [a, b]) => format!("ccx q[{a}],q[{b}],q[{t}];"),
|
||||
(BaseGate::Z, [a, b]) => format!("ccz q[{a}],q[{b}],q[{t}];"),
|
||||
(base, ctrls) => format!("// unmapped: {base:?} controls={ctrls:?} target={t}"),
|
||||
// A controlled identity is the identity.
|
||||
(BaseGate::I, _) => continue,
|
||||
// Every other controlled gate in OpenQASM 3's modifier form, which
|
||||
// `from_qasm` reads back: nothing is left out.
|
||||
(base, ctrls) => {
|
||||
let name = match base {
|
||||
BaseGate::P => format!("p(2*pi/{})", 1u64 << g.param),
|
||||
b => format!("{b:?}").to_lowercase(),
|
||||
};
|
||||
let m = if ctrls.len() == 1 { "ctrl".to_string() } else { format!("ctrl({})", ctrls.len()) };
|
||||
let qs: Vec<String> = ctrls.iter().chain([&t]).map(|q| format!("q[{q}]")).collect();
|
||||
format!("{m} @ {name} {};", qs.join(","))
|
||||
}
|
||||
};
|
||||
s.push_str(&line);
|
||||
s.push('\n');
|
||||
|
|
@ -689,6 +702,11 @@ pub fn from_qasm(src: &str) -> Result<Circuit, String> {
|
|||
if line.contains("measure") {
|
||||
continue;
|
||||
}
|
||||
// Gate modifiers: `ctrl @`, `ctrl(n) @`, `negctrl @`, `negctrl(n) @`.
|
||||
if let Some((polarity, call)) = modifiers(&line)? {
|
||||
push_modified(&mut c, &polarity, &call)?;
|
||||
continue;
|
||||
}
|
||||
let head: String = line.chars().take_while(|c| c.is_ascii_alphanumeric() || *c == '_').collect();
|
||||
let head = head.as_str();
|
||||
match head {
|
||||
|
|
@ -805,6 +823,92 @@ pub fn from_qasm(src: &str) -> Result<Circuit, String> {
|
|||
Ok(c)
|
||||
}
|
||||
|
||||
/// The polarities of a line's gate modifiers, in order (`true`: a control on
|
||||
/// |1⟩, `false`: on |0⟩), and the gate call after them; `None` when the line
|
||||
/// has none.
|
||||
fn modifiers(line: &str) -> Result<Option<(Vec<bool>, String)>, String> {
|
||||
let mut rest = line.trim_start();
|
||||
let mut polarity = Vec::new();
|
||||
loop {
|
||||
let (on, after) = if let Some(a) = rest.strip_prefix("negctrl") {
|
||||
(false, a)
|
||||
} else if let Some(a) = rest.strip_prefix("ctrl") {
|
||||
(true, a)
|
||||
} else {
|
||||
break;
|
||||
};
|
||||
let after = after.trim_start();
|
||||
let (count, after) = match after.strip_prefix('(') {
|
||||
Some(a) => {
|
||||
let end = a.find(')').ok_or("unterminated modifier argument")?;
|
||||
let n: usize = a[..end].trim().parse().map_err(|_| format!("the modifier's count `{}` is not a whole number", a[..end].trim()))?;
|
||||
(n, a[end + 1..].trim_start())
|
||||
}
|
||||
None => (1, after),
|
||||
};
|
||||
let after = after.strip_prefix('@').ok_or("expected `@` after a gate modifier")?;
|
||||
if count == 0 {
|
||||
return Err("a control modifier adds at least one control".into());
|
||||
}
|
||||
polarity.extend(std::iter::repeat_n(on, count));
|
||||
rest = after.trim_start();
|
||||
}
|
||||
Ok((!polarity.is_empty()).then(|| (polarity, rest.to_string())))
|
||||
}
|
||||
|
||||
/// A gate under control modifiers: the added controls are the first qubits
|
||||
/// named, in the modifiers' order, and a control on |0⟩ is a control between
|
||||
/// two `x`. Only gates read here exactly, global phase and all, may be
|
||||
/// controlled, since a control makes a global phase observable: `rz`, `rx`,
|
||||
/// `ry`, `sx` and `u`, exact here up to their global phase, are refused.
|
||||
fn push_modified(c: &mut Circuit, polarity: &[bool], call: &str) -> Result<(), String> {
|
||||
let head: String = call.chars().take_while(|c| c.is_ascii_alphanumeric() || *c == '_').collect();
|
||||
// (the base gate, the controls the name already has, whether it takes an angle)
|
||||
let (base, inner, angle) = match head.as_str() {
|
||||
"x" => (BaseGate::X, 0, false),
|
||||
"y" => (BaseGate::Y, 0, false),
|
||||
"z" => (BaseGate::Z, 0, false),
|
||||
"h" => (BaseGate::H, 0, false),
|
||||
"s" => (BaseGate::S, 0, false),
|
||||
"sdg" => (BaseGate::Sdg, 0, false),
|
||||
"t" => (BaseGate::T, 0, false),
|
||||
"tdg" => (BaseGate::Tdg, 0, false),
|
||||
"p" | "u1" | "phase" => (BaseGate::P, 0, true),
|
||||
"cx" | "cnot" => (BaseGate::X, 1, false),
|
||||
"cy" => (BaseGate::Y, 1, false),
|
||||
"cz" => (BaseGate::Z, 1, false),
|
||||
"ch" => (BaseGate::H, 1, false),
|
||||
"cp" | "cu1" | "cphase" => (BaseGate::P, 1, true),
|
||||
"ccx" | "toffoli" => (BaseGate::X, 2, false),
|
||||
"ccz" => (BaseGate::Z, 2, false),
|
||||
other => {
|
||||
return Err(format!(
|
||||
"a control modifier on '{other}' is not supported (supported under one: x y z h s sdg t tdg p cx cy cz ch cp ccx ccz)"
|
||||
));
|
||||
}
|
||||
};
|
||||
let qs = parse_qubits(call);
|
||||
let want = polarity.len() + inner + 1;
|
||||
if qs.len() != want {
|
||||
return Err(format!("'{head}' under {} added control(s) takes {want} qubits, given {}", polarity.len(), qs.len()));
|
||||
}
|
||||
let (controls, target) = (qs[..want - 1].to_vec(), qs[want - 1]);
|
||||
distinct_qubits(&controls, target)?;
|
||||
let flip: Vec<u8> = controls.iter().zip(polarity).filter(|(_, on)| !**on).map(|(q, _)| *q).collect();
|
||||
for &q in &flip {
|
||||
c.x(q);
|
||||
}
|
||||
if angle {
|
||||
push_turn(c, &controls, target, parse_angle(&arg_n(call, 0)?)?)?;
|
||||
} else {
|
||||
c.ops.push(Gate { base, controls, target, param: 0 });
|
||||
}
|
||||
for &q in &flip {
|
||||
c.x(q);
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Refuse a gate that names one qubit twice, a control repeated or equal to
|
||||
/// the target, as an invalid circuit in every gateset.
|
||||
fn distinct_qubits(controls: &[u8], target: u8) -> Result<(), String> {
|
||||
|
|
@ -1220,7 +1324,9 @@ mod tests {
|
|||
|
||||
#[test]
|
||||
fn qasm_round_trip() {
|
||||
for c in [bell(), ghz(3), Circuit::qft(4), grover(3, 5)] {
|
||||
// Grover's search past three qubits uses a Z with three or more
|
||||
// controls, which the export once wrote as a comment and lost.
|
||||
for c in [bell(), ghz(3), Circuit::qft(4), grover(3, 5), grover(5, 0b10110), grover(7, 77)] {
|
||||
let qasm = to_qasm(&c);
|
||||
let back = from_qasm(&qasm).expect("reparse");
|
||||
// same reconstructed state ⇒ same circuit semantics
|
||||
|
|
@ -1232,6 +1338,46 @@ mod tests {
|
|||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn every_controlled_gate_is_exported_and_read_back_as_it_was() {
|
||||
use BaseGate::*;
|
||||
for base in [X, Y, Z, H, S, Sdg, T, Tdg, P] {
|
||||
for nc in 0..=5u8 {
|
||||
let mut c = Circuit::new(7);
|
||||
c.ops.push(Gate { base, controls: (0..nc).collect(), target: 6, param: if base == P { 3 } else { 0 } });
|
||||
let qasm = to_qasm(&c);
|
||||
assert!(!qasm.contains("//"), "{qasm}");
|
||||
assert_eq!(from_qasm(&qasm).unwrap().ops, c.ops, "{qasm}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn control_modifiers_are_read_with_their_polarity() {
|
||||
let c = from_qasm("qubit[4] q;\nnegctrl @ ctrl(2) @ z q[0], q[1], q[2], q[3];\nctrl @ cp(pi/4) q[1], q[2], q[0];").unwrap();
|
||||
let x0 = Gate { base: BaseGate::X, controls: vec![], target: 0, param: 0 };
|
||||
assert_eq!(
|
||||
c.ops,
|
||||
vec![
|
||||
x0.clone(),
|
||||
Gate { base: BaseGate::Z, controls: vec![0, 1, 2], target: 3, param: 0 },
|
||||
x0,
|
||||
Gate { base: BaseGate::P, controls: vec![1, 2], target: 0, param: 3 },
|
||||
]
|
||||
);
|
||||
// A global phase a control would expose, a count of nothing, and a
|
||||
// wrong number of qubits are refused.
|
||||
for (src, needle) in [
|
||||
("qubit[2] q;\nctrl @ rz(pi/2) q[0], q[1];", "not supported"),
|
||||
("qubit[2] q;\nctrl(0) @ x q[0], q[1];", "at least one"),
|
||||
("qubit[3] q;\nctrl @ x q[0], q[1], q[2];", "takes 2 qubits"),
|
||||
("qubit[2] q;\nctrl @ x q[1], q[1];", "equals target"),
|
||||
] {
|
||||
let e = from_qasm(src).unwrap_err();
|
||||
assert!(e.contains(needle), "{src}: {e}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn qasm_imports_foreign_circuit() {
|
||||
let src = "OPENQASM 3.0;\ninclude \"stdgates.inc\";\nqubit[2] q;\nh q[0];\ncx q[0], q[1];\n";
|
||||
|
|
|
|||
|
|
@ -38,28 +38,28 @@ pub use trig::sin_cos;
|
|||
// Natural logarithm
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LG1: f64 = 6.666_666_666_666_735_130e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LG2: f64 = 3.999_999_999_940_941_908e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LG3: f64 = 2.857_142_874_366_239_149e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LG4: f64 = 2.222_219_843_214_978_396e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LG5: f64 = 1.818_357_216_161_805_012e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
const LG6: f64 = 1.531_383_769_920_937_332e-1;
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
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", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
pub fn ln(x: f64) -> f64 {
|
||||
let mut bits = x.to_bits();
|
||||
let mut hx = (bits >> 32) as u32;
|
||||
|
|
@ -359,7 +359,7 @@ mod tests {
|
|||
})
|
||||
}
|
||||
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy"))]
|
||||
#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel", feature = "quantum_resource", feature = "quantum_frame", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_arch", feature = "quantum_qio", feature = "quantum_dmrg", feature = "quantum_noisy", feature = "quantum_tableau"))]
|
||||
#[test]
|
||||
fn ln_is_within_an_ulp_of_the_platform() {
|
||||
for x in inputs().filter(|&x| x > 0.0) {
|
||||
|
|
|
|||
2
wai/wai-rs/Cargo.lock
generated
2
wai/wai-rs/Cargo.lock
generated
|
|
@ -4310,7 +4310,7 @@ version = "0.1.0"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-web/Cargo.lock
generated
2
wai/wai-web/Cargo.lock
generated
|
|
@ -1241,7 +1241,7 @@ version = "0.1.0"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.6.3"
|
||||
version = "0.6.4"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
Loading…
Reference in a new issue