From d9219ca07e69ca08c6cc7ec9dc805955fe87d857 Mon Sep 17 00:00:00 2001 From: Nick Petalas Date: Wed, 11 Nov 2026 13:26:56 +0000 Subject: [PATCH] feat(research): attain exact unit Stein coefficient bounds --- research/exact-unit-stein-splitting/PROOF.md | 131 +++++++++++ research/exact-unit-stein-splitting/README.md | 86 ++++++++ .../exact-unit-stein-splitting/SOURCES.json | 45 ++++ research/exact-unit-stein-splitting/TASK.json | 96 ++++++++ .../exact-unit-stein-splitting/VERIFY.json | 34 +++ .../exact-unit-stein-splitting/algebra.py | 158 +++++++++++++ research/exact-unit-stein-splitting/verify.py | 207 ++++++++++++++++++ 7 files changed, 757 insertions(+) create mode 100644 research/exact-unit-stein-splitting/PROOF.md create mode 100644 research/exact-unit-stein-splitting/README.md create mode 100644 research/exact-unit-stein-splitting/SOURCES.json create mode 100644 research/exact-unit-stein-splitting/TASK.json create mode 100644 research/exact-unit-stein-splitting/VERIFY.json create mode 100644 research/exact-unit-stein-splitting/algebra.py create mode 100644 research/exact-unit-stein-splitting/verify.py diff --git a/research/exact-unit-stein-splitting/PROOF.md b/research/exact-unit-stein-splitting/PROOF.md new file mode 100644 index 0000000000..9eaba89aec --- /dev/null +++ b/research/exact-unit-stein-splitting/PROOF.md @@ -0,0 +1,131 @@ +# Exact unit bounds by splitting orthogonal residuals + +This contribution was produced by Nick Petalas. It is a scoped written +proof over Q(i), with a bounded exact implementation and finite +corroboration. It is not a complete multiplication theorem. + +## Statement + +For every E in Q(i)^(m by n), of rank r, there are s<=2r columns g_j in +Q(i)^m and rows a_j in Q(i)^n such that + + E = sum_j g_j a_j, + ||g_j||_2 <= 1, + ||a_j||_2 <= ||E||_2. + +In fact, with G having columns g_j and A having rows a_j, + + GG* = P_col(E), A = G*E. + +Zero E has s=0. All nonzero emitted terms are retained. The bound 2r is +an upper bound, not an asserted global minimum. The coefficient row uses +an adjoint exactly once, in G*E. The standard coordinate factorization +uses m terms with the same unit bounds. The present bound depends on r +and is useful for low-rank displacements in a larger ambient space. + +## Construction and proof + +Use the exact nonzero orthogonal residuals v_1,...,v_r of E's columns, +as in Jose Manuel / jmbarrancoml's [PR4198](https://github.com/CrocSwap/integer-mult-bounds/pull/4198), +`research/rational-stein-generators/PROOF.md`, section 1. Write +n_k=v_k*v_k>0. The inherited projection identity is + + sum_k v_k v_k* / n_k = P_col(E). + +For each k, choose z_k1,z_k2 in Q(i) with + + |z_k1|^2 + |z_k2|^2 = n_k. + +Such a choice always exists. For example, clear a common positive integer +denominator D from v_k, and write v_k=w/D with Gaussian-integer w. +The positive integer N=w*w is a sum of four integer squares a^2+b^2+c^2+d^2 +by Lagrange's theorem. Then z_k1=(a+bi)/D and z_k2=(c+di)/D work. +The classical theorem is not a new claim here; a primary formal source +is [Mathlib's SumFourSquares module](https://leanprover-community.github.io/mathlib4_docs/Mathlib/NumberTheory/SumFourSquares.html). +No Lean replay of this package is claimed. + +Define, omitting each zero z, + + g_kl = v_k z_kl / n_k, a_kl = g_kl* E. + +There are at most two terms per residual, and + + sum_l g_kl g_kl* + = v_k v_k* (sum_l |z_kl|^2) / n_k^2 + = v_k v_k* / n_k. + +Thus GG*=P_col(E), and GA=GG*E=E. Also + + ||g_kl||_2^2 = |z_kl|^2/n_k <= 1, + ||a_kl||_2 = ||g_kl* E||_2 <= ||g_kl||_2 ||E||_2 <= ||E||_2. + +This proves the stated arbitrary-dimensional rational construction. +It also proves the stronger positive-semidefinite inequality + + E*E - a_kl* a_kl = E*(I-g_kl g_kl*)E >= 0. + +If N is a sum of two integer squares, the implementation uses one term +for that residual. Its finite search finds such a representation whenever +one exists within the complete range 0..floor(sqrt(N)). Otherwise the +four-square representation gives two nonzero terms. No optimality of the +aggregate count is inferred from this per-residual choice. + +## Stein and dirty-state interface + +For a square E=T-ZTZ*, with Z the nilpotent lower shift, finite telescoping +gives the credited identity + + T = sum_(h=0)^(m-1) Z^h E (Z*)^h + = sum_j L(g_j) U(a_j). + +The second equality is the coefficient-row convention of PR4198. Each +new term can be passed to its existing compensated dirty-source macro, +with that macro's own legal near-identity normalizer and all four calls +paid. Summing complete macros restores the arbitrary source bank and adds +Tx to the arbitrary target bank. This is a logical exact interface only. +We have not run those child macros or implemented physical coefficient +pages in this package. + +PR4198 uses r terms with row norms below 2||E||. Nick Petalas's +[PR6228](https://github.com/CrocSwap/integer-mult-bounds/pull/6228) +already gives r terms below (1+2^-p)||E|| and an exact-factor-one +obstruction in the rank-minimal class. Those results are credited +background, not replay targets. The new step here changes the term-count +contract and attains the exact factor one. + +## Encoding and implementation scope + +For a residual with rational components of height H, one common +denominator D uses at most 2mH bits. The Gaussian-integer coordinates of w +then use at most (2m+1)H+1 bits. Each four-square coordinate has magnitude +at most sqrt(N), so its bit length is bounded by half the bit length of N, +rounded upward. Forming g and a uses rational products and finite sums. +The construction therefore has finite rational output of polynomial +encoding length, using the inherited rational residual height bound. +This statement does not assert polynomial bit running time for our +exhaustive four-square enumerator. + +The implementation first removes the common integer gcd of w. This only +changes D's rational scale and leaves N smaller. It accepts at most eight +rows and eight columns, input component height at most 256, and primitive +N<=65,536 at every residual. Capacity exhaustion rejects explicitly. +The theorem has no such finite limits; the implementation and its finite +checks do. + +Every square loop is bounded by floor(sqrt(N))+1, its pair table has +N+1 entries, and its split loop has N+1 iterations. Matrix operations +have at most sixteen term columns and eight data rows or columns. There +is no recursion, external solver or predecessor execution. Managed Python +allocation is a Power of Ten rule 3 deviation, bounded by those capacities. +Warnings are treated as errors in every registered command. The source +and manifest checks are separate from mathematical verification. + +## Limits + +The construction can double the number of compensated generator macros. +It proves no decrease in total gate count, precision, memory, runtime or +recursive moment. Gaussian-rational coefficients may have odd denominators. +Exact coefficient preparation, rounded children, physical delivery, +restoration under rounding, compiler/native admission, stock/tape costs +and the outer theorem retain their existing obligations. The selected +frontier and all five framework assumptions are unchanged. diff --git a/research/exact-unit-stein-splitting/README.md b/research/exact-unit-stein-splitting/README.md new file mode 100644 index 0000000000..af45308272 --- /dev/null +++ b/research/exact-unit-stein-splitting/README.md @@ -0,0 +1,86 @@ +# Exact factor-one rational Stein generators + +Scoped written proof with exact finite corroboration. For every rectangular +Gaussian-rational matrix E of rank r, the new construction gives + + E = sum_(j=1)^s g_j a_j, s <= 2r, + ||g_j||_2 <= 1, ||a_j||_2 <= ||E||_2. + +The columns form a rational Parseval frame for col(E), and a_j=g_j*E. +This improves the coefficient norm contract of Jose Manuel's +[PR4198](https://github.com/CrocSwap/integer-mult-bounds/pull/4198) +from a factor below two to exact factor one by permitting extra terms. +Nick Petalas's [PR6228](https://github.com/CrocSwap/integer-mult-bounds/pull/6228) +already approaches one with r terms. Its fixed-term normalizer is not +used here. No new kappa is claimed. The standard coordinate factorization already +gives the unit bounds with m terms; this construction bounds the count +by 2r independently of the ambient row dimension m. + +The new step represents each rational residual norm as the sum of two +Gaussian-rational squared norms, then splits its orthogonal projector. +The rational residual interface and Stein/dirty-shear identities are +credited to PR4198; Lagrange's four-square theorem is classical. +[PROOF.md](PROOF.md) gives the arbitrary-dimensional argument and limits. + +## Evidence + +The verifier passes 414 matrix fixtures, containing total rank 698 and +901 emitted terms. Exactly 123 fixtures use extra terms. It checks 2,048 +positive integer square decompositions and seventeen rejection controls. +For E=(1,1,1)^T, it emits two terms with squared column norms 1/3 and 2/3 +and squared row norms 1 and 2. The exact matrix norm squared is three. + +The standalone checker independently row-reduces every input, reconstructs +every matrix entry, checks the full Parseval projector and verifies +E*E-a_j*a_j is positive semidefinite by exact Hermitian elimination. +It also checks the complete two-square preference against a separate +small enumeration. Counts and the explicit rank-one witness are emitted +by the verifier, with hashes of every package source file. No receipt is +handwritten and no predecessor scientific verifier is executed. + +## Reproduce + +From the repository root, using CPython 3.12 and its standard library: + +```sh +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -I -W error -B research/exact-unit-stein-splitting/verify.py +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -I -O -W error -B research/exact-unit-stein-splitting/verify.py +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -B scripts/research_package.py --package research/exact-unit-stein-splitting +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -B scripts/check_research_commit.py --package research/exact-unit-stein-splitting --check-source-history +``` + +Any CPython 3.12 path can replace the local interpreter. Expected markers +are PASS_EXACT_UNIT_STEIN_SPLITTING, PASS_REGISTERED_RESEARCH_PACKAGE and +PASS_RESEARCH_COMMIT_BYTES. The default replay is read-only. Optional +`--receipt /absolute/new/path.json` writes a source-bound receipt and +refuses to overwrite an existing file. + +## Limitations + +The written theorem has arbitrary finite dimensions. The executable +accepts dimensions at most eight, component height at most 256, and +primitive residual squared norm at most 65,536. It rejects exhaustion. +Its exhaustive integer-square search is not a polynomial bit-time claim. +The 2r term bound is not claimed optimal. Additional compensated macros +may outweigh the smaller coefficients. No physical precision or supplier +saving, complete multiplier, human review or formal-kernel acceptance is +claimed. All five framework assumptions remain. + +## Open next questions + +- Can a rational Parseval construction always use fewer than 2r terms? + Start with the embedding of a rank-two rational Hermitian form into a + three-coordinate standard form. This could reduce the added macros; + arbitrary rational-form embedding and its preparation costs are open here. +- Can an actual rounded child benefit after all extra terms are charged? + Start with one emitted tuple and its full compensated execution, including + coefficient preparation, rounding and delivery. The physical contract + and full caller invoice are missing. AAAmirinu's + [PR6239](https://github.com/CrocSwap/integer-mult-bounds/pull/6239) + already cautions that a smaller coefficient row alone need not reduce + child precision. + +This contribution was produced by Nick Petalas. Credit Jose Manuel / +jmbarrancoml and Vlado Vojdanovski for the residual and dirty-shear compiler, +Cinna Davis and the retained Gaussian/phase/caller authors, and classical +Gram-Schmidt and four-square theory. Apache-2.0. diff --git a/research/exact-unit-stein-splitting/SOURCES.json b/research/exact-unit-stein-splitting/SOURCES.json new file mode 100644 index 0000000000..32f4851378 --- /dev/null +++ b/research/exact-unit-stein-splitting/SOURCES.json @@ -0,0 +1,45 @@ +{ + "schema_version": 1, + "base_commit": "e54523517d0f03e58352e647c403a5af3e525a57", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "75d77b328ea02ded535e7f23aee5b1e55d0bb1a8", + "sha256": "154c3d30455d778fd3343b6b4851f0fa4665f0f7625b0adb3675e2bafbbf6bf1" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "3404fbdd5e7529a0674eef26114b5071e4586cec", + "sha256": "ead6cdc9d0af6d58e8a49744d8ec8c7c15e13d2ba93c276c434bcf9b51b93a81" + }, + { + "path": "research/rational-stein-generators/PROOF.md", + "usage": "provenance", + "git_blob": "075e828b68dca080ab20461fba57047e46b3b271", + "sha256": "90e5681d0e32be1926df7f3fd2715aa43f983fa7e4c01ee42a15322f9d19e639" + } + ], + "references": [ + { + "pr": 4198, + "author": "jmbarrancoml", + "role": "Main predecessor: exact orthogonal residual and Stein/dirty-shear interface." + }, + { + "pr": 6228, + "author": "petalas", + "role": "Background only: near-unit fixed-term normalization; no construction or executable reused." + }, + { + "pr": 6239, + "author": "AAAmirinu", + "role": "Background only: coefficient norm need not reduce rounded-child precision." + }, + { + "url": "https://leanprover-community.github.io/mathlib4_docs/Mathlib/NumberTheory/SumFourSquares.html", + "role": "Classical Lagrange four-square theorem; no Lean replay claimed." + } + ] +} diff --git a/research/exact-unit-stein-splitting/TASK.json b/research/exact-unit-stein-splitting/TASK.json new file mode 100644 index 0000000000..ff83aed511 --- /dev/null +++ b/research/exact-unit-stein-splitting/TASK.json @@ -0,0 +1,96 @@ +{ + "schema_version": 1, + "task_id": "uniform-transfer-audit", + "priority_revision": 8, + "base_commit": "e54523517d0f03e58352e647c403a5af3e525a57", + "original_card": { + "id": "uniform-transfer-audit", + "question": "Can one of the selected record's five inherited assumptions be proved at its actual caller?", + "starting_points": [ + "docs/research/current-priorities.md", + "docs/research/research-family-status.md", + "certificates/selected-result.json", + "docs/research/pr-intake-evidence/pr5599-completion/PROOF-REVIEW.md" + ], + "methods": [ + "Construct and independently check one source-bound instance of the stated first experiment.", + "Seek an exact counterexample or target-specific complete-budget obstruction to that same proposed change." + ], + "acceptance_gate": "Identify the proposition and what would change downstream; do not equate a finite replay with a universal proof.", + "first_step": { + "command": [ + "python3", + "-B", + "scripts/research_selected_baseline.py" + ], + "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY", + "inspect": "Read revision 8, the selected review and live exact-source handoffs. Then follow this card's bounded first_experiment.", + "scope": "Source-bound orientation only, not a supplier replay or theorem proof." + }, + "first_experiment": "Choose one assumption ID (native_all_columns, owned_erasure, compiler, precision_tape, analytic_recovery), one quantified implication and exact dependency closure. Check current handoffs, including #5625 row reentry and #5707 root-prefix lowering, without treating their partial scope as the whole theorem. Prove or falsify that implication with a caller-bound argument; stop at the chosen obligation. Finite tests or root digits do not discharge an all-size statement.", + "bottleneck_link": "Judge the actual current supplier contract and complete signed target/paid bill; frozen-profile ceilings require an applicability proof.", + "allowed_variations": [ + "One exact source head, helper/cut region or quantified proposition and changed hypothesis; check live handoffs first. Different seeds are not different scopes." + ], + "scope_change_triggers": [ + "Changing source head, helper/cut domain, operation alphabet or caller boundary without recording the new scope.", + "Treating an archived proposal or local numerical pass as complete supplier admission." + ], + "overlap_terms": [ + "native_all_columns", + "owned_erasure", + "compiler", + "precision_tape", + "analytic_recovery" + ], + "existing_checks": [ + { + "scope": "PR #4100: proposed unconditional construction and finite checks, without independent acceptance.", + "evidence": "research/unconditional-batched-crt/README.md", + "remaining": "Independent proof audit of one caller-bound all-size implication; all other premises remain explicit.", + "status": "author-checked", + "source_commit": "337d04939571a637aeab54a0f58f54fde811663a", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "Selected PR5599: complete helper-subspace construction and independent caller replay; five framework assumptions remain.", + "evidence": "docs/research/pr-intake-evidence/pr5599-completion/PROOF-REVIEW.md", + "remaining": "The next experiment must identify uncovered source-specific work. Successor archive intake does not transfer this admission.", + "status": "reviewed", + "source_commit": "09dd0965711ca7eeb169057a420abab1c9052137", + "selected_record_sha256": "2a0220e866417774606e0f991e41b9c1d26b1f43b97857238cf478dbc7d330b5" + } + ] + }, + "scope_status": "changed", + "actual_scope": "Exact unit-norm Gaussian-rational Stein generators with at most twice the rank-many terms.", + "scope_change_reason": "The donor requests an original result, not an audit. This changes the coefficient term-count contract at PR4198 rather than reviewing a selected framework assumption.", + "bottleneck_link": "Side result on the coefficient preparation boundary. Exact smaller norms could affect a rounded child, but extra generator macros and physical costs must be paid first.", + "budget_minutes": 60, + "stop_rule": "Ship the first exact factor-one construction and its bounded verification; do not expand into rounded/native implementation.", + "outcome": "completed", + "evidence": "PROOF.md proves the general construction; verify.py checks exact rational finite fixtures and emits a source-bound receipt.", + "author": "Nick Petalas", + "ai_assistance": "Independent human and formal proof-kernel review are unclaimed.", + "open_questions": [ + { + "question": "Can a universal rational Parseval construction use fewer than 2r terms?", + "why_it_matters": "Fewer added compensated generator macros could reduce overhead.", + "first_step": "Study a rank-two rational Hermitian form embedded into three standard coordinates.", + "blocked_by": "A rational embedding theorem and exact coefficient preparation algorithm are not supplied here." + }, + { + "question": "Does an actual rounded child benefit after all extra terms are charged?", + "why_it_matters": "Coefficient norm alone does not determine physical precision or cost.", + "first_step": "Bind one emitted tuple to the complete compensated rounded child and charge preparation, delivery and every extra call.", + "blocked_by": "Physical coefficient delivery and complete caller costs remain missing." + } + ], + "environment": { + "implementation": "CPython", + "version": "3.12.15", + "system": "Darwin", + "machine": "arm64", + "dependencies": "standard library only" + } +} diff --git a/research/exact-unit-stein-splitting/VERIFY.json b/research/exact-unit-stein-splitting/VERIFY.json new file mode 100644 index 0000000000..975d9ae294 --- /dev/null +++ b/research/exact-unit-stein-splitting/VERIFY.json @@ -0,0 +1,34 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [], + "commands": [ + { + "argv": [ + "{python}", + "-I", + "-W", + "error", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 120, + "checks": "Complete exact reconstruction, Parseval projector, independent rank and Hermitian PSD checks, square enumeration and rejection controls." + }, + { + "argv": [ + "{python}", + "-I", + "-O", + "-W", + "error", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 120, + "checks": "Same full contract with assertions disabled." + } + ] +} diff --git a/research/exact-unit-stein-splitting/algebra.py b/research/exact-unit-stein-splitting/algebra.py new file mode 100644 index 0000000000..94be17e65a --- /dev/null +++ b/research/exact-unit-stein-splitting/algebra.py @@ -0,0 +1,158 @@ +"""Bounded exact Gaussian-rational arithmetic for the coefficient construction.""" +from fractions import Fraction as Q +from math import gcd, isqrt, lcm + +C = tuple[Q, Q] +V = tuple[C, ...] +M = tuple[V, ...] +ZERO: C = (Q(0), Q(0)) +ONE: C = (Q(1), Q(0)) +MAX_DIM = 8 +MAX_HEIGHT = 256 +MAX_SQUARE_SUM = 65536 + + +def add(a: C, b: C) -> C: + return a[0] + b[0], a[1] + b[1] + + +def neg(a: C) -> C: + return -a[0], -a[1] + + +def mul(a: C, b: C) -> C: + return a[0]*b[0] - a[1]*b[1], a[0]*b[1] + a[1]*b[0] + + +def conj(a: C) -> C: + return a[0], -a[1] + + +def scale(a: C, q: Q) -> C: + return a[0]*q, a[1]*q + + +def total(values) -> C: + result = ZERO + for value in values: + result = add(result, value) + return result + + +def dot(a: V, b: V) -> C: + assert len(a) == len(b) <= MAX_DIM + return total(mul(conj(x), y) for x, y in zip(a, b)) + + +def norm2(v: V) -> Q: + value = dot(v, v) + assert value[1] == 0 and value[0] >= 0 + return value[0] + + +def adj(a: M) -> M: + return tuple(tuple(conj(x) for x in col) for col in zip(*a)) + + +def product(a: M, b: M) -> M: + assert a and b and len(a[0]) == len(b) + assert max(len(a), len(b), len(b[0])) <= 2*MAX_DIM + return tuple(tuple(total(mul(x, y) for x, y in zip(row, col)) + for col in zip(*b)) for row in a) + + +def outer(v: V) -> M: + return tuple(tuple(mul(x, conj(y)) for y in v) for x in v) + + +def matrix_sum(a: M, b: M) -> M: + assert len(a) == len(b) and len(a[0]) == len(b[0]) + return tuple(tuple(add(x, y) for x, y in zip(r, s)) for r, s in zip(a, b)) + + +def zero_matrix(rows: int, cols: int) -> M: + assert 1 <= rows <= 2*MAX_DIM and 1 <= cols <= 2*MAX_DIM + return tuple(tuple(ZERO for _ in range(cols)) for _ in range(rows)) + + +def validate(a: M) -> None: + if not isinstance(a, tuple) or not 1 <= len(a) <= MAX_DIM: + raise ValueError("matrix row capacity") + if not isinstance(a[0], tuple) or not 1 <= len(a[0]) <= MAX_DIM: + raise ValueError("matrix column capacity") + for row in a: + if not isinstance(row, tuple) or len(row) != len(a[0]): + raise ValueError("ragged matrix") + for value in row: + if not isinstance(value, tuple) or len(value) != 2: + raise ValueError("Gaussian-rational entry required") + for part in value: + if type(part) is not Q: + raise ValueError("exact Fraction components required") + if max(abs(part.numerator).bit_length(), part.denominator.bit_length()) > MAX_HEIGHT: + raise ValueError("input height capacity") + + +def square_pairs(n: int) -> tuple[int, int, int, int]: + """Two-square representations are preferred, then a bounded four-square search.""" + if type(n) is not int or not 1 <= n <= MAX_SQUARE_SUM: + raise ValueError("primitive squared-norm capacity") + pairs = [None]*(n+1) + for x in range(isqrt(n)+1): + for y in range(isqrt(n-x*x)+1): + pairs[x*x+y*y] = (x, y) + if pairs[n] is not None: + x, y = pairs[n] + return x, y, 0, 0 + for k in range(n+1): + if pairs[k] is not None and pairs[n-k] is not None: + x, y = pairs[k] + z, t = pairs[n-k] + assert x*x+y*y+z*z+t*t == n + return x, y, z, t + raise ArithmeticError("four-square exhaustion contradicts Lagrange") + + +def split_norm(v: V) -> tuple[C, C]: + """Represent ||v|| squared as two Gaussian-rational squared norms.""" + assert 1 <= len(v) <= MAX_DIM and norm2(v) > 0 + denominator = lcm(*(q.denominator for x in v for q in x)) + integers = tuple(int(q*denominator) for x in v for q in x) + divisor = gcd(*integers) + primitive = tuple(x//divisor for x in integers) + n = sum(x*x for x in primitive) + a, b, c, d = square_pairs(n) + factor = Q(divisor, denominator) + result = ((factor*a, factor*b), (factor*c, factor*d)) + assert sum(norm2((z,)) for z in result) == norm2(v) + return result + + +def residuals(a: M) -> tuple[V, ...]: + validate(a) + accepted = [] + for column in zip(*a): + v = tuple(column) + for previous in accepted: + coefficient = scale(dot(previous, v), 1/norm2(previous)) + v = tuple(add(x, neg(mul(coefficient, y))) for x, y in zip(v, previous)) + if norm2(v): + assert all(dot(old, v) == ZERO for old in accepted) + accepted.append(v) + return tuple(accepted) + + +def construct(a: M) -> tuple[tuple[V, ...], tuple[V, ...], tuple[V, ...]]: + basis = residuals(a) + columns = [] + rows = [] + for v in basis: + for z in split_norm(v): + if z != ZERO: + g = tuple(scale(mul(x, z), 1/norm2(v)) for x in v) + row = tuple(dot(g, tuple(col)) for col in zip(*a)) + assert 0 < norm2(g) <= 1 + columns.append(g) + rows.append(row) + assert len(basis) <= len(columns) <= 2*len(basis) + return tuple(columns), tuple(rows), basis diff --git a/research/exact-unit-stein-splitting/verify.py b/research/exact-unit-stein-splitting/verify.py new file mode 100644 index 0000000000..da9bed2c35 --- /dev/null +++ b/research/exact-unit-stein-splitting/verify.py @@ -0,0 +1,207 @@ +"""Check the new coefficient contract; never execute predecessor code.""" +import argparse +from fractions import Fraction as Q +from hashlib import sha256 +from itertools import product as cartesian +from math import isqrt +import json +from pathlib import Path +import sys + +PACKAGE = Path(__file__).resolve().parent +sys.path.insert(0, str(PACKAGE)) +from algebra import (C, V, M, ZERO, ONE, MAX_DIM, MAX_SQUARE_SUM, add, adj, + conj, construct, dot, matrix_sum, mul, neg, norm2, + outer, product, scale, square_pairs, validate, zero_matrix) + + +def require(condition: bool, message: str) -> None: + if not condition: + raise ArithmeticError(message) + + +def rank(a: M) -> int: + """Independent row elimination, not the column residual constructor.""" + rows = [list(row) for row in a] + pivot_row = 0 + for j in range(len(a[0])): + pivot = next((i for i in range(pivot_row, len(rows)) if rows[i][j] != ZERO), None) + if pivot is None: + continue + rows[pivot_row], rows[pivot] = rows[pivot], rows[pivot_row] + value = rows[pivot_row][j] + inverse = scale(conj(value), 1/norm2((value,))) + rows[pivot_row] = [mul(x, inverse) for x in rows[pivot_row]] + for i in range(len(rows)): + if i != pivot_row: + factor = rows[i][j] + rows[i] = [add(x, neg(mul(factor, y))) + for x, y in zip(rows[i], rows[pivot_row])] + pivot_row += 1 + if pivot_row == len(rows): + break + return pivot_row + + +def positive_semidefinite(a: M) -> bool: + """Exact Hermitian Schur criterion, including singular zero pivots.""" + if a != adj(a): + return False + h = [list(row) for row in a] + for k in range(len(a)): + pivot = h[k][k] + if pivot[1] or pivot[0] < 0: + return False + if pivot == ZERO: + if any(h[k][j] != ZERO for j in range(k+1, len(a))): + return False + else: + for i in range(k+1, len(a)): + for j in range(k+1, len(a)): + correction = scale(mul(h[i][k], h[k][j]), 1/pivot[0]) + h[i][j] = add(h[i][j], neg(correction)) + return True + + +def contract(a: M, columns: tuple[V, ...], rows: tuple[V, ...]) -> dict: + validate(a) + r = rank(a) + require(r <= len(columns) <= 2*r and len(columns) == len(rows), "term count") + projection = zero_matrix(len(a), len(a)) + rebuilt = zero_matrix(len(a), len(a[0])) + gram = product(adj(a), a) + for g, row in zip(columns, rows): + require(len(g) == len(a) and len(row) == len(a[0]), "factor shape") + require(0 < norm2(g) <= 1, "column norm") + require(row == tuple(dot(g, tuple(c)) for c in zip(*a)), "adjoint row") + projection = matrix_sum(projection, outer(g)) + term = tuple(tuple(mul(x, y) for y in row) for x in g) + rebuilt = matrix_sum(rebuilt, term) + row_gram = outer(tuple(conj(x) for x in row)) + difference = matrix_sum(gram, tuple(tuple(neg(x) for x in rr) for rr in row_gram)) + require(positive_semidefinite(difference), "row contraction") + require(rebuilt == a, "reconstruction") + require(projection == adj(projection), "Hermitian projection") + require(product(projection, projection) == projection, "Parseval projection") + require(product(projection, a) == a and rank(projection) == r, "column-space projection") + return {"rank": r, "terms": len(rows), "split_surplus": len(rows)-r} + + +def real(value) -> C: + return Q(value), Q(0) + + +def fixtures() -> tuple[M, ...]: + result = [] + for alphabet in [(real(-1), ZERO, ONE), (ZERO, ONE, (Q(0), Q(1)), (Q(1), Q(1)))]: + for values in cartesian(alphabet, repeat=4): + result.append((tuple(values[:2]), tuple(values[2:]))) + for m in range(1, 6): + for n in range(1, 6): + for denominator in (1, 3): + result.append(tuple(tuple(real(Q((i+1)*(j%2+1), denominator)) + for j in range(n)) for i in range(m))) + basis = ((real(1), real(1), real(1)), + (real(1), real(-1), ZERO), + (real(1), real(1), real(-2))) + for exponent in (0, 3, 120): + for shift in range(8): + a = tuple(tuple(add(scale(basis[j][i], Q(1, 2**exponent) if j == 2 else Q(1)), + scale(basis[0][i], Q(shift) if j else Q(0))) + for j in range(3)) for i in range(3)) + result.append(a) + result.extend((zero_matrix(1, 1), zero_matrix(3, 5), + ((ONE, ZERO), (ONE, ZERO), (ONE, ZERO)))) + return tuple(result) + + +def controls() -> int: + caught = 0 + invalid = [(), ((ONE,), (ONE, ZERO)), ((1.0,),), + tuple((ONE,) for _ in range(MAX_DIM+1)), + (((Q(1, 2**300), Q(0)),),)] + for matrix in invalid: + try: + construct(matrix) + except ValueError: + caught += 1 + else: + raise ArithmeticError("invalid matrix accepted") + invalid.append(((ONE,), ((Q(256), Q(0)),))) + try: + construct(invalid[-1]) + except ValueError: + caught += 1 + else: + raise ArithmeticError("residual capacity accepted") + for value in (0, -1, MAX_SQUARE_SUM+1, True): + try: + square_pairs(value) + except ValueError: + caught += 1 + else: + raise ArithmeticError("invalid square capacity accepted") + a = ((ONE,), (ONE,), (ONE,)) + columns, rows, _ = construct(a) + faults = [(columns[:1], rows[:1]), + (columns, tuple(tuple(conj(x) for x in row) for row in rows)), + (tuple(tuple(scale(x, Q(2)) for x in col) for col in columns), rows), + (columns, tuple(tuple(scale(x, Q(3)) for x in row) for row in rows))] + for cols, rr in faults: + try: + contract(a, cols, rr) + except ArithmeticError: + caught += 1 + else: + raise ArithmeticError("semantic fault accepted") + for bad in [((real(-1),),), ((ZERO, ONE), (ONE, ZERO)), + ((ONE, (Q(0), Q(1))), (ZERO, ONE))]: + require(not positive_semidefinite(bad), "PSD negative control") + caught += 1 + require(positive_semidefinite(((ZERO, ZERO), (ZERO, ONE))), "PSD singular control") + require(square_pairs(MAX_SQUARE_SUM) == (256, 0, 0, 0), "capacity endpoint") + return caught + + +def run() -> dict: + four_square_cases = 0 + for n in range(1, 2049): + a, b, c, d = square_pairs(n) + require(a*a+b*b+c*c+d*d == n, "four-square identity") + one_term_possible = any(x*x+y*y == n for x in range(isqrt(n)+1) + for y in range(isqrt(n)+1)) + require((c == d == 0) == one_term_possible, "two-square preference") + four_square_cases += 1 + stats = [contract(a, *construct(a)[:2]) for a in fixtures()] + example = ((ONE,), (ONE,), (ONE,)) + columns, rows, _ = construct(example) + require(len(columns) == 2, "rank-one exact-unit two-term witness") + require(sum(norm2(g) for g in columns) == 1, "rank-one Parseval trace") + result = {"status": "PASS_EXACT_UNIT_STEIN_SPLITTING", + "four_square_cases": four_square_cases, + "matrix_cases": len(stats), + "total_rank": sum(s["rank"] for s in stats), + "total_terms": sum(s["terms"] for s in stats), + "cases_using_extra_terms": sum(s["split_surplus"] > 0 for s in stats), + "rejection_controls": controls(), + "example_column_norms_squared": [str(norm2(g)) for g in columns], + "example_row_norms_squared": [str(norm2(row)) for row in rows]} + names = ("algebra.py", "verify.py", "README.md", "PROOF.md", "TASK.json", "SOURCES.json", "VERIFY.json") + result["source_sha256"] = {name: sha256((PACKAGE/name).read_bytes()).hexdigest() for name in names} + return result + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--receipt", type=Path) + args = parser.parse_args() + result = run() + rendered = json.dumps(result, sort_keys=True, indent=2)+"\n" + if args.receipt: + with args.receipt.open("x", encoding="utf-8") as output: + require(output.write(rendered) == len(rendered), "receipt write") + print(rendered, end="") + + +if __name__ == "__main__": + main()