diff --git a/research/bitset-gauss-elimination/BENCHMARK.json b/research/bitset-gauss-elimination/BENCHMARK.json new file mode 100644 index 0000000000..6f1c366116 --- /dev/null +++ b/research/bitset-gauss-elimination/BENCHMARK.json @@ -0,0 +1,137 @@ +{ + "errors": 0, + "executable": "/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12", + "platform": "macOS-27.0.1-arm64-arm-64bit", + "python": "3.12.15 (main, Oct 14 2026, 00:55:37) [Clang 23.1.3 ]", + "samples": [ + { + "mode": "baseline", + "returncode": 0, + "seconds": 0.44511516694910824, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "baseline" + } + }, + { + "mode": "bitset", + "returncode": 0, + "seconds": 0.2640917918179184, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "bitset" + } + }, + { + "mode": "baseline", + "returncode": 0, + "seconds": 0.46423900010995567, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "baseline" + } + }, + { + "mode": "bitset", + "returncode": 0, + "seconds": 0.2641322501003742, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "bitset" + } + }, + { + "mode": "baseline", + "returncode": 0, + "seconds": 0.43559599993750453, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "baseline" + } + }, + { + "mode": "bitset", + "returncode": 0, + "seconds": 0.2659395409282297, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "bitset" + } + }, + { + "mode": "baseline", + "returncode": 0, + "seconds": 0.4352738750167191, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "baseline" + } + }, + { + "mode": "bitset", + "returncode": 0, + "seconds": 0.3204354590270668, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "bitset" + } + }, + { + "mode": "baseline", + "returncode": 0, + "seconds": 0.5051918339449912, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "baseline" + } + }, + { + "mode": "bitset", + "returncode": 0, + "seconds": 0.26453037513419986, + "stdout_sha256": "9e9ca672563db4c1d3a281be16c2fdcd798ae96de4fd9cd8c2e9ca1ce624effa", + "work": { + "gauss_calls": 25, + "maximum_dimension": 128, + "mode": "bitset" + } + } + ], + "source_hashes": { + "benchmark.py": "01b81ef23003aba719a611bc312efa44a3f955bc3118afc6bac9c252c975a94f", + "fast_gauss.py": "85bcce68ff3233b63a0d1f5cc0bae7efd596ec9bf7973cf1b8f278d43851fb34", + "run_projection.py": "053091c8b84e703ed622440a9b84af100171ebe8a7cc578e28a40631e94ba50c", + "run_replay.py": "c7e2d139c26bef270f9bc547aa0a5b6b5dd9faaf732123498e8ddd117889203e" + }, + "summary": { + "baseline": { + "maximum": 0.5051918339449912, + "median": 0.44511516694910824, + "minimum": 0.4352738750167191 + }, + "bitset": { + "maximum": 0.3204354590270668, + "median": 0.26453037513419986, + "minimum": 0.2640917918179184 + } + }, + "workload": "projection" +} diff --git a/research/bitset-gauss-elimination/PINS.json b/research/bitset-gauss-elimination/PINS.json new file mode 100644 index 0000000000..390ebce785 --- /dev/null +++ b/research/bitset-gauss-elimination/PINS.json @@ -0,0 +1,22 @@ +{ + ".gitattributes": "e8de51f40e73a262bbd8959a4230e3dce0b2b15ad68a19888750b7fc26593cf3", + "INPUTS.json": "c62bd55fd36c136ba2ec34f6d91ed481497932fadc7134d0b7e5c10009fd588c", + "PROOF.md": "b3834a57ba82f4fe316669910693b6c2399f50af8e7a73722d61cddc00ec31aa", + "README.md": "3a01996ff6da7402442842db70c45ea6aaceab1ae03c17dad1fcb2e11fee1f89", + "RESULT.json": "ce586c74a8633e2907b3227c0cbcff069e08e2d2afe4695c8d56abc844dff97f", + "SOURCES.json": "c8e43472b247d62b17ca74898d05945d36bd590db650c7a2101ded6bc48116c7", + "TASK.json": "3c302a0c814c47ecc0c7fb64b50872f4f3f5ed18dfdb65976e8f19ae6daadb20", + "VERIFY.json": "3940455d0e2e8195a7a8cfcd6504d07f76f25ace824b1d77d6f374d2ce3e3348", + "contractions.py": "a2fc065d6f293ddcffe5b6eec7fff0fdf57c8e3d4c6be355bd92e87d3f1865ef", + "source/CONSTRUCTION.json.gz": "00d13d7774adf087f93bc58d4c5860e296f8408c0ee25a4db254a056f8da7dec", + "source/LICENSE": "c71d239df91726fc519c6eb72d318ec65820627232b2f796219e87dcf35d0ab4", + "source/PHYSICAL.json": "f78ba861d23ff33f2c9ffeca8c3ab6a7f1ec906f38adff5d325c13efb187c7c3", + "source/PR5526-PROOF.md": "45483e4b5f604a5e70602c5caaa4394d4a0d0d1a6d0b6e8a93b836d50b23f764", + "source/PR5698-NOTICE.md": "1f55b3e4bb80a60611b1fca0622da5a99fe7877446f03d0edd1dc8f3a5abadd2", + "source/PR5704-NOTICE.md": "a77adc71861e3b9c73c703c5db5b671f61a1101d85d14cca08547e21d6d8dce4", + "source/PR5704-PROOF.md": "e5039414a02e36cdeab82fb9803fc4de0c3e0f3997e8812ef9d227e5d417ec91", + "source/PR5704-quadratic.py": "d406584aed1479b75c387bc964b9e797b2f4f75e5f8c6bbea00df07f7d590a91", + "source/clifford_exact.py": "2982ead07721e83358dc67c6f49c2b4c16a236a6bebdfd72f58221be1f5d0962", + "source/prune_native.py": "f2b37c183a10922eabc8e0436b426a08866430e96bd0adc7e178951c7fbb5257", + "verify.py": "be69f731a2b0ff7a4d47c90e49a9466ed2ab38a571235eff6953891f9707e53c" +} diff --git a/research/bitset-gauss-elimination/PROOF.md b/research/bitset-gauss-elimination/PROOF.md new file mode 100644 index 0000000000..7148214509 --- /dev/null +++ b/research/bitset-gauss-elimination/PROOF.md @@ -0,0 +1,107 @@ +# Equivalence of adjacency-bitset elimination + +This implementation is produced by Nick Petalas. Quadratic Gauss-sum +elimination is classical. Ali Romman's PR5727 retains its contraction, +phase, state and finite-span contracts. The contribution here is this +implementation and its consumer and performance evidence. + +The admitted input is a canonical polynomial on at most 256 Boolean variables, + +```text +q(x) = c + sum a_i x_i + 2 sum_(i> k & 1: + diagonal[k] = (diagonal[k] - coefficient) % 4 + adjacency[k] ^= neighbors ^ (1 << k) + + +def pair_step(diagonal, adjacency, left, right, coefficient, target, n): + """Substitute y=target XOR parity(left) after summing the paired x.""" + for k in range(n): + in_left, in_right = left >> k & 1, right >> k & 1 + diagonal[k] = (diagonal[k] + coefficient * (1 - 2 * target) * in_left + + 2 * target * in_right + 2 * (in_left & in_right)) % 4 + delta = (right if in_left else 0) ^ (left if in_right else 0) + if coefficient & 1 and in_left: + delta ^= left ^ (1 << k) + adjacency[k] ^= delta & ~(1 << k) + + +def gauss(n, poly): + """Return sum_x i**q(x), without affine-polynomial substitutions.""" + constant, diagonal, adjacency = unpack(n, poly) + active, real, imaginary = (1 << n) - 1, 1, 0 + for _ in range(n): + if not active: + break + i = active.bit_length() - 1 + active ^= 1 << i + coefficient = diagonal[i] + neighbors = adjacency[i] & active + if coefficient & 1: + sign = 1 if coefficient == 1 else -1 + real, imaginary = real - sign * imaginary, imaginary + sign * real + odd_step(diagonal, adjacency, neighbors, coefficient, n) + elif neighbors: + j = neighbors.bit_length() - 1 + active ^= 1 << j + left, right = neighbors & active, adjacency[j] & active + target, other = coefficient // 2, diagonal[j] + constant = (constant + other * target) % 4 + pair_step(diagonal, adjacency, left, right, other, target, n) + real, imaginary = 2 * real, 2 * imaginary + elif coefficient: + return 0, 0 + else: + real, imaginary = 2 * real, 2 * imaginary + require(active == 0, "Every variable was eliminated") + return ((real, imaginary), (-imaginary, real), + (-real, -imaginary), (imaginary, -real))[constant] diff --git a/research/bitset-gauss-elimination/run_projection.py b/research/bitset-gauss-elimination/run_projection.py new file mode 100644 index 0000000000..1bcc93f04d --- /dev/null +++ b/research/bitset-gauss-elimination/run_projection.py @@ -0,0 +1,58 @@ +"""A fixed complete 64-axis, four-source Clifford projection calculation.""" +import argparse +import json +import random +import runpy +import sys +from pathlib import Path + +ROOT = Path(__file__).resolve().parent +sys.path.insert(0, str(ROOT)) +from fast_gauss import gauss, require +from run_replay import authenticate, SOURCE + + +def circuits(n=64): + require(type(n) is int and 1 <= n <= 64, "Workload width in 1..64") + rng, result = random.Random(932 + n), [] + for _ in range(5): + rows = [0] * n + for i in range(n): + for j in range(i + 1): + if rng.randrange(2): + rows[i] ^= 1 << j + if i != j: + rows[j] ^= 1 << i + linear = tuple((rows[i] >> i & 1) + 2 * rng.randrange(2) for i in range(n)) + result.append([("C", 1 << i, 1) for i in range(n)] + + [("D", tuple(rows), linear)]) + return result + + +def run(mode): + authenticate() + source = runpy.run_path(str(SOURCE / "verify.py"), run_name="projection_source") + engine, original = source["C"], source["C"].gauss + counts = {"gauss_calls": 0, "maximum_dimension": 0} + + def replacement(n, poly): + counts["gauss_calls"] += 1 + counts["maximum_dimension"] = max(counts["maximum_dimension"], n) + value = original(n, poly) if mode == "baseline" else gauss(n, poly) + if mode == "differential": + require(value == original(n, poly), "Projection Gauss differential") + return value + + engine.gauss = replacement + states = [engine.choi(64, circuit) for circuit in circuits()] + result = engine.best_approximation(states[:4], states[4], 1) + require(result["span_rank"] == 4, "Four independent source operators") + print(json.dumps(result, default=source["encode"], sort_keys=True)) + print("PASS_BITSET_GAUSS_PROJECTION") + print(json.dumps(dict(mode=mode, **counts), sort_keys=True), file=sys.stderr) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--mode", choices=("baseline", "bitset", "differential"), required=True) + run(parser.parse_args().mode) diff --git a/research/bitset-gauss-elimination/run_replay.py b/research/bitset-gauss-elimination/run_replay.py new file mode 100644 index 0000000000..576b8ca721 --- /dev/null +++ b/research/bitset-gauss-elimination/run_replay.py @@ -0,0 +1,57 @@ +"""Run the unchanged complete contraction replay with one Gauss-sum replacement.""" +import argparse +from hashlib import sha256 +import json +from pathlib import Path +import runpy +import sys + +ROOT = Path(__file__).resolve().parent +sys.path.insert(0, str(ROOT)) +from fast_gauss import gauss, require + +SOURCE = ROOT.parent / "clifford-exact-contractions" + + +def authenticate(): + pin_file = ROOT / "PINS.json" + require(pin_file.stat().st_size <= 65536, "Bounded pin manifest") + records = json.loads(pin_file.read_text()) + require(1 <= len(records) <= 64, "Bounded source inventory") + for name, digest in records.items(): + relative = Path(name) + require(not relative.is_absolute() and ".." not in relative.parts, + "Repository-local source name") + path = SOURCE / name + require(path.is_file() and not path.is_symlink(), "Regular source input") + require(path.stat().st_size <= 1048576, "Source input at most 1 MiB") + require(sha256(path.read_bytes()).hexdigest() == digest, "Changed source: " + name) + + +def run(mode): + authenticate() + source = runpy.run_path(str(SOURCE / "verify.py"), run_name="contraction_replay") + original = source["C"].gauss + counts = {"gauss_calls": 0, "maximum_dimension": 0} + + def replacement(n, poly): + counts["gauss_calls"] += 1 + counts["maximum_dimension"] = max(counts["maximum_dimension"], n) + value = original(n, poly) if mode == "baseline" else gauss(n, poly) + if mode == "differential": + require(value == original(n, poly), "Gauss differential mismatch") + return value + + source["C"].gauss = replacement + result = source["run"]() + expected = json.loads((SOURCE / "RESULT.json").read_text()) + require(result == expected, "Entire committed predecessor result") + print(json.dumps(result, indent=2, sort_keys=True)) + print("PASS_EXACT_CLIFFORD_CONTRACTIONS") + print(json.dumps(dict(mode=mode, **counts), sort_keys=True), file=sys.stderr) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--mode", choices=("baseline", "bitset", "differential"), required=True) + run(parser.parse_args().mode) diff --git a/research/bitset-gauss-elimination/verify.py b/research/bitset-gauss-elimination/verify.py new file mode 100644 index 0000000000..082ef82d44 --- /dev/null +++ b/research/bitset-gauss-elimination/verify.py @@ -0,0 +1,126 @@ +"""Bounded independent literal oracles and complete consumer conformance.""" +import argparse +from hashlib import sha256 +from itertools import product +import json +from math import comb +from pathlib import Path +import random +import subprocess +import sys + +ROOT = Path(__file__).resolve().parent +sys.path.insert(0, str(ROOT)) +from fast_gauss import gauss, require + +UNITS = ((1, 0), (0, 1), (-1, 0), (0, -1)) + + +def literal(n, poly): + real, imaginary = 0, 0 + for x in range(1 << n): + phase = sum(c for mask, c in poly.items() if x & mask == mask) % 4 + a, b = UNITS[phase] + real, imaginary = real + a, imaginary + b + return real, imaginary + + +def exhaustive(): + count = 0 + for n in range(5): + masks = [1 << i for i in range(n)] + [ + (1 << i) | (1 << j) for i in range(n) for j in range(i)] + domains = [range(4)] * n + [(0, 2)] * (len(masks) - n) + for values in product(*domains): + poly = {mask: value for mask, value in zip(masks, values) if value} + base = literal(n, poly) + for constant in range(4): + if constant: + poly[0] = constant + else: + poly.pop(0, None) + a, b = base + expected = ((a, b), (-b, a), (-a, -b), (b, -a))[constant] + require(gauss(n, poly) == expected, "Complete literal quadratic sum") + count += 1 + return count + + +def larger(): + rng, count = random.Random(102932), 0 + for n in range(5, 13): + for _ in range(8): + poly = {1 << i: rng.randrange(1, 4) for i in range(n)} + for i in range(n): + for j in range(i): + if rng.randrange(2): + poly[(1 << i) | (1 << j)] = 2 + require(gauss(n, poly) == literal(n, poly), "Larger literal sum") + count += 1 + for n in (0, 1, 2, 4, 8, 16, 32, 64, 128, 256): + for a in range(4): + poly = {(1 << i) | (1 << j): 2 for i in range(n) for j in range(i)} + if a: + poly.update({1 << i: a for i in range(n)}) + real = imaginary = 0 + for weight in range(n + 1): + x, y = UNITS[(a * weight + weight * (weight - 1)) % 4] + real += comb(n, weight) * x + imaginary += comb(n, weight) * y + require(gauss(n, poly) == (real, imaginary), "Independent binomial dense sum") + count += 1 + return count + + +def controls(): + invalid = [(-1, {}), (257, {}), (True, {}), (1, []), + (1, {2: 2}), (1, {True: 1}), (1, {1: 0}), + (1, {1: True}), (2, {3: 1}), (3, {7: 2})] + for n, poly in invalid: + try: + gauss(n, poly) + except ValueError: + pass + else: + raise ValueError("Invalid polynomial was accepted") + witnesses = [(1, {1: 2}, (0, 0)), (1, {1: 3}, (1, -1)), + (2, {3: 2}, (2, 0)), (0, {0: 3}, (0, -1))] + for n, poly, expected in witnesses: + require(gauss(n, poly) == expected, "Cancellation, odd sign, pair and phase controls") + return len(invalid), len(witnesses) + + +def consumer(program, mode): + command = [sys.executable, "-I", "-W", "error", "-B"] + if sys.flags.optimize: + command.append("-O") + command += [str(ROOT / program), "--mode", mode] + result = subprocess.run(command, check=False, capture_output=True, timeout=120) + require(result.returncode == 0, result.stderr.decode()) + return dict(stdout_sha256=sha256(result.stdout).hexdigest(), work=json.loads(result.stderr)) + + +def main(quick): + result = dict(exhaustive=exhaustive(), larger=larger(), controls=controls()) + result["projection"] = consumer("run_projection.py", "differential") + require(result["projection"]["work"]["gauss_calls"] == 25, "Complete projection work") + if not quick: + result["complete_replay"] = consumer("run_replay.py", "differential") + require(result["complete_replay"]["work"]["gauss_calls"] == 82237, + "Complete differential replay work") + benchmark = json.loads((ROOT / "BENCHMARK.json").read_text()) + for name, digest in benchmark["source_hashes"].items(): + require(sha256((ROOT / name).read_bytes()).hexdigest() == digest, + "Changed benchmark implementation: " + name) + require(all(row["stdout_sha256"] == result["projection"]["stdout_sha256"] + for row in benchmark["samples"]), "Timed complete projection output") + result["package_sha256"] = {p.name: sha256(p.read_bytes()).hexdigest() + for p in sorted(ROOT.iterdir()) if p.is_file()} + print(json.dumps(result, indent=2, sort_keys=True)) + print("PASS_BITSET_GAUSS_ELIMINATION") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--quick", action="store_true") + main(parser.parse_args().quick)