open-standards/wai/quantum-conformance/tools/referee.py
Transaction Science e1f3b6faa0 wai-quantum circuit2: the number range decided on the literal's digits in both readers
A third review upheld six findings.

- std's f64 parse stops adding exponent digits once the exponent reaches
  0x10000, so a literal whose exponent offsets about 655,000 zeros was judged
  wrong both ways: exactly 1 refused, 10^309 accepted. The Rust reader now
  decides the range rule on the digits, as the referee does: the significant
  digits and the decimal exponent, with saturating exponent arithmetic, against
  the 309 digits of 2^1024 − 2^970. The limit is one constant, RANGE_LIMIT, and
  a test derives it in decimal from 2^1024 and 2^970 and checks it against
  std's parse where that parse is exact.
- The referee passed a zero-padded exponent of over 4300 digits to int(),
  which Python refuses; it now converts only the significant exponent digits.
- quantum-sim §8, the CHANGELOG and SPEC §5 now list both directions of the
  range change: numbers at or above the limit that earlier builds let through
  are refused, numbers below it that some builds refused are accepted. §2 says
  the rule is decided on the digits, at any length.
- The corpus README and the wasip2 job comment were stale.

The corpus holds both long literals and the padded exponent: an accept file
(exactly 1 and exactly 10, written long) and a reject file (10^309).

Co-Authored-By: Claude Opus 5.5 <[email protected]>
2026-10-09 18:50:35 -04:00

1037 lines
42 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/usr/bin/env python3
"""Referee for quantum-conformance/: both capabilities of the quantum-sim extension
(extensions/quantum-sim.md), `wai.quantum.circuit` (WQC1) and `wai.quantum.circuit2`
(WQC2), recomputed with integers only.
Written from the extension's rules. It shares no code with wai-rs or wai-quantum,
which generate the corpus: the phase tables, the container writer and reader, the
fixed-point simulator, the shot sampler and BLAKE3 are all written here. Standard
library only, and no floating point anywhere.
It checks, in this order, and reports every failure rather than stopping at one:
1. BLAKE3 against the published test vectors.
2. The phase tables. circuit/1's is derived by its recurrence (§3.1) and circuit2's
by an interval-certified procedure in exact integers (§3.2): each step keeps a
true enclosure of 2^62·cos(2π/2^k) and 2^62·sin(2π/2^k), and both ends of every
enclosure round to the same integer far from a tie, so that integer is the
nearest one. Each derived table must equal its literal table (restated below) and
its BLAKE3 digest, and the extension's Appendix B must state those same tables and
digests.
3. The vector trees: `./` (circuit/1) and `v2/` (circuit2). Every file is rebuilt byte
for byte from the circuit definitions below and read back by the reader here, and
every manifest entry must equal the entry computed here: n_qubits, n_ops, the
circuit, statevector and histogram hashes, the measurement, and in `v2/` the
capability. The four original circuit/1 vectors must also keep their frozen hashes.
4. The reject tree: `reject.json` must list exactly the cases below, the directory
must hold exactly the listed files, and the reader here must refuse each one with
its listed class and by the rule its name states.
5. The accept tree: `accept.json` must list exactly the cases below, the directory
must hold exactly the listed files, and the reader here must accept each one and
read it as listed: its capability, n_qubits, op count, circuit hash and
measurement.
6. Mutation completeness: under each table, a change of ±1 to either component of any
P(1..32), with H following re P(3), must change dyadic_ladder's statevector hash.
It was written from the same reading of the text as the Rust reference, so a misreading
the two share would pass both.
python3 -I referee.py [quantum-conformance dir [quantum-sim.md]]
The extension defaults to ../extensions/quantum-sim.md beside the corpus directory.
Exit status 0 only when every check passes.
"""
import json
import os
import struct
import sys
from math import isqrt
# ---------------------------------------------------------------- BLAKE3 (hash mode)
_IV = [0x6A09E667, 0xBB67AE85, 0x3C6EF372, 0xA54FF53A, 0x510E527F, 0x9B05688C, 0x1F83D9AB, 0x5BE0CD19]
_PERM = [2, 6, 3, 10, 7, 0, 4, 13, 1, 11, 12, 5, 9, 14, 15, 8]
_CHUNK_START, _CHUNK_END, _PARENT, _ROOT = 1, 2, 4, 8
_M32 = 0xFFFFFFFF
_CHUNK = 1024
_BLOCK = 64
def _rotr(x, n):
return ((x >> n) | (x << (32 - n))) & _M32
def _g(s, a, b, c, d, x, y):
s[a] = (s[a] + s[b] + x) & _M32
s[d] = _rotr(s[d] ^ s[a], 16)
s[c] = (s[c] + s[d]) & _M32
s[b] = _rotr(s[b] ^ s[c], 12)
s[a] = (s[a] + s[b] + y) & _M32
s[d] = _rotr(s[d] ^ s[a], 8)
s[c] = (s[c] + s[d]) & _M32
s[b] = _rotr(s[b] ^ s[c], 7)
def _compress(cv, block, counter, block_len, flags):
m = list(struct.unpack("<16I", block))
s = list(cv) + _IV[:4] + [counter & _M32, (counter >> 32) & _M32, block_len, flags]
for r in range(7):
_g(s, 0, 4, 8, 12, m[0], m[1])
_g(s, 1, 5, 9, 13, m[2], m[3])
_g(s, 2, 6, 10, 14, m[4], m[5])
_g(s, 3, 7, 11, 15, m[6], m[7])
_g(s, 0, 5, 10, 15, m[8], m[9])
_g(s, 1, 6, 11, 12, m[10], m[11])
_g(s, 2, 7, 8, 13, m[12], m[13])
_g(s, 3, 4, 9, 14, m[14], m[15])
if r < 6:
m = [m[i] for i in _PERM]
return [s[i] ^ s[i + 8] for i in range(8)] + [s[i + 8] ^ cv[i] for i in range(8)]
def _chunk_output(chunk, counter):
"""The last block of a chunk, ready to compress: (cv, block, counter, len, flags)."""
cv = _IV
n_blocks = max(1, (len(chunk) + _BLOCK - 1) // _BLOCK)
for b in range(n_blocks):
block = chunk[b * _BLOCK:(b + 1) * _BLOCK]
flags = _CHUNK_START if b == 0 else 0
if b == n_blocks - 1:
return (cv, block.ljust(_BLOCK, b"\0"), counter, len(block), flags | _CHUNK_END)
cv = _compress(cv, block, counter, _BLOCK, flags)[:8]
def _cv(output):
cv, block, counter, block_len, flags = output
return _compress(cv, block, counter, block_len, flags)[:8]
def _parent(left_cv, right_cv):
return (_IV, struct.pack("<8I", *left_cv) + struct.pack("<8I", *right_cv), 0, _BLOCK, _PARENT)
def blake3(data):
"""The 32-byte BLAKE3 hash of `data`."""
chunks = [data[i:i + _CHUNK] for i in range(0, len(data), _CHUNK)] or [b""]
stack = []
for i, chunk in enumerate(chunks[:-1]):
cv = _cv(_chunk_output(chunk, i))
total = i + 1
while total & 1 == 0:
cv = _cv(_parent(stack.pop(), cv))
total >>= 1
stack.append(cv)
node = _chunk_output(chunks[-1], len(chunks) - 1)
while stack:
node = _parent(stack.pop(), _cv(node))
cv, block, counter, block_len, flags = node
return struct.pack("<8I", *_compress(cv, block, counter, block_len, flags | _ROOT)[:8])
# The published BLAKE3 test vectors (input byte i is i mod 251): input length, hash.
BLAKE3_VECTORS = [
(0, "af1349b9f5f9a1a6a0404dea36dcc9499bcb25c9adc112b7cc9a93cae41f3262"),
(1, "2d3adedff11b61f14c886e35afa036736dcd87a74d27b5c1510225d0f592e213"),
(63, "e9bc37a594daad83be9470df7f7b3798297c3d834ce80ba85d6e207627b7db7b"),
(64, "4eed7141ea4a5cd4b788606bd23f46e212af9cacebacdc7d1f4c6dc7f2511b98"),
(65, "de1e5fa0be70df6d2be8fffd0e99ceaa8eb6e8c93a63f2d8d1c30ecb6b263dee"),
(1023, "10108970eeda3eb932baac1428c7a2163b0e924c9a9e25b35bba72b28f70bd11"),
(1024, "42214739f095a406f3fc83deb889744ac00df831c10daa55189b5d121c855af7"),
(1025, "d00278ae47eb27b34faecf67b4fe263f82d5412916c1ffd97c8cb7fb814b8444"),
(2048, "e776b6028c7cd22a4d0ba182a8bf62205d2ef576467e838ed6f2529b85fba24a"),
(2049, "5f4d72f40d7a5f82b15ca2b2e44b1de3c2ef86c426c95c1af0b6879522563030"),
(3072, "b98cb0ff3623be03326b373de6b9095218513e64f1ee2edd2525c7ad1e5cffd2"),
(3073, "7124b49501012f81cc7f11ca069ec9226cecb8a2c850cfe644e327d22d3e1cd3"),
(4096, "015094013f57a5277b59d8475c0501042c0b642e531b0a1c8f58d2163229e969"),
(4097, "9b4052b38f1c5fc8b1f9ff7ac7b27cd242487b3d890d15c96a1c25b8aa0fb995"),
(8193, "bab6c09cb8ce8cf459261398d2e7aef35700bf488116ceb94a36d0f5f1b7bc3b"),
(16384, "f875d6646de28985646f34ee13be9a576fd515f76b5b0a26bb324735041ddde4"),
(31744, "62b6960e1a44bcc1eb1a611a8d6235b6b4b78f32e7abc4fb4c6cdcce94895c47"),
]
# ---------------------------------------------------------------- the phase tables
FRAC = 30
ONE = 1 << FRAC
HALF = 1 << (FRAC - 1)
DYADIC_MAX = 32
# Appendix B.1: circuit/1, (re, im) for k = 1..32. P(19..32) apply no phase at all.
PHASE_V1 = [None] + [
(-1073741824, 0), (0, 1073741824), (759250124, 759250124), (992008094, 410903207),
(1053110175, 209476638), (1068571463, 105245102), (1072448454, 52686012),
(1073418432, 26350953), (1073660972, 13176485), (1073721610, 6588405),
(1073736770, 3294284), (1073740560, 1647223), (1073741507, 823774),
(1073741744, 411887), (1073741803, 207243), (1073741818, 103621),
(1073741822, 56755), (1073741823, 32768),
] + [(1073741823, 0)] * 14
PHASE_V1_DIGEST = "3fa9805ff96feb819a1cf371fedcd3b5228f93eddf933d0edcf401cbf8c28466"
# Appendix B.2: circuit2, the integers nearest 2^30·cos(2π/2^k) and 2^30·sin(2π/2^k).
PHASE_V2 = [None] + [
(-1073741824, 0), (0, 1073741824), (759250125, 759250125), (992008094, 410903207),
(1053110176, 209476638), (1068571464, 105245103), (1072448455, 52686014),
(1073418433, 26350943), (1073660973, 13176464), (1073721611, 6588356),
(1073736771, 3294193), (1073740561, 1647099), (1073741508, 823550),
(1073741745, 411775), (1073741804, 205887), (1073741819, 102944),
(1073741823, 51472),
] + [(ONE, im) for im in (25736, 12868, 6434, 3217, 1608, 804, 402, 201, 101, 50, 25, 13, 6, 3, 2)]
PHASE_V2_DIGEST = "d34cb793715348e7c92db027c0bc3dddb6008f5730470bafaed978ea7fb23510"
def derive_v1():
"""circuit/1's table by its recurrence: with c = re P(k-1),
re P(k) = isqrt(floor((2^30 + c) / 2) · 2^30), im P(k) = isqrt(floor((2^30 - c) / 2) · 2^30)."""
t = [None, (-ONE, 0), (0, ONE)]
for k in range(3, DYADIC_MAX + 1):
c = t[k - 1][0]
t.append((isqrt(((ONE + c) // 2) << FRAC), isqrt(((ONE - c) // 2) << FRAC)))
return t
def derive_v2():
"""circuit2's table by enclosures at scale 2^62, and the reasons it is not certified.
cos(θ/2) = √((1 + cos θ)/2) and sin(θ/2) = sin θ / (2 cos(θ/2)), each bound rounded
outward, starting from (cos, sin)(π/2) = (0, 1), which is exact."""
s = 1 << 62
def ceil_isqrt(x):
r = isqrt(x)
return r if r * r == x else r + 1
def ceil_div(a, b):
return -((-a) // b)
def nearest(x):
# From 2^-62 to 2^-30, half up.
return (x + (1 << 31)) >> 32
problems = []
c_lo, c_hi, s_lo, s_hi = 0, 0, s, s
t = [None, (-ONE, 0), (0, ONE)]
for k in range(3, DYADIC_MAX + 1):
c_lo2 = isqrt(((s + c_lo) * s) >> 1)
c_hi2 = ceil_isqrt(ceil_div((s + c_hi) * s, 2))
s_lo2 = (s_lo * s) // (2 * c_hi2)
s_hi2 = ceil_div(s_hi * s, 2 * c_lo2)
if nearest(c_lo2) != nearest(c_hi2) or nearest(s_lo2) != nearest(s_hi2):
problems.append("k=%d: the enclosure's ends round apart" % k)
for x in (c_lo2, c_hi2, s_lo2, s_hi2):
if abs((x % (1 << 32)) - (1 << 31)) <= (1 << 20):
problems.append("k=%d: an end lies within 2^20 of a tie" % k)
t.append((nearest(c_lo2), nearest(s_lo2)))
c_lo, c_hi, s_lo, s_hi = c_lo2, c_hi2, s_lo2, s_hi2
return t, problems
def table_digest(t):
"""BLAKE3 over i64le(re) ‖ i64le(im) for k = 1..32, no domain prefix."""
return blake3(b"".join(struct.pack("<qq", *t[k]) for k in range(1, DYADIC_MAX + 1))).hex()
# ---------------------------------------------------------------- the capabilities
class Gateset:
def __init__(self, label, magic, capability, ext, numeric, name, table):
self.label, self.magic, self.capability = label, magic, capability
self.ext, self.numeric, self.name, self.table = ext, numeric, name, table
def contract(self, n_qubits):
"""The canonical contract bytes (§2.1)."""
return ('{"ext":"%s","n_qubits":%d,"numeric":"%s","frac":%d,"gateset":"%s"}'
% (self.ext, n_qubits, self.numeric, FRAC, self.name)).encode()
V1 = Gateset("circuit/1", b"WQC1", "wai.quantum.circuit", "wai.quantum.circuit/1",
"wai.det.fixed64", "cliffordT+dyadicP", PHASE_V1)
V2 = Gateset("circuit2", b"WQC2", "wai.quantum.circuit2", "wai.quantum.circuit2/1",
"wai.det.amp64", "cliffordT+dyadicP.rn", PHASE_V2)
BY_MAGIC = {V1.magic: V1, V2.magic: V2}
SECT_CONTRACT, SECT_OPLOG, SECT_MEASURE = 0x01, 0x02, 0x03
I, X, Y, Z, H, S, SDG, T, TDG, P = range(10)
def gate(op, target, controls=(), param=0):
"""One op as (opcode, controls, target, param)."""
return (op, tuple(controls), target, param)
# ---------------------------------------------------------------- writer (§2)
def oplog_bytes(ops):
out = bytearray(struct.pack("<I", len(ops)))
for op, controls, target, param in ops:
out += bytes([op, len(controls), target]) + struct.pack("<H", param) + bytes(controls)
return bytes(out)
def measure_bytes(seed, shots):
return ('{"seed":%d,"shots":%d,"basis":"computational"}' % (seed, shots)).encode()
def write_wqc(gs, n_qubits, ops, measure):
sections = [(SECT_CONTRACT, gs.contract(n_qubits)), (SECT_OPLOG, oplog_bytes(ops))]
if measure is not None:
sections.append((SECT_MEASURE, measure_bytes(*measure)))
out = bytearray(gs.magic) + struct.pack("<H", len(sections))
offset = 0
for kind, data in sections:
out += bytes([kind]) + struct.pack("<II", offset, len(data))
offset += len(data)
for _, data in sections:
out += data
return bytes(out)
# ---------------------------------------------------------------- reader (§2, in its check order)
class Refused(Exception):
"""A container the reader refuses: its class (malformed, contract, invalid) and the rule."""
def __init__(self, kind, rule):
super().__init__("%s: %s" % (kind, rule))
self.kind, self.rule = kind, rule
class NotIJson(Exception):
"""A section that is not I-JSON (§2), and the rule it breaks: utf-8, syntax, depth,
duplicate member, surrogate, noncharacter or number range."""
def __init__(self, rule):
super().__init__(rule)
self.rule = rule
MAX_DEPTH = 127 # nested arrays and objects a section may hold, the outermost first (§2)
# 2^1024 − 2^970, the least magnitude that rounds to an infinite double (§2), as digits.
RANGE_LIMIT = str((1 << 1024) - (1 << 970))
def out_of_range(literal):
"""Whether a JSON number's magnitude is at least 2^1024 − 2^970, decided on its
digits: the significant digits D and the exponent K with value 0.D × 10^K."""
sign_free = literal.lstrip("-")
mantissa, _, exponent = sign_free.replace("E", "e").partition("e")
whole, _, frac = mantissa.partition(".")
digits = whole + frac
significant = digits.lstrip("0").rstrip("0")
if not significant:
return False
e = exponent.lstrip("+-").lstrip("0") # never the raw text: int() limits its length
if len(e) > 18: # a section holds at most 2^32 digits, so the sign alone decides
return not exponent.startswith("-")
e = -int(e or "0") if exponent.startswith("-") else int(e or "0")
k = len(whole) - (len(digits) - len(digits.lstrip("0"))) + e
if k != len(RANGE_LIMIT):
return k > len(RANGE_LIMIT)
n = max(len(significant), len(RANGE_LIMIT))
return significant.ljust(n, "0") >= RANGE_LIMIT.ljust(n, "0")
def strict_json(raw):
"""`raw` read as I-JSON (§2), or NotIJson naming the rule it breaks. An integer reads
as a Python int, so `-0` and any number with a fraction or exponent read as
("number", literal): never an unsigned integer, and never a float."""
def number(text):
if out_of_range(text):
raise NotIJson("number range")
return ("number", text)
def integer(text):
if out_of_range(text):
raise NotIJson("number range")
return ("number", text) if text == "-0" else int(text)
def no_constant(name):
raise NotIJson("syntax")
def unique(pairs):
names = [k for k, _ in pairs]
if len(set(names)) != len(names):
raise NotIJson("duplicate member")
return dict(pairs)
def text(s):
for ch in s:
o = ord(ch)
if 0xD800 <= o <= 0xDFFF:
raise NotIJson("surrogate")
if 0xFDD0 <= o <= 0xFDEF or o & 0xFFFE == 0xFFFE:
raise NotIJson("noncharacter")
def walk(v, depth):
if isinstance(v, (dict, list)):
if depth + 1 > MAX_DEPTH:
raise NotIJson("depth")
items = v.items() if isinstance(v, dict) else enumerate(v)
for k, x in items:
if isinstance(k, str):
text(k)
walk(x, depth + 1)
elif isinstance(v, str):
text(v)
try:
decoded = raw.decode("utf-8")
except UnicodeDecodeError:
raise NotIJson("utf-8")
try:
v = json.loads(decoded, parse_constant=no_constant, parse_float=number,
parse_int=integer, object_pairs_hook=unique)
except RecursionError:
raise NotIJson("depth")
except ValueError: # json.JSONDecodeError among them
raise NotIJson("syntax")
walk(v, 0)
return v
def unsigned(v, below):
"""Whether `v` is a JSON unsigned integer less than `below` (not a float, not a boolean)."""
return type(v) is int and 0 <= v < below
def read_contract(gs, raw):
"""The contract's n_qubits, under the capability the magic names."""
try:
cv = strict_json(raw)
except NotIJson as e:
raise Refused("contract", "not JSON: " + e.rule)
if not isinstance(cv, dict):
raise Refused("contract", "not an object")
def n_qubits():
n = cv.get("n_qubits")
if not unsigned(n, 256):
raise Refused("malformed", "n_qubits")
return n
if gs is V1:
# The four values, in any formatting, with any other members.
for key, value in (("ext", gs.ext), ("numeric", gs.numeric)):
if cv.get(key) != value:
raise Refused("contract", key)
if not (unsigned(cv.get("frac"), 1 << 64) and cv["frac"] == FRAC):
raise Refused("contract", "frac")
if cv.get("gateset") != gs.name:
raise Refused("contract", "gateset")
return n_qubits()
n = n_qubits()
if raw != gs.contract(n):
raise Refused("contract", "canonical contract")
return n
def read_oplog(raw):
"""The ops and the bytes they took; bytes after the last op are the caller's to judge."""
if len(raw) < 4:
raise Refused("malformed", "op-log header")
n_ops = struct.unpack_from("<I", raw, 0)[0]
ops, pos = [], 4
for _ in range(n_ops):
if len(raw) - pos < 5:
raise Refused("malformed", "op header")
op, n_ctrl, target = raw[pos], raw[pos + 1], raw[pos + 2]
if op > P:
raise Refused("malformed", "opcode")
param = struct.unpack_from("<H", raw, pos + 3)[0]
pos += 5
if len(raw) - pos < n_ctrl:
raise Refused("malformed", "op controls")
ops.append((op, tuple(raw[pos:pos + n_ctrl]), target, param))
pos += n_ctrl
return ops, pos
def check_valid(gs, n_qubits, ops):
"""§2.2: indices in range, no control on its target, 1 ≤ k ≤ 32 for P; circuit2 also
requires distinct controls and param 0 on every gate other than P."""
for op, controls, target, param in ops:
if target >= n_qubits:
raise Refused("invalid", "target")
if len(controls) > 255:
raise Refused("invalid", "control count")
for c in controls:
if c >= n_qubits:
raise Refused("invalid", "control")
if c == target:
raise Refused("invalid", "control is target")
if gs is V2 and len(set(controls)) != len(controls):
raise Refused("invalid", "repeated control")
if op == P:
if not 1 <= param <= DYADIC_MAX:
raise Refused("invalid", "phase index")
elif gs is V2 and param != 0:
raise Refused("invalid", "param")
def read_measure(gs, raw):
"""circuit/1: JSON, a missing or non-u64 seed or shots reads as 0. circuit2: u64 seed
and shots, and exactly the canonical bytes."""
try:
mv = strict_json(raw)
except NotIJson as e:
raise Refused("malformed", "measure JSON: " + e.rule)
if gs is V1:
def field(key):
v = mv.get(key) if isinstance(mv, dict) else None
return v if unsigned(v, 1 << 64) else 0
return (field("seed"), field("shots"))
for key in ("seed", "shots"):
if not (isinstance(mv, dict) and unsigned(mv.get(key), 1 << 64)):
raise Refused("malformed", "measure " + key)
m = (mv["seed"], mv["shots"])
if raw != measure_bytes(*m):
raise Refused("malformed", "canonical measure")
return m
def read_wqc(data):
"""(gateset, n_qubits, ops, measure or None), or Refused."""
if len(data) < 6 or data[:4] not in BY_MAGIC:
raise Refused("malformed", "magic")
gs = BY_MAGIC[data[:4]]
n_sections = struct.unpack_from("<H", data, 4)[0]
blob = 6 + 9 * n_sections
if blob > len(data):
raise Refused("malformed", "section table")
sections = {}
for i in range(n_sections):
p = 6 + 9 * i
kind = data[p]
offset, length = struct.unpack_from("<II", data, p + 1)
start = blob + offset
end = start + length
if end > len(data):
raise Refused("malformed", "section bounds")
if kind not in (SECT_CONTRACT, SECT_OPLOG, SECT_MEASURE):
continue # unknown kinds are ignored
if kind in sections:
raise Refused("malformed", "duplicate section")
sections[kind] = data[start:end]
if SECT_CONTRACT not in sections:
raise Refused("malformed", "missing contract")
if SECT_OPLOG not in sections:
raise Refused("malformed", "missing op-log")
n_qubits = read_contract(gs, sections[SECT_CONTRACT])
ops, used = read_oplog(sections[SECT_OPLOG])
if gs is V2 and used != len(sections[SECT_OPLOG]):
raise Refused("malformed", "trailing bytes")
check_valid(gs, n_qubits, ops)
measure = read_measure(gs, sections[SECT_MEASURE]) if SECT_MEASURE in sections else None
return gs, n_qubits, ops, measure
# ---------------------------------------------------------------- simulator (§3.0)
def fx(a, b):
"""(a·b + 2^29) >> 30, exact, with an arithmetic shift."""
return (a * b + HALF) >> FRAC
def cmul(a, b):
return (fx(a[0], b[0]) - fx(a[1], b[1]), fx(a[0], b[1]) + fx(a[1], b[0]))
def cadd(a, b):
return (a[0] + b[0], a[1] + b[1])
def matrix(op, param, table, h):
z, one, i, minus_i = (0, 0), (ONE, 0), (0, ONE), (0, -ONE)
p3 = table[3]
return {
I: ((one, z), (z, one)),
X: ((z, one), (one, z)),
Y: ((z, minus_i), (i, z)),
Z: ((one, z), (z, (-ONE, 0))),
H: (((h, 0), (h, 0)), ((h, 0), (-h, 0))),
S: ((one, z), (z, i)),
SDG: ((one, z), (z, minus_i)),
T: ((one, z), (z, p3)),
TDG: ((one, z), (z, (p3[0], -p3[1]))),
P: ((one, z), (z, table[param] if op == P else z)),
}[op]
def simulate(n_qubits, ops, table):
"""The statevector from |0…0⟩. H is re P(3) of the same table."""
h = table[3][0]
dim = 1 << n_qubits
amps = [(0, 0)] * dim
amps[0] = (ONE, 0)
for op, controls, target, param in ops:
(m00, m01), (m10, m11) = matrix(op, param, table, h)
tbit = 1 << target
cmask = 0
for c in controls:
cmask |= 1 << c
for i in range(dim):
if i & tbit == 0 and i & cmask == cmask:
j = i | tbit
a, b = amps[i], amps[j]
amps[i] = cadd(cmul(m00, a), cmul(m01, b))
amps[j] = cadd(cmul(m10, a), cmul(m11, b))
return amps
# ---------------------------------------------------------------- the hashes (§4)
def circuit_hash(gs, n_qubits, ops):
return blake3(b"wai:quantum-circuit\x01" + gs.contract(n_qubits) + oplog_bytes(ops)).hex()
def statevector_hash(n_qubits, amps):
out = bytearray(b"wai:quantum-statevector\x01") + bytes([n_qubits]) + struct.pack("<I", FRAC)
for re, im in amps:
out += struct.pack("<qq", re, im)
return blake3(bytes(out)).hex()
def splitmix64(state):
state = (state + 0x9E3779B97F4A7C15) & 0xFFFFFFFFFFFFFFFF
z = state
z = ((z ^ (z >> 30)) * 0xBF58476D1CE4E5B9) & 0xFFFFFFFFFFFFFFFF
z = ((z ^ (z >> 27)) * 0x94D049BB133111EB) & 0xFFFFFFFFFFFFFFFF
return state, z ^ (z >> 31)
def histogram_hash(n_qubits, amps, seed, shots):
"""Each shot draws r = splitmix64 mod W, W the sum of the weights re² + im², and
counts the first basis state whose cumulative weight exceeds r."""
cumulative, total = [], 0
for re, im in amps:
total += re * re + im * im
cumulative.append(total)
counts = [0] * len(amps)
state = seed
for _ in range(shots):
state, r = splitmix64(state)
r %= total
lo, hi = 0, len(cumulative)
while lo < hi:
mid = (lo + hi) // 2
if cumulative[mid] > r:
hi = mid
else:
lo = mid + 1
counts[min(lo, len(counts) - 1)] += 1
out = bytearray(b"wai:quantum-histogram\x01") + bytes([n_qubits]) + struct.pack("<QQ", seed, shots)
for c in counts:
out += struct.pack("<Q", c)
return blake3(bytes(out)).hex()
# ---------------------------------------------------------------- the corpus
def qft_ops(n):
"""H on qubit j, then CP(k = l - j + 1; control l, target j) for l > j; then the
reversal swaps, each as CX(a→b), CX(b→a), CX(a→b)."""
ops = []
for j in range(n):
ops.append(gate(H, j))
for l in range(j + 1, n):
ops.append(gate(P, j, (l,), l - j + 1))
for j in range(n // 2):
a, b = j, n - 1 - j
ops += [gate(X, b, (a,)), gate(X, a, (b,)), gate(X, b, (a,))]
return ops
def dyadic_ladder_ops():
ops = [gate(X, 0)] + [gate(op, 1) for op in (Y, Z, S, SDG, T, TDG, I)]
chain = [gate(P, 1, (0,), k) for k in range(1, DYADIC_MAX + 1)]
return ops + chain + chain
# name: (n_qubits, ops, measure (seed, shots) or None). Both trees hold all of them.
VECTORS = {
"bell": (2, [gate(H, 0), gate(X, 1, (0,))], (12648430, 100000)),
"ghz3": (3, [gate(H, 0), gate(X, 1, (0,)), gate(X, 2, (1,))], None),
"qft3": (3, qft_ops(3), None),
"qft5": (5, qft_ops(5), None),
"qft5_x22": (5, [gate(X, 1), gate(X, 2), gate(X, 4)] + qft_ops(5), None),
"dyadic_ladder": (2, dyadic_ladder_ops(), None),
"qft8_ones": (8, [gate(X, q) for q in range(8)] + qft_ops(8), None),
}
# The four original circuit/1 vectors are frozen: circuit, statevector, histogram hash.
FROZEN_V1 = {
"bell": ("894d9a9582b19f27b8d08d06938357237fc00f7f30bf4bfc2be0c04fa5e47a4e",
"bd5df69662f4e61878d17baee1f3cf5705c3726d6740d1a4bf7be5133c099efb",
"0893f48ef32ca9f2865f6617537f0a4a5334f67e3c944bba845beb2f377acdb3"),
"ghz3": ("34e44f7d68037cbb3aee677f56b894dc48df6d408aa068a16d3b18076e5e7c5c",
"08c5e6fb50b91a00d28795506cc5d99ce1e9a3a8b4c41e4bfd05f5d829af9542", None),
"qft3": ("7e8414be4410557a05b52bcc1084e1c5f93b6caa0f28fa45f4d2ab661efd2a60",
"0fe1e5b8e2ac8e9afbfa926b3877967ebfead92efb7ed8e1ff6b117668926d6e", None),
"qft5": ("f36a9af99a251cae6f42144b2be3c57523f15e1e2ab3492cf15cac709fc048be",
"c268d77f2ddf84fecb028e80ce989123912870f04e5b9a1e3d885c9ce235c27d", None),
}
# The reject corpus: file stem → (class, the rule that refuses it here).
REJECTS = {
"magic_wqc3": ("malformed", "magic"),
"v1_ext_other": ("contract", "ext"),
"v1_numeric_other": ("contract", "numeric"),
"v1_frac_31": ("contract", "frac"),
"v1_frac_float": ("contract", "frac"),
"v1_gateset_rn": ("contract", "gateset"),
"v1_contract2_in_wqc1": ("contract", "ext"),
"v1_fields_missing": ("contract", "ext"),
"v1_contract_not_object": ("contract", "not an object"),
"v1_n_qubits_256": ("malformed", "n_qubits"),
"v1_n_qubits_258": ("malformed", "n_qubits"),
"v1_contract_twice": ("malformed", "duplicate section"),
"v1_oplog_twice": ("malformed", "duplicate section"),
"v1_contract_duplicate_member": ("contract", "not JSON: duplicate member"),
"v1_contract_lone_surrogate": ("contract", "not JSON: surrogate"),
"v1_contract_too_deep": ("contract", "not JSON: depth"),
"v1_n_qubits_minus0": ("malformed", "n_qubits"),
"v1_contract_lone_low_surrogate": ("contract", "not JSON: surrogate"),
"v1_contract_noncharacter": ("contract", "not JSON: noncharacter"),
"v1_contract_raw_noncharacter": ("contract", "not JSON: noncharacter"),
"v1_contract_integer_rounds_to_infinity": ("contract", "not JSON: number range"),
"v1_contract_number_rounds_to_infinity": ("contract", "not JSON: number range"),
"v1_contract_long_literal": ("contract", "not JSON: number range"),
"v1_contract_bom": ("contract", "not JSON: syntax"),
"v1_contract_nan": ("contract", "not JSON: syntax"),
"v1_contract_control_character": ("contract", "not JSON: syntax"),
"v1_measure_duplicate_member": ("malformed", "measure JSON: duplicate member"),
"v1_contract_not_utf8": ("contract", "not JSON: utf-8"),
"v1_contract_number_out_of_range": ("contract", "not JSON: number range"),
"v1_n_qubits_fraction": ("malformed", "n_qubits"),
"v1_n_qubits_missing": ("malformed", "n_qubits"),
"v1_measure_not_json": ("malformed", "measure JSON: syntax"),
"short_payload": ("malformed", "magic"),
"table_past_end": ("malformed", "section table"),
"section_past_end": ("malformed", "section bounds"),
"unknown_section_past_end": ("malformed", "section bounds"),
"no_contract": ("malformed", "missing contract"),
"no_oplog": ("malformed", "missing op-log"),
"v1_oplog_short": ("malformed", "op-log header"),
"v1_op_truncated_controls": ("malformed", "op controls"),
"v1_opcode_10": ("malformed", "opcode"),
"v1_target_out_of_range": ("invalid", "target"),
"v1_control_out_of_range": ("invalid", "control"),
"v1_control_is_target": ("invalid", "control is target"),
"v1_phase_0": ("invalid", "phase index"),
"v1_phase_33": ("invalid", "phase index"),
"v1_n_ops_u32max": ("malformed", "op header"),
"v1_offset_wrap": ("malformed", "section bounds"),
"v2_contract_whitespace": ("contract", "canonical contract"),
"v2_contract_key_order": ("contract", "canonical contract"),
"v2_contract_extra_member": ("contract", "canonical contract"),
"v2_contract1_in_wqc2": ("contract", "canonical contract"),
"v2_n_qubits_256": ("malformed", "n_qubits"),
"v2_contract_not_object": ("contract", "not an object"),
"v2_n_qubits_rounds_to_infinity": ("contract", "not JSON: number range"),
"v2_n_qubits_largest_finite": ("malformed", "n_qubits"),
"v2_phase_0": ("invalid", "phase index"),
"v2_measure_twice": ("malformed", "duplicate section"),
"v2_measure_seed_signed": ("malformed", "measure seed"),
"v2_duplicate_controls": ("invalid", "repeated control"),
"v2_param_on_h": ("invalid", "param"),
"v2_oplog_trailing_byte": ("malformed", "trailing bytes"),
"v2_opcode_10": ("malformed", "opcode"),
"v2_measure_spaced": ("malformed", "canonical measure"),
"v2_measure_no_shots": ("malformed", "measure shots"),
"v2_measure_basis_x": ("malformed", "canonical measure"),
}
# ---------------------------------------------------------------- the checks
class Report:
def __init__(self):
self.failures = 0
def check(self, label, problem):
"""`problem` is None for a pass, else why it failed."""
if problem is None:
print("PASS", label)
else:
self.failures += 1
print("FAIL %s: %s" % (label, problem))
def check_blake3(rep):
bad = [n for n, want in BLAKE3_VECTORS if blake3(bytes(i % 251 for i in range(n))).hex() != want]
rep.check("BLAKE3 against the published test vectors", "lengths %s differ" % bad if bad else None)
def check_tables(rep):
v1 = derive_v1()
rep.check("circuit/1 table: the recurrence gives Appendix B.1",
None if v1 == PHASE_V1 else "k = %s differ" % [k for k in range(1, 33) if v1[k] != PHASE_V1[k]])
d1 = table_digest(PHASE_V1)
rep.check("circuit/1 table digest", None if d1 == PHASE_V1_DIGEST else "got " + d1)
v2, problems = derive_v2()
rep.check("circuit2 table: every enclosure certifies its nearest integer", "; ".join(problems) or None)
rep.check("circuit2 table: the certified integers are Appendix B.2",
None if v2 == PHASE_V2 else "k = %s differ" % [k for k in range(1, 33) if v2[k] != PHASE_V2[k]])
d2 = table_digest(PHASE_V2)
rep.check("circuit2 table digest", None if d2 == PHASE_V2_DIGEST else "got " + d2)
APPENDIX_BEGIN = "<!-- generated: wai_quantum_vectors appendix; do not edit by hand -->"
APPENDIX_END = "<!-- end generated -->"
def read_appendix(path):
"""Appendix B of the extension as {subsection: (capability, table, digest)}, or the
reason it cannot be read. A subsection is a `### B.n `capability`` heading, rows
`| k | re | im |`, and a `Digest: `hex`` line."""
try:
with open(path, encoding="utf-8") as f:
text = f.read()
except OSError as e:
return None, str(e)
begin, end = text.find(APPENDIX_BEGIN), text.find(APPENDIX_END)
if begin < 0 or end < begin:
return None, "no generated Appendix B"
out, current = {}, None
for line in text[begin + len(APPENDIX_BEGIN):end].splitlines():
line = line.strip()
if line.startswith("### B."):
sub, _, cap = line[4:].partition(" ")
current = sub
out[sub] = (cap.strip("`"), [None], None)
elif current and line.startswith("| ") and line[2:3].isdigit():
try:
k, re_, im = (int(c.strip()) for c in line.strip("|").split("|"))
except ValueError:
return None, "%s: unreadable row %r" % (current, line)
table = out[current][1]
if k != len(table):
return None, "%s: row k = %d out of order" % (current, k)
table.append((re_, im))
elif current and line.startswith("Digest: `"):
out[current] = (out[current][0], out[current][1], line[len("Digest: `"):].rstrip("`"))
return out, None
def check_appendix(rep, path):
"""The extension's Appendix B states the literal tables and digests above."""
appendix, why = read_appendix(path)
if appendix is None:
rep.check("extension Appendix B (%s)" % path, why)
return
for sub, gs, digest in (("B.1", V1, PHASE_V1_DIGEST), ("B.2", V2, PHASE_V2_DIGEST)):
got = appendix.get(sub)
if got is None:
rep.check("extension Appendix %s" % sub, "missing")
continue
cap, table, stated = got
problems = []
if cap != gs.capability:
problems.append("names %s" % cap)
if table != gs.table:
problems.append("k = %s differ" % [k for k in range(1, DYADIC_MAX + 1)
if k >= len(table) or table[k] != gs.table[k]])
if stated != digest:
problems.append("digest %s" % stated)
rep.check("extension Appendix %s is the %s table and digest" % (sub, gs.label), "; ".join(problems) or None)
extra = sorted(set(appendix) - {"B.1", "B.2"})
rep.check("extension Appendix B has no other table", "has %s" % extra if extra else None)
def read_object(path):
"""A JSON object from `path`, or the reason it is not one."""
try:
with open(path, "rb") as f:
v = strict_json(f.read())
except OSError as e:
return None, str(e)
if not isinstance(v, dict):
return None, "not a JSON object"
return v, None
def listed_and_present(listed, directory):
"""None when `directory` holds exactly the `.wqc` files `listed` names."""
try:
present = {f for f in os.listdir(directory) if f.endswith(".wqc")}
except OSError as e:
return str(e)
listed = set(listed)
if listed == present:
return None
return "listed but missing %s, present but unlisted %s" % (sorted(listed - present), sorted(present - listed))
def expected_entry(gs, name, label_capability):
"""The manifest entry for vector `name` in `gs`, computed here, and the hashes."""
n, ops, measure = VECTORS[name]
amps = simulate(n, ops, gs.table)
entry = {"n_qubits": n, "n_ops": len(ops), "circuit_hash": circuit_hash(gs, n, ops),
"statevector_hash": statevector_hash(n, amps)}
if label_capability:
entry["capability"] = gs.capability
if measure is not None:
entry["measure"] = {"seed": measure[0], "shots": measure[1],
"histogram_hash": histogram_hash(n, amps, *measure)}
return entry
def check_tree(rep, root, sub, gs, label_capability):
tree = os.path.join(root, sub) if sub else root
prefix = sub + "/" if sub else ""
manifest, why = read_object(os.path.join(tree, "expected_hashes.json"))
if manifest is None:
rep.check(prefix + "expected_hashes.json", why)
return
names = ["%s.wqc" % v for v in VECTORS]
rep.check(prefix + "expected_hashes.json lists every vector and no other",
None if set(manifest) == set(names) else "lists %s" % sorted(manifest))
rep.check(prefix + "vectors/ holds exactly the listed files", listed_and_present(names, os.path.join(tree, "vectors")))
for name, (n, ops, measure) in VECTORS.items():
label = "%s%s.wqc (%s)" % (prefix, name, gs.label)
try:
with open(os.path.join(tree, "vectors", name + ".wqc"), "rb") as f:
data = f.read()
except OSError as e:
rep.check(label, str(e))
continue
problems = []
if data != write_wqc(gs, n, ops, measure):
problems.append("the file is not the bytes rebuilt here")
try:
if read_wqc(data) != (gs, n, ops, measure):
problems.append("the file reads back as another circuit")
except Refused as r:
problems.append("refused (%s)" % r)
want = expected_entry(gs, name, label_capability)
got = manifest.get(name + ".wqc")
if got != want:
differ = sorted(k for k in set(want) | set(got or {}) if (got or {}).get(k) != want.get(k))
problems.append("manifest members %s differ from the values computed here" % differ)
if gs is V1 and name in FROZEN_V1:
frozen = FROZEN_V1[name]
computed = (want["circuit_hash"], want["statevector_hash"],
want["measure"]["histogram_hash"] if "measure" in want else None)
if computed != frozen:
problems.append("a frozen circuit/1 hash changed")
rep.check(label, "; ".join(problems) or None)
def check_rejects(rep, root):
tree = os.path.join(root, "reject")
listing, why = read_object(os.path.join(tree, "reject.json"))
if listing is None:
rep.check("reject/reject.json", why)
return
want = {"%s.wqc" % stem: kind for stem, (kind, _) in REJECTS.items()}
rep.check("reject/reject.json lists every reject case with its class",
None if listing == want else "differs in %s" % sorted(
f for f in set(listing) | set(want) if listing.get(f) != want.get(f)))
rep.check("reject/ holds exactly the listed files", listed_and_present(listing, tree))
for stem, (kind, rule) in REJECTS.items():
label = "reject/%s.wqc" % stem
try:
with open(os.path.join(tree, stem + ".wqc"), "rb") as f:
data = f.read()
except OSError as e:
rep.check(label, str(e))
continue
try:
read_wqc(data)
rep.check(label, "accepted; it must be refused as %s (%s)" % (kind, rule))
except Refused as r:
rep.check(label + " refused as %s (%s)" % (kind, rule),
None if (r.kind, r.rule) == (kind, rule) else "refused as %s (%s)" % (r.kind, r.rule))
# The accept corpus: the file stems `accept.json` must list. Each file's reading is
# recomputed here and compared with the listing.
ACCEPTS = [
"v1_members_reordered", "v1_whitespace", "v1_further_members", "v1_escapes_and_characters",
"v1_depth_127", "v1_numbers_largest_finite", "v1_numbers_long_literals", "v1_numbers_toward_zero",
"v1_n_qubits_255",
"v1_measure_lax", "v1_measure_missing_shots", "v1_measure_u64_max", "v1_measure_not_u64",
"v1_measure_other_types", "v1_measure_not_object", "v1_oplog_trailing_bytes", "v1_param_on_h",
"v1_repeated_control", "v1_unknown_sections", "v2_unknown_sections", "v2_n_qubits_255",
"v2_measure_u64_max",
]
def check_accepts(rep, root):
tree = os.path.join(root, "accept")
listing, why = read_object(os.path.join(tree, "accept.json"))
if listing is None:
rep.check("accept/accept.json", why)
return
want = sorted("%s.wqc" % stem for stem in ACCEPTS)
rep.check("accept/accept.json lists every accept case",
None if sorted(listing) == want else "differs in %s" % sorted(set(listing) ^ set(want)))
rep.check("accept/ holds exactly the listed files", listed_and_present(listing, tree))
for file in want:
label = "accept/" + file
try:
with open(os.path.join(tree, file), "rb") as f:
data = f.read()
gs, n_qubits, ops, measure = read_wqc(data)
except OSError as e:
rep.check(label, str(e))
continue
except Refused as r:
rep.check(label, "refused as %s (%s); it must be accepted" % (r.kind, r.rule))
continue
got = {
"capability": gs.capability,
"n_qubits": n_qubits,
"n_ops": len(ops),
"circuit_hash": circuit_hash(gs, n_qubits, ops),
"measure": None if measure is None else {"seed": measure[0], "shots": measure[1]},
}
rep.check(label + " accepted and read as listed",
None if got == listing.get(file) else "read as %s, listed as %s" % (got, listing.get(file)))
def check_mutations(rep):
n, ops, _ = VECTORS["dyadic_ladder"]
for gs in (V1, V2):
base = statevector_hash(n, simulate(n, ops, gs.table))
missed = []
for k in range(1, DYADIC_MAX + 1):
for component in (0, 1):
for delta in (-1, 1):
table = list(gs.table)
entry = list(table[k])
entry[component] += delta
table[k] = tuple(entry)
if statevector_hash(n, simulate(n, ops, table)) == base:
missed.append("P(%d).%s%+d" % (k, "re" if component == 0 else "im", delta))
rep.check("dyadic_ladder sees every ±1 change of P(1..32) in the %s table" % gs.label,
"missed %s" % missed if missed else None)
def main():
root = sys.argv[1] if len(sys.argv) > 1 else os.path.join(os.path.dirname(os.path.abspath(__file__)), "..")
extension = sys.argv[2] if len(sys.argv) > 2 else os.path.join(root, "..", "extensions", "quantum-sim.md")
rep = Report()
check_blake3(rep)
check_tables(rep)
check_appendix(rep, extension)
check_tree(rep, root, "", V1, False)
check_tree(rep, root, "v2", V2, True)
check_rejects(rep, root)
check_accepts(rep, root)
check_mutations(rep)
if rep.failures:
print("%d failed — NOT conformant" % rep.failures)
sys.exit(1)
print("referee: every check passes")
if __name__ == "__main__":
main()