From 0dfca71b00de4bc059f7ce3e9f320f741f7dbc01 Mon Sep 17 00:00:00 2001 From: Nick Petalas Date: Tue, 3 Nov 2026 00:47:52 +0000 Subject: [PATCH] feat(research): reduce bulk precision for cyclic-band solves --- research/row-adaptive-band-precision/PROOF.md | 201 ++++++++++++ .../row-adaptive-band-precision/README.md | 110 +++++++ .../row-adaptive-band-precision/SOURCES.json | 24 ++ .../row-adaptive-band-precision/TASK.json | 89 ++++++ .../row-adaptive-band-precision/VERIFY.json | 34 ++ .../row-adaptive-band-precision/verify.py | 294 ++++++++++++++++++ 6 files changed, 752 insertions(+) create mode 100644 research/row-adaptive-band-precision/PROOF.md create mode 100644 research/row-adaptive-band-precision/README.md create mode 100644 research/row-adaptive-band-precision/SOURCES.json create mode 100644 research/row-adaptive-band-precision/TASK.json create mode 100644 research/row-adaptive-band-precision/VERIFY.json create mode 100644 research/row-adaptive-band-precision/verify.py diff --git a/research/row-adaptive-band-precision/PROOF.md b/research/row-adaptive-band-precision/PROOF.md new file mode 100644 index 0000000000..08e2fe3605 --- /dev/null +++ b/research/row-adaptive-band-precision/PROOF.md @@ -0,0 +1,201 @@ +# Row-dependent precision in cyclic-band elimination + +This contribution was produced by Nick Petalas. Apache-2.0. + +## Claim and predecessor boundary + +This is a scoped written mathematical proof, with finite exact corroboration. +It changes the precision allocation in Ali Romman's +[PR5284](https://github.com/CrocSwap/integer-mult-bounds/pull/5284), +[cyclic-band proof](../cyclic-band-gap-solver/PROOF.md). That predecessor supplies +one-round fixed-integer elimination, the cyclic fill pattern, variable-gap +stability and its common-grid solve. Its row residual uses `N K h` in every +row. The new step permits monotone, unequal row grids and replaces that +residual by a row-support bound. The arithmetic method and classical +row-dominance estimates retain their predecessor credit. + +Let A be a real N by N dyadic matrix on the p-bit fractional grid. Suppose + + a_ii - sum_(j != i) |a_ij| >= delta > 0, + max_i sum_j |a_ij| <= M, + a_ij = 0 when min(|i-j|, N-|i-j|) > W, 0 <= W < N. + +No symmetry or sign condition on off-diagonal entries is required. Let b be +real or complex on a supplied dyadic grid. In the complex case use the norm +`max_i (|Re b_i| + |Im b_i|)` and assume it is at most Y>0. Choose a desired +absolute solution error t>0. Input formation and final-grid errors are outside +this contract. + +Take a symmetric symbolic elimination pattern containing the cyclic band, +closed under natural-order Schur updates. Let m_i count the strictly lower +factor positions and the diagonal/upper positions in row i, including positions +whose numerical value is zero. Write + + K = max(1,2M), + p <= q_0 <= ... <= q_(N-1), S_i = 2^q_i, + h_i = 1/S_i, e = max_i m_i K h_i. + +The following construction gives factors B=LU with `||B-A||_infinity <= e`. +If e<=delta/2, all its pivots are positive and + + ||U|| <= 3M/2, gap(U) >= delta/2, + ||L|| <= 2M/delta, ||B^-1|| <= 2/delta, + ||L^-1|| <= 3M/delta. + +Every norm on a real matrix is its absolute row-sum norm. + +## Exact mixed-grid construction + +Store unfinished Schur row i at scale S_i squared. Initialize +`C_ij = a_ij S_i^2`, exactly. At pivot k perform + + U_int[k,j] = trunc(C[k,j]/S_k), j>=k, + L_int[i,k] = trunc(C[i,k]/(U_int[k,k] 2^(q_i-q_k))), i>k, + C[i,j] -= L_int[i,k] U_int[k,j] 2^(q_i-q_k), i,j>k. + +Visit the declared symbolic positions only. Here trunc rounds toward zero, +`U[k,j]=U_int[k,j]/S_k`, `L[i,k]=L_int[i,k]/S_i`, and L has unit diagonal. +The nonnegative exponent `q_i-q_k` is essential. It lets each update be an +exact integer subtraction on the row-i grid without increasing that grid. +It also makes the lower formula precisely rounding `C_real[i,k]/U[k,k]` +to row i's factor grid. No multiplication or Schur update is rounded. + +Expanding row i's exact updates shows that a completed upper position has +residual less than h_i, and a completed lower position has residual less than +`h_i |U[j,j]|`. The same identities apply to a partial factorization: put +these residuals at completed positions of A and zero at unfinished positions. +Its completed factors are exact factors of this partially perturbed matrix. +Outside the symbolic factor pattern the residual is zero, since the pattern +is closed under elimination. Consequently a partial row residual is at most +`m_i K h_i` whenever completed pivots have magnitude at most 2M. + +For completeness, the needed classical exact row-dominance facts are as +follows. Eliminating multiplier l from a row of a positive-diagonal matrix +with row gap g and row norm R increases the surviving row gap by at least +`|l| g` and decreases its absolute row sum by at least `|l| g`. Pivots stay +positive. If s is the accumulated absolute multiplier sum in one row, its +final gap is at least `g(1+s)` while its final row norm is at most `R-gs`. +Thus `||L|| <= (R+g)/(2g)`, `||U||<=R`, `gap(U)>=g`, and +`||B^-1||<=1/g`. These are the estimates used in PR5284, cited here as +predecessor mathematics rather than an audit target. + +Apply these facts to each partially perturbed matrix with +`R=M+e`, `g=delta-e`. Since e<=delta/2<=M/2, a newly completed pivot has +magnitude at most `M+e<=3M/2<2M` and is at least delta/2. This closes the +induction from the first pivot, proves every division is defined and proves +the factor residual claim. The L bound follows from +`(M+delta)/(2(delta-e))<=2M/delta`. Finally +`L^-1=U B^-1` gives `||L^-1||<=3M/delta`. The same arguments hold for every +leading principal factor pair, since deleting columns cannot decrease a +row-dominance margin. + +## An independently chosen solve grid + +Choose a common solve grid step `h_s=2^(-q_s)` containing b's input grid. +Accumulate each triangular dot product exactly, then round the solved +coordinate once. A lower or upper row-i product has grid `S_i 2^q_s`, so the +mixed factor grids introduce no extra numerical rounding. Execute separately +on the real and imaginary components. The triangular equations obey + + L z = b + f, ||f|| <= 2 h_s, + U x = z + g, ||g|| <= 3 M h_s. + +Hence `B x=b+r`, with + + ||r|| <= R_s = 2 h_s (1+3M^2/delta). + +The factor two is a sufficient common bound for the complex component norm; +a real solve is covered too. The resolvent identity and `||A^-1||<=1/delta` +give + + ||x-A^-1 b|| <= (eY + delta R_s)/(delta(delta-e)). + +In particular it suffices to choose + + e <= min(delta/2, t delta^2/(4Y)), + h_s <= t delta^2/(8delta+24M^2). + +The first contribution is at most `2eY/delta^2<=t/2`, and the second is at +most `2R_s/delta<=t/2`. Factor precision and solution precision can therefore +be chosen separately. If b has a finer input grid, increase q_s to include it. + +## The cyclic-band improvement + +For N>2W put the last W indices in the border. Natural-order elimination of +a bulk index connects only its next W bulk indices and the W border indices. +A bulk row has at most W earlier bulk neighbors and at most 2W later +neighbors. Therefore + + m_i <= min(N,3W+1) for i=N-W. + +This is also the usual fill argument in PR5284, now used to price residuals +row by row. If N<=2W, use the dense count N everywhere. The formulas above +then still give the same safe allocation because `min(N,3W+1)=N`. If W=0, +every row has one position and no border is needed. + +For a funded e-budget define + + q_bulk = max(p, ceil(log2(min(N,3W+1) K / e_budget))), + q_high = max(p, ceil(log2(N K / e_budget))). + +Assign q_bulk to the first N-W rows and q_high to the final W rows. These +precisions are nondecreasing, so the mixed-grid theorem applies. The full +length enters only the border factor precision. For fixed W, delta, M, Y, +p and t, the required precision on every bulk factor row is independent of N. +This removes a logarithmic length term from those rows' sufficient guards. +The common-grid allocation using the same e-budget requires q_high on every +row. No assertion is made that either sufficient allocation is optimal. + +The arithmetic schedule still has at most `4 N W^2` factor products and +`8 N W` products per complex solve. There are at most `(4W+1)N` stored factor +positions. The new shifts are exact rescalings, not additional oracle calls. +This is an arithmetic/precision statement; it does not certify a native +instruction tariff, variable-record tape layout, wall-time gain or complete +recursive bill. A machine choosing one common word width still pays its +maximum required width. + +## Integer ranges + +One can retain the predecessor's dimension-free magnitude bounds. Put + + L_0=2M/delta, U_0=3M/2, + Z=(3M/delta)(Y+2h_s), X=(2/delta)(Y+R_s), + C=max(1,M+L_0 U_0,Y+L_0 Z,Z+U_0 X,L_0,U_0,Y,Z,X). + +Absolute partial dot sums are bounded by their row sums. Thus factor row i's +Schur entries, update products and division numerators need magnitude at +most `C S_i^2`. Its divisor `U_int[k,k] 2^(q_i-q_k)` has magnitude at most +`U_0 S_i`. A signed width `2q_i+ceil(log2 C)+2` suffices for factor work. +Solve row i's accumulators need at most `C S_i 2^q_s`, so width +`q_i+q_s+ceil(log2 C)+2` suffices there. Forward prefixes use the leading +principal estimates above. Backward intermediates are already coordinates +of the final x, so the final bound covers them. Finite Python execution uses +unbounded integers; these are written sufficient range bounds, not an emitted +bounded-word implementation certificate. + +## First strict fixture and finite evidence + +Take N=65, W=1, diagonal 1, and cyclic adjacent entries -63/128. Its row gap +is delta=1/64. Use M=2, Y=1, alternating right-hand side 1,-1, and target +`t=2^-16`. The chosen e-budget is `2^-30`. The allocation is + + factor rows 0..63: q=34, factor row 64: q=39, + common solve grid: q_s=35, matched common factor grid: q=39. + +There are 319 symbolic factor positions, 254 in bulk rows. Their declared +fraction-bit total falls from 12,441 to 11,171, a saving of 1,270. This counts +fractional precision per symbolic factor position, including zeros; it is +not a measurement of Python memory or a physical encoded representation. +The constructor performs 253 factor update products. The exact original- +matrix residual is `1443/4398046511104`, and the general sufficient error +bound is `1793/268435440 < 2^-16`. Independently, residual divided by delta +also certifies the requested solution accuracy. + +The verifier additionally checks the support bounds on 1,125 symbolic patterns +with 1,171,245 row positions, and checks 252 signed and unsigned, generally +nonsymmetric small dyadic systems against an independent pivoted rational +solver. These are finite corroboration of the written all-size result, not a +replacement for it. No predecessor scientific suite is replayed. Complex +vectors are covered by the written two-real-solve argument; the executable +fixtures use real right-hand sides. diff --git a/research/row-adaptive-band-precision/README.md b/research/row-adaptive-band-precision/README.md new file mode 100644 index 0000000000..0987004eeb --- /dev/null +++ b/research/row-adaptive-band-precision/README.md @@ -0,0 +1,110 @@ +# Smaller bulk guards for cyclic-band elimination + +This contribution was produced by Nick Petalas. Apache-2.0. + +Claim type: scoped written mathematical proof with exact finite corroboration. +For a real dyadic cyclic-band matrix of half-width W and row-dominance gap +delta>0, natural-order elimination can use lower factor precision on its first +N-W rows. Its final W border rows retain the full-length guard: + +```text +predecessor factor residual: N K 2^(-q) in every row +new factor residual: max_i m_i K 2^(-q_i) +bulk support count: m_i <= min(N,3W+1) +border support count: m_i <= N +``` + +[PROOF.md](PROOF.md) proves the mixed-grid update identity, stability induction, +separate solve precision, integer ranges and all-size cyclic support bound. +The precisions q_i must be nondecreasing. This is an original precision +allocation step on Ali Romman's [PR5284](https://github.com/CrocSwap/integer-mult-bounds/pull/5284) +algorithm, not an audit of that result. + +On the first strict 65-row tridiagonal fixture, bulk factor precision is 34 +fractional bits instead of the matched uniform allocation's 39; the final row +uses 39. The abstract factor fraction-bit total is 11,171 instead of 12,441. +Exact residuals certify the same 2^-16 solution target. No optimality, +wall-time, native/tape bill, Gaussian-family admission, new kappa or selected- +frontier change is claimed. The five selected framework assumptions remain. + +## Reproduce + +From the repository root, using CPython 3.12 and its standard library: + +```sh +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -B scripts/research_package.py --package research/row-adaptive-band-precision +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -B scripts/check_research_commit.py --package research/row-adaptive-band-precision --check-source-history +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -W error -B research/row-adaptive-band-precision/verify.py +``` + +The absolute interpreter path records the local environment. A different +CPython 3.12 executable can replace it. VERIFY.json registers normal and +optimized warning-strict runs. Each must print +`PASS_ROW_ADAPTIVE_BAND_PRECISION`, `errors=0`, 252 small exact solve comparisons, +1,125 symbolic patterns, and the strict 65-row precision reduction. JSON stdout +also binds all package-file hashes. Receipts are generated, not handwritten. +No output files or bytecode caches are created by the verifier. + +The all-size theorem is written mathematics. Finite fixtures are real, +N<=129 for pattern checks and N<=65 for numerical checks; the implementation's +explicit capacity is N<=257, W<=8, q<=512. Small numerical systems are checked +against an independently transcribed pivoted rational solver. This does not +claim independent human or formal proof-kernel review. The complex extension +is a written two-real-solve argument. No full predecessor, selected supplier, +Lean, Docker, native backend or physical-tape replay is required or claimed. + +## Scope and credit + +Ali Romman / PR5284 supplies the variable-gap solver, one-round factor-entry +method, diagonal-dominance stability estimates and cyclic fill. The new work +is the monotone row-grid theorem, its support-priced residual and reduced bulk +allocation. Dugongue / PR4554 supplies the earlier band/Woodbury interface; +Ryanson Jonathan / PR4556 supplies its independent check, cited as context and +not replayed. PR5260 and PR5326 are prior fixed-source and Gaussian Reference +consumers. Classical Gaussian elimination and row dominance retain their +classical status. No worldwide novelty claim is made for those methods. +Frozen predecessor files and their credit lines are unchanged. + +The assigned uniform-transfer-audit card is explicitly changed to this +numerical-interface result. It addresses precision costs that a future paid +Gaussian caller must fund; it does not discharge the selected precision_tape +assumption. SOURCES.json pins the exact predecessor proof and originating +catalogue. The verifier imports no predecessor code. + +The advance [CLAIM](https://github.com/CrocSwap/integer-mult-bounds/discussions/1166#discussioncomment-19363375) +reserved this question. Live overlap covered 430 open PR bodies and 6,404 +Discussion records with complete remote pagination at the recorded snapshot. +The initial sparse local search covered 115 of 8,917 tracked notes. Relevant +proofs were then read directly, and main/merged searches covered mixed row +precision and cyclic-band terms. Arbitrary diffs, unmaterialized notes and +unpublished work remain outside that coverage. No novelty certificate is +claimed. + +## Open next questions + +Can an actual Gaussian production caller use the row allocation to reduce its +paid precision/storage bill? First freeze its coefficient error, row gap, +symbolic fill and record representation, then account for the W high-precision +border rows and solve buffers. Source generation, fixed-tape layout, native +operations and the complete recursive invoice remain blockers. A single common +machine-word representation may obtain no storage benefit. + +Can a different elimination order lower the number of high-precision rows +without increasing fill or the required row-residual budget? First derive its +symbolic row counts and a valid mixed-grid schedule. The monotone grid rule in +this result may prevent a naive permutation from applying directly. No such +optimization is included here. + +## Verifier constraints + +Every loop is bounded by explicit matrix, bandwidth, precision or fixture +limits. Functions are under 60 lines. Invalid shapes, grids and divisors fail +with explicit errors retained under optimized Python. The verifier uses no +recursion, dynamic code generation, subprocesses or network calls. Python's +managed allocation is a Power of Ten rule 3 deviation, bounded here by N<=257, +W<=8, q<=512, finite fixture counts and a registered 60-second timeout. +Package hashing is bounded to eight files of at most one MiB each. There is no +claim of a fixed physical Python heap. Reassess storage if those limits change. +Warnings-as-errors in both interpreter modes are the available diagnostics; +no additional static analyzer is configured for this isolated standard-library +research verifier. diff --git a/research/row-adaptive-band-precision/SOURCES.json b/research/row-adaptive-band-precision/SOURCES.json new file mode 100644 index 0000000000..f5ef5c00df --- /dev/null +++ b/research/row-adaptive-band-precision/SOURCES.json @@ -0,0 +1,24 @@ +{ + "schema_version": 1, + "base_commit": "7c49fe32c55de02c2c873655c804b16e4c3fef0a", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "d59b5035e7b4cbaa9d688a1d32dd0d6eb79c466a", + "sha256": "9b8d0eea95b3bf0dfcb7938a732a6439b87a2685dec2a5ab9e525b3b34debbfa" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "3404fbdd5e7529a0674eef26114b5071e4586cec", + "sha256": "ead6cdc9d0af6d58e8a49744d8ec8c7c15e13d2ba93c276c434bcf9b51b93a81" + }, + { + "path": "research/cyclic-band-gap-solver/PROOF.md", + "usage": "input", + "git_blob": "924929376c304d743a15e67aed2ac7712fd89892", + "sha256": "508fd707fcfe1490ee3de0047b7d4d699cf341d2d22dfce1ea297e16e9beebb9" + } + ] +} diff --git a/research/row-adaptive-band-precision/TASK.json b/research/row-adaptive-band-precision/TASK.json new file mode 100644 index 0000000000..ea9169bdb2 --- /dev/null +++ b/research/row-adaptive-band-precision/TASK.json @@ -0,0 +1,89 @@ +{ + "schema_version": 1, + "task_id": "uniform-transfer-audit", + "priority_revision": 8, + "base_commit": "7c49fe32c55de02c2c873655c804b16e4c3fef0a", + "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": "Monotone row-dependent factor precision for variable-gap real dyadic cyclic-band elimination, with an independently chosen solve grid and a smaller bulk support residual.", + "scope_change_reason": "The audit card is changed to an original numerical-interface construction. No selected assumption is discharged; the result reduces sufficient factor guards in a reusable cyclic-band consumer.", + "bottleneck_link": "Precision and retained storage are costs a Gaussian caller must fund. This side lemma removes the full-length residual count on bulk factor rows, but supplies neither native/tape realization nor a favorable complete signed invoice.", + "budget_minutes": 60, + "stop_rule": "Ship the first exact strict factor-precision improvement. Do not optimize elimination order or build a physical Gaussian caller.", + "outcome": "completed", + "evidence": "PROOF.md gives the all-size row-grid theorem. verify.py independently constructs factors and checks exact residuals, 252 pivoted rational solve comparisons, 1125 symbolic patterns, a matched common-grid control and eight rejection controls.", + "author": "Nick Petalas", + "ai_assistance": "Independent human and formal proof-kernel review are unclaimed.", + "open_questions": [ + { + "question": "Can an actual Gaussian production caller lower its paid precision/storage bill using these row grids?", + "why_it_matters": "A numerical guard reduction only improves the complete supplier when its coefficient and record costs decrease.", + "first_step": "Freeze one caller row gap, symbolic fill, coefficient error and physical record layout, then charge border rows and solve buffers.", + "blocked_by": "Source generation; Native operators; Fixed-tape layout; Complete recursive invoice" + }, + { + "question": "Can an elimination order reduce high-precision rows without increasing fill or invalidating the grid schedule?", + "why_it_matters": "The final W border rows still pay the full-length precision count.", + "first_step": "Derive row supports for one alternative order and establish an exact mixed-grid update invariant.", + "blocked_by": "Monotone row-grid requirement; Fill growth and full residual accounting" + } + ] +} diff --git a/research/row-adaptive-band-precision/VERIFY.json b/research/row-adaptive-band-precision/VERIFY.json new file mode 100644 index 0000000000..9f42d91a2d --- /dev/null +++ b/research/row-adaptive-band-precision/VERIFY.json @@ -0,0 +1,34 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [ + "research/cyclic-band-gap-solver/PROOF.md" + ], + "commands": [ + { + "argv": [ + "{python}", + "-W", + "error", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 60, + "checks": "Exact mixed-grid LU residuals, row-support bounds, independent rational small solutions, first strict precision gain, matched uniform allocation and rejection controls." + }, + { + "argv": [ + "{python}", + "-O", + "-W", + "error", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 60, + "checks": "The identical finite certificate and controls with interpreter assertions disabled." + } + ] +} diff --git a/research/row-adaptive-band-precision/verify.py b/research/row-adaptive-band-precision/verify.py new file mode 100644 index 0000000000..f62ad74d57 --- /dev/null +++ b/research/row-adaptive-band-precision/verify.py @@ -0,0 +1,294 @@ +"""Exact finite evidence for row-dependent cyclic-band precision. + +This contribution was produced by Nick Petalas. Apache-2.0. +""" +from fractions import Fraction as F +from hashlib import sha256 +from itertools import islice +from pathlib import Path +import json + +MAX_N = 257 +MAX_W = 8 +MAX_Q = 512 + + +def need(condition, message): + if not condition: + raise ValueError(message) + + +def trunc_ratio(numerator, denominator): + need(denominator > 0, "nonpositive divisor") + magnitude = abs(numerator) // denominator + return -magnitude if numerator < 0 else magnitude + + +def precision(bound, minimum): + need(bound > 0 and 0 <= minimum <= MAX_Q, "precision arguments") + for q in range(minimum, MAX_Q + 1): + if F(1, 1 << q) <= bound: + return q + raise ValueError("precision capacity") + + +def symbolic(n, w): + need(1 <= n <= MAX_N and 0 <= w <= min(MAX_W, n - 1), "shape") + adjacency = [set() for _ in range(n)] + for i in range(n): + for j in range(i + 1, n): + if min(j - i, n - j + i) <= w: + adjacency[i].add(j) + adjacency[j].add(i) + lower = [[] for _ in range(n)] + upper = [] + products = 0 + for k in range(n): + later = sorted(j for j in adjacency[k] if j > k) + upper.append([k] + later) + for i in later: + lower[i].append(k) + for j in later: + adjacency[i].add(j) + products += 1 + counts = [len(lower[i]) + len(upper[i]) for i in range(n)] + need(products <= 4 * n * w * w, "symbolic product bound") + return lower, upper, counts, products + + +def validate(a, p, delta, maximum, rhs, y_bound): + n = len(a) + need(1 <= n <= MAX_N and len(rhs) == n, "matrix/RHS shape") + need(0 <= p <= 64 and delta > 0 and maximum >= delta, "parameters") + need(y_bound > 0 and max(map(abs, rhs)) <= y_bound, "RHS bound") + for i, row in enumerate(a): + need(len(row) == n, "matrix row") + need(all((v * (1 << p)).denominator == 1 for v in row), "input grid") + off = sum(abs(row[j]) for j in range(n) if j != i) + need(row[i] - off >= delta, "row gap") + need(sum(map(abs, row)) <= maximum, "row norm") + + +def allocate(n, w, p, delta, maximum, y_bound, target): + need(target > 0 and delta > 0 and y_bound > 0, "error budgets") + k_bound = max(F(1), 2 * maximum) + e_budget = min(delta / 2, target * delta * delta / (4 * y_bound)) + solve_step = target * delta * delta / (8 * delta + 24 * maximum * maximum) + q_solve = precision(solve_step, p) + q_global = precision(e_budget / (n * k_bound), p) + bulk = min(n, 3 * w + 1) + q_bulk = precision(e_budget / (bulk * k_bound), p) + grids = [q_bulk if i < n - w else q_global for i in range(n)] + need(all(grids[i] <= grids[i + 1] for i in range(n - 1)), "monotone grids") + return grids, q_solve, q_global, e_budget + + +def factor(a, lower, upper, grids): + n = len(a) + need(len(grids) == n and all(0 <= q <= MAX_Q for q in grids), "factor grids") + need(all(grids[i] <= grids[i + 1] for i in range(n - 1)), "grid order") + scales = [1 << q for q in grids] + c = [{} for _ in range(n)] + l = [{} for _ in range(n)] + u = [{} for _ in range(n)] + for i in range(n): + for j in lower[i] + upper[i]: + exact = a[i][j] * scales[i] * scales[i] + need(exact.denominator == 1, "Schur input grid") + c[i][j] = exact.numerator + operations = 0 + for k in range(n): + for j in upper[k]: + u[k][j] = trunc_ratio(c[k][j], scales[k]) + need(u[k][k] > 0, "positive pivot") + for i in upper[k][1:]: + shift = 1 << (grids[i] - grids[k]) + l[i][k] = trunc_ratio(c[i][k], u[k][k] * shift) + for j in upper[k][1:]: + need(j in c[i], "symbolic update closure") + c[i][j] -= l[i][k] * u[k][j] * shift + operations += 1 + return l, u, operations + + +def solve(l, u, grids, q_solve, rhs): + n = len(l) + s_solve = 1 << q_solve + need(all((v * s_solve).denominator == 1 for v in rhs), "RHS solve grid") + z = [0] * n + x = [0] * n + for i in range(n): + scale = 1 << grids[i] + total = (rhs[i] * scale * s_solve).numerator + total -= sum(value * z[j] for j, value in l[i].items()) + z[i] = trunc_ratio(total, scale) + for i in range(n - 1, -1, -1): + scale = 1 << grids[i] + total = z[i] * scale + total -= sum(value * x[j] for j, value in u[i].items() if j != i) + x[i] = trunc_ratio(total, u[i][i]) + return [F(value, s_solve) for value in x] + + +def factor_residual(a, l, u, grids): + n = len(a) + residual = [[-v for v in row] for row in a] + for i in range(n): + for k, value in {**l[i], i: 1 << grids[i]}.items(): + left = F(value, 1 << grids[i]) + for j, right in u[k].items(): + residual[i][j] += left * F(right, 1 << grids[k]) + return max(sum(map(abs, row)) for row in residual) + + +def exact_solve(a, rhs): + """Independent pivoted rational reference, used only for n <= 12.""" + n = len(a) + need(n <= 12, "reference capacity") + rows = [list(a[i]) + [rhs[i]] for i in range(n)] + for k in range(n): + pivot = k + for i in range(k + 1, n): + if abs(rows[i][k]) > abs(rows[pivot][k]): + pivot = i + need(rows[pivot][k] != 0, "reference singularity") + rows[k], rows[pivot] = rows[pivot], rows[k] + value = rows[k][k] + rows[k] = [v / value for v in rows[k]] + for i in range(n): + if i != k: + value = rows[i][k] + rows[i] = [rows[i][j] - value * rows[k][j] for j in range(n + 1)] + return [row[-1] for row in rows] + + +def fixture(n, w, gap, signed): + a = [[F(0) for _ in range(n)] for _ in range(n)] + for i in range(n): + for j in range(n): + if i != j and min(abs(i - j), n - abs(i - j)) <= w: + coefficient = ((17 * i + 11 * j) % 4 + 1) + sign = -1 if not signed or (i + 3 * j) % 3 else 1 + a[i][j] = F(sign * coefficient, 32) + a[i][i] = gap + sum(map(abs, a[i])) + rhs = [F((i % 5) - 2, 2) for i in range(n)] + return a, rhs + + +def run_case(a, rhs, w, p, delta, maximum, target): + validate(a, p, delta, maximum, rhs, F(1)) + lower, upper, counts, expected_ops = symbolic(len(a), w) + for i, row in enumerate(a): + need(all(v == 0 or j in lower[i] + upper[i] for j, v in enumerate(row)), "input pattern") + grids, q_solve, q_global, e_budget = allocate(len(a), w, p, delta, maximum, F(1), target) + l, u, products = factor(a, lower, upper, grids) + need(products == expected_ops, "complete factor work") + actual = factor_residual(a, l, u, grids) + bound = max(counts[i] * max(F(1), 2 * maximum) / (1 << grids[i]) for i in range(len(a))) + need(actual <= bound <= e_budget <= delta / 2, "factor certificate") + x = solve(l, u, grids, q_solve, rhs) + residual = max(abs(sum(a[i][j] * x[j] for j in range(len(a))) - rhs[i]) for i in range(len(a))) + h_solve = F(1, 1 << q_solve) + r_solve = 2 * h_solve * (1 + 3 * maximum * maximum / delta) + error_bound = (bound + delta * r_solve) / (delta * (delta - bound)) + need(error_bound <= target and residual <= delta * target, "solution certificate") + if len(a) <= 12: + reference = exact_solve(a, rhs) + need(max(abs(v - r) for v, r in zip(x, reference)) <= error_bound, "independent solution") + return {"n": len(a), "w": w, "bulk_q": grids[0], "border_q": grids[-1], + "solve_q": q_solve, "global_q": q_global, "entries": sum(counts), + "fraction_bits": sum(c * q for c, q in zip(counts, grids)), + "global_fraction_bits": sum(counts) * q_global, + "factor_products": products, "residual": str(residual), "error_bound": str(error_bound)} + + +def pattern_checks(): + cases = 0 + positions = 0 + for n in range(1, 130): + for w in range(min(MAX_W, n - 1) + 1): + lower, upper, counts, _ = symbolic(n, w) + for i, count in enumerate(counts): + need(count <= (min(n, 3 * w + 1) if i < n - w else n), "row support bound") + need(sum(counts) <= (4 * w + 1) * n, "storage bound") + positions += sum(counts) + cases += 1 + return cases, positions + + +def rejection_checks(): + count = 0 + for denominator in (0, -1): + try: + trunc_ratio(1, denominator) + except ValueError: + count += 1 + else: + raise ValueError("bad divisor accepted") + for shape in ((0, 0), (MAX_N + 1, 1), (4, 4), (20, MAX_W + 1)): + try: + symbolic(*shape) + except ValueError: + count += 1 + else: + raise ValueError("bad shape accepted") + a, rhs = fixture(4, 1, F(1, 8), False) + lower, upper, _, _ = symbolic(4, 1) + try: + factor(a, lower, upper, [20, 19, 20, 20]) + except ValueError: + count += 1 + else: + raise ValueError("descending grid accepted") + try: + validate(a, 2, F(1, 8), F(2), rhs, F(1)) + except ValueError: + count += 1 + else: + raise ValueError("off-grid input accepted") + return count + + +def main(): + pattern_cases, positions = pattern_checks() + small = 0 + for n in range(1, 13): + for w in range(min(3, n - 1) + 1): + for exponent in (1, 4, 10): + for signed in (False, True): + gap = F(1, 1 << exponent) + a, rhs = fixture(n, w, gap, signed) + maximum = max(sum(map(abs, row)) for row in a) + # The summary is discarded; run_case checks every invariant. + _ = run_case(a, rhs, w, 10, gap, maximum, F(1, 1 << 12)) + small += 1 + n = 65 + a = [[F(0) for _ in range(n)] for _ in range(n)] + for i in range(n): + a[i][i] = F(1) + a[i][(i - 1) % n] = a[i][(i + 1) % n] = F(-63, 128) + rhs = [F(1 if i % 2 == 0 else -1) for i in range(n)] + improvement = run_case(a, rhs, 1, 7, F(1, 64), F(2), F(1, 1 << 16)) + need(improvement["fraction_bits"] < improvement["global_fraction_bits"], "strict precision gain") + lower, upper, counts, expected_ops = symbolic(n, 1) + uniform = [improvement["global_q"]] * n + l, u, products = factor(a, lower, upper, uniform) + need(products == expected_ops == improvement["factor_products"], "matched common-grid work") + common_x = solve(l, u, uniform, improvement["solve_q"], rhs) + common_residual = max(abs(sum(a[i][j] * common_x[j] for j in range(n)) - rhs[i]) for i in range(n)) + need(common_residual <= F(1, 64 * (1 << 16)), "matched common-grid accuracy") + rejected = rejection_checks() + root = Path(__file__).resolve().parent + entries = list(islice(root.iterdir(), 9)) + need(len(entries) <= 8, "package entry capacity") + files = sorted(p for p in entries if p.is_file()) + need(len(files) <= 8 and all(p.stat().st_size <= 1048576 for p in files), "receipt capacity") + hashes = {p.name: sha256(p.read_bytes()).hexdigest() for p in files} + print(json.dumps({"status": "PASS_ROW_ADAPTIVE_BAND_PRECISION", "errors": 0, + "small_solves": small, "pattern_cases": pattern_cases, + "pattern_positions": positions, "rejections": rejected, + "improvement": improvement, "files": hashes}, sort_keys=True)) + + +if __name__ == "__main__": + main()