diff --git a/research/pr4554-lowering-check/.gitattributes b/research/pr4554-lowering-check/.gitattributes new file mode 100644 index 0000000000..0415142f74 --- /dev/null +++ b/research/pr4554-lowering-check/.gitattributes @@ -0,0 +1 @@ +PR4554-PHASE-TABLES.md -text diff --git a/research/pr4554-lowering-check/PR4554-PHASE-TABLES.md b/research/pr4554-lowering-check/PR4554-PHASE-TABLES.md new file mode 100644 index 0000000000..2dab8f3e6d --- /dev/null +++ b/research/pr4554-lowering-check/PR4554-PHASE-TABLES.md @@ -0,0 +1,221 @@ +# Direct bounded-word phase preparation and linear-program lowering + +Prepared for Dugongue with substantial OpenAI Codex assistance. Apache-2.0. +This packet contains constructive native arithmetic components. No Turing-machine +simulation, transition table, or Turing-machine multiplication theorem is used. +It does not by itself establish a new multiplication exponent. + +## 1. Paid arbitrary root powers without a per-query logarithmic factor + +**Theorem.** Let D>=2, R=2^t>=2, d=ceil(log_R D), E=R^d, and Q=2^p satisfy +Q>=64(E+d). Suppose a supplied signed Gaussian dyadic z0 has each-coordinate +error less than 1/Q from a unit complex number z. There is a direct bounded-word +algorithm constructing d tables of R Gaussian dyadics such that: + +* setup takes O(d(R+t)) word operations, plus the separately paid seed; +* every query 0<=e=2p+8 bits; tables require O(dR) words; +* the algorithm uses only bounded-word MUL, ADD/SUB, Boolean operations, + shifts, indirect memory accesses, and control. Neither an input-sized + multiplication nor a transcendental instruction is present. + +The unit complex number may be z=exp(2*pi*i/D). The theorem does not assume that +root acquisition is free. The recovery worker supplies `root_seed.py`, a native +Machin/Taylor construction with a proved coordinate error <1/Q and O(p) +bounded-word operations at a larger, still O(p+log D), physical word width. +The common machine must use the maximum width required by seed and tables. + +### Construction + +Define z_j by t repeated rounded squarings per level, starting with z_0=z0, +so it approximates z^(R^j). Store T_j[0]=1 and generate + + T_j[a] = rounded_product(T_j[a-1], z_j), 1<=a=64(E+d) gives +(2E+d)*eta <= 1/16. The same bound at all intermediate table/base nodes is +less than 1/4. Their exact targets have modulus 1, so all actual coordinates +have magnitude <2. The integer coordinate magnitude is therefore <=2Q. +Each individual product has magnitude <=4Q^2 and each real/imaginary sum at +most 8Q^2, which fits the stated signed-word envelope. Rounding is performed +with sign handling, addition of Q/2 and a right shift; it uses no unbounded +integer division. This proves the error, word and cost claims. ∎ + +### Uniform polynomial-order specialization + +Suppose D<=N^C for one fixed C. Choose + + t=max(1,floor(log2(N+2))/16), R=2^t. + +Then d=O_C(1), E<=RD=O(N^(C+1/16)), and table setup is +O_C(N^(1/16)+log N). To request error at most N^-A, it suffices to choose +p >= (A+C+1)*ceil(log2(N+2)) plus a fixed constant depending on C. +All coordinates and products therefore use a fixed number of Theta(log N)-bit +words; choosing physical word width large enough makes each primitive product +constant cost. Including the native seed costs O_C(N^(1/16)+polylog N), and +each requested root power costs O_C(1). + +The parameter R, each table exponent, and all table addresses fit the same +polynomial address envelope. A caller allocates disjoint table/data arenas and +includes all IO. No amortization over multiple multiplication inputs is used. + +## 2. A constructive sufficient finite-precision lowering theorem + +A linear straight-line program has input coordinates and gates ADD, SUB and +SCALE by prepared input-independent complex constants. Its scalar DAG depth +counts gates on the longest data path. Outputs are explicitly named registers. +Assume the following *checked properties of the actual supplied program*, not +an unproved stable-transform oracle: + +* the gate count and prepared/index work are at most N^u; +* every scale constant has modulus <=C0, with C0 fixed independently of N; +* the data DAG depth is at most a*log2(N+2), for fixed a; +* inputs have modulus <=N^b, and each prepared scalar has a certified error + <=epsilon supplied by an actual preparation algorithm. + +Then rounding all data gates with absolute local error <=epsilon produces +outputs with error <=epsilon*N^A for a fixed effectively computable A, and all +intermediate magnitudes are <=N^B for a fixed effectively computable B. +Thus p=O(log N) fractional and integer bits suffice for any prescribed inverse +polynomial output error. The rounded program has O(1) bounded-word operations +per original gate on words of a sufficiently large constant times log N bits. + +**Proof.** Let C=max(2,C0), and let h be a node's data depth. Exact gate values +are bounded by N^b*C^h by induction. With depth bound a log2(N+2), this is a +polynomial V=N^b*(N+2)^(a log2 C). For epsilon<=1, add/sub errors satisfy +E_new<=E_left+E_right+epsilon. A scale error satisfies + + E_new <= (C0+1)*E_input + epsilon*V + epsilon. + +The extra term comes from using an approximate prepared scalar on an +approximate data value. Induction gives + + E_at_depth_h <= epsilon*(V+1)*(C0+3)^(h+1). + +Its right side is polynomial in N times epsilon. The same estimate bounds +actual intermediate magnitudes. Fixing p at a sufficiently large constant +times log2(N+2), and allowing the proved polynomial integer headroom, keeps +all multiplications, sums, shifts and signed coordinates in a fixed number of +bounded words. ADD/SUB/SCALE then lower by explicit fixed-point arithmetic, +not by a simulation theorem. ∎ + +The theorem extends to a fixed number of such linear phases with a pointwise +bilinear product between them: if magnitudes are polynomially bounded, the +additional product error is bounded by |a|*E_b+|b|*E_a+E_a*E_b plus one local +rounding error. The recovery worker separately supplies the <1/2 rounding and +carry proof for integer convolution. + +### What this does and does not admit + +The native Gaussian quarter-Fourier C=((1+i)/2)I+((1-i)/2)X is unitary in +ell2; so are its Gaussian-unit monomial wrappers and tensor powers. A fixed +finite rational scalar motif therefore has a finite explicit norm factor K. +Composing it through L=O(log N) levels yields K^L=poly(N), giving precisely +the type of growth certificate above when its actual schedule/depth is bounded. +This is a sufficient component, not evidence that every exact-scalar provider +satisfies the certificate. In particular: + +* unrestricted scalar heights do not automatically have logarithmic-word + implementations; +* a division by a nonzero but extremely small prepared scalar can violate the + bound, and must be analyzed from its actual source; +* all-length sector/CRT/chirp preparation must be charged and bounded; +* a polynomially sized schedule may have linear depth, so a gate-count bound + alone is insufficient; +* the phase-table theorem supplies root powers, not arbitrary root-extraction + denominators or an algorithm for recovering huge arbitrary constants. + +PR #647 explicitly excludes numerical stability and permits nonzero +scalar-preparation divisions. PR #956 likewise declares unrestricted scalar +heights. Their native Gaussian providers motivate this lifting route, but their +full all-length algorithms are not admitted here by assumption. The root team +must bind the sufficient certificate to a particular actual source before +claiming a new native Word RAM multiplication bound. + +## 3. Evidence and reused sources + +`phase_tables.py` uses a direct signed-Gaussian-dyadic operational model and +explicit bounded-word overflow checks. Its 45 generated table configurations +and 774 queries are compared with exact rational unit-circle arithmetic; finite +checks include actual table corruption and insufficient-precision controls. +The rational exact targets are validation devices, not runtime root oracles. +All essential checks use explicit raises and survive optimized Python. + +Project precedents read in the shared canonical tree: +`Problem18RAMLax13Ram20261001`, `MpLimbs20261007`, `MpSparse20261007`, +`Problem18RAMRootLog20261007`, and the previous Word RAM packet. `MpLimbs.lv` +and `lv_dig` supply the established radix semantics; `add_run`, `map_run`, +`copyIn_run` and `copyOut_run` are the existing framed limb routines. Their +source/kernel status is revision-specific; this packet does not rerun Lean. + +Primary literature/model comparison: Fürer, *How Fast Can We Multiply Large +Integers on an Actual Computer?* (2014), Theorem 4.1, already gives O(n)-time +multiplication on its logarithmic-word RAM with word multiplication. Our +component is not that result's historical replacement or a new complexity +record. https://arxiv.org/html/1402.1811v1 + +The OpenAI exact-DFT companion explicitly permits unrestricted complex scalar +magnitudes and distinguishes its exact-register result from stable numerical +and bit-complexity algorithms. Its model contract must be honored when using +it as a precedent. https://github.com/openai/math/blob/main/preprints/Integer-multiplication-below-n-log-n-September-23-2026/build/sections/01-history.tex + +Prepared phases, guard precision, arithmetic, indexing, allocated buffers and IO +must all enter the complete caller invoice. No new numerical kappa is claimed +by this component packet. + +## 4. A quantitative bound for the actual six-positive-C shear denominators + +PR #956's all-length transfer names t1=1+mu*conjugate(mu) and t2=mu-t1. +For every complex mu, these satisfy + + |t1| >= 1, |t2| >= 3/4. + +In particular their inverses have modulus at most 1 and 4/3 respectively. +To prove this, write mu=a+ib. Then t1=1+a^2+b^2 is real and at least 1, while + + Re(t2)=a-1-a^2-b^2 + =-((a-1/2)^2+b^2+3/4) <= -3/4. + +Therefore |t2|>=|Re(t2)|>=3/4. These lower bounds are universal and effective; +there is no near-zero-denominator obstruction in these two named divisions. +If an implementation's denominator approximation has error at most 1/4, the +corresponding approximate denominators remain bounded away from zero, with +absolute values at least 3/4 and 1/2. Standard reciprocal perturbation then +bounds their inverse error by a constant times their denominator error. + +This does not bound the magnitudes of mu or of the other Newton/Toeplitz +coefficients. Those quantities require their own source-specific estimates. +It improves the sufficient finite-precision contract for this specific actual +shear decomposition without silently admitting the whole provider. diff --git a/research/pr4554-lowering-check/README.md b/research/pr4554-lowering-check/README.md new file mode 100644 index 0000000000..678cae6d8c --- /dev/null +++ b/research/pr4554-lowering-check/README.md @@ -0,0 +1,38 @@ +# Independent check of #4554: finite-precision lowering of linear programs (`PHASE-TABLES.md` §2, §4) + +**Task:** independent-moment (independent source-bound coverage). **Base:** canonical main `f626db5c`. **No new κ.** + +This is the last component of Dugongue's #4554 (head `4f46c6a9`) that had no independent check; #4556–#4562 cover the rest. `PHASE-TABLES.md` §2 proves a generic sufficient lowering. A linear straight-line program with scale constants `|c| ≤ C0`, data depth `≤ a·log₂(N+2)`, inputs `≤ N^b` and prepared-scalar errors `≤ ε`, with every data gate rounded within `ε`, has node errors at most `ε(V+1)(C0+3)^(h+1)` and magnitudes at most `N^b C^h`. §4 bounds the perturbed shear denominators. The text is vendored as `PR4554-PHASE-TABLES.md` (blob `2dab8f3e`, SHA256-pinned). + +## Results + +All checked claims are confirmed. + +1. **Induction steps, exact rationals.** For every `h ≤ 200` and `C0 ∈ {1/2, 1, 3/2, 2, 7/3, 10}` (1,200 cases): + - `2C^(h−1)` and `C0·C^(h−1)` are both at most `C^h`, with `C = max(2, C0)`; + - the closed form `B_h = ε(V+1)(C0+3)^(h+1)` dominates both the add step `2B_(h−1) + ε` and the scale step `(C0+1)B_(h−1) + εV + ε`. +2. **Random programs.** 60 random layered ADD/SUB/SCALE programs (`N ≤ 32`, `b ≤ 1`, `C0 ≤ 2`, `a ≤ 2`, `P ∈ {16, 24, 32}`) ran in exact rationals against a copy with perturbed prepared scalars and gates rounded toward zero. At every node the error is at most `ε(V+1)(C0+3)^(h+1)` (worst ratio 0.0094), and the exact magnitude is at most `N^b C^h`. +3. **FFT instantiation.** On radix-2 FFT DAGs (depth `2 log₂ N`, `C0 = 1`, `N` up to 512), the observed rounded-twiddle error lies below the sharper `FFT-ERROR.md` chain (checked in #4557), which in turn lies below the §2 bound. **Observation:** the §2 bound is valid but very loose here, 2·10^9 against an observed 3·10^−6 at `N = 512`. That fits its role as a generic sufficient certificate: it shows `O(log N)` bits suffice, not how many. +4. **Bilinear extension and §4.** + - `|âb̂ − ab| ≤ |a|E_b + |b|E_a + E_aE_b` holds in 2,000 random cases. + - Shear denominators perturbed by up to 1/4 keep `|t1| ≥ 3/4` and `|t2| ≥ 1/2`. The exact identities behind this are checked in #4562. +5. **Binding.** 5 statements are bound verbatim to the vendored text. + +## Verify + +```text +python3 -B research/pr4554-lowering-check/verify.py +``` + +This runs in about 1 s with the standard library only. + +## Limitations + +- **Finite checks.** The random programs and FFT DAGs are finite. The universal statement rests on the written induction, whose steps are checked exactly. +- **Floating point.** FFT comparisons use double-precision floats, whose rounding is far below the compared bounds. + +## Open next questions + +- None new for #4554's constructive components; with #4556–#4562, each now has an independent check. The PR's own stated gap remains: a faster native cyclic transform satisfying `NATIVE-OUTER-ASSEMBLY.md`'s contract. + +Credit: Dugongue (#4554). Prepared by Ryanson Jonathan with substantial Claude Code (Claude Opus 5.5) assistance. diff --git a/research/pr4554-lowering-check/SOURCES.json b/research/pr4554-lowering-check/SOURCES.json new file mode 100644 index 0000000000..a2a5ea8bbd --- /dev/null +++ b/research/pr4554-lowering-check/SOURCES.json @@ -0,0 +1,18 @@ +{ + "schema_version": 1, + "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "ed33caef5c5827693e9572d252a66f99cbb73d40", + "sha256": "5732707cc1e460b090b53f2463162749ad860de878a18f7408938671d295a9d9" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "64bfc5af6709ed9811171e0a2108f05e1219cb65", + "sha256": "352a6184713b8843d96fd2c2071eb7aeaa655906b7c51d45ad444912666c8a96" + } + ] +} diff --git a/research/pr4554-lowering-check/TASK.json b/research/pr4554-lowering-check/TASK.json new file mode 100644 index 0000000000..146c582cc1 --- /dev/null +++ b/research/pr4554-lowering-check/TASK.json @@ -0,0 +1,100 @@ +{ + "schema_version": 1, + "task_id": "independent-moment", + "priority_revision": 7, + "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1", + "original_card": { + "id": "independent-moment", + "question": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?", + "starting_points": [ + "docs/research/research-family-status.md", + "certificates/selected-result.json", + "research/staged544-mixed-fourth-entry/README.md", + "research/staged-fourth-mixture-ceiling/README.md", + "research/r437-two-window-carrier/README.md" + ], + "methods": [ + "Independently reconstruct one changed paid ledger or pooled dual with exact source and model binding.", + "Check a genuinely uncovered moment implication with a different enclosure method and semantic controls." + ], + "acceptance_gate": "Demonstrate independence, source binding, positive cases and semantic failure controls. Name the existing checks and the missing coverage. More digits or another equivalent evaluation of an unchanged root is not by itself new verification coverage.", + "first_step": { + "command": [ + "python3", + "-B", + "scripts/research_baseline.py" + ], + "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY", + "inspect": "Choose one uncovered implication in #4088 literal mixed fill, #4091 pooled dual/model binding or #4067 changed frame/host bill. #4091 pins an earlier #4088 head: validate word, capacity and stock equivalence before applying that ceiling to the later fill.", + "scope": "Source-binding orientation only; no scientific replay, theorem acceptance or supplier admission. Budget any heavy replay separately." + }, + "first_experiment": "Independently reconstruct one changed ledger or pooled certificate with a positive case and a semantic failure control. Compare exact sources and name what existing checks miss. Keep this separate from supplier-boundary-audit (ordinary/outer dispatch) and uniform-transfer-audit (all-size proof). Stop after the chosen implication; more digits of an unchanged root are not additional mechanism coverage.", + "bottleneck_link": "Independent pricing prevents unsupported frontier claims and can reveal genuine accounting slack.", + "allowed_variations": [ + "Reconstruct either a source ledger or a numeric enclosure with a different method." + ], + "scope_change_triggers": [ + "Importing the primary function and calling its output independent verification.", + "Calling an aggregate histogram or source archive an independently accepted physical supplier." + ], + "overlap_terms": [ + "independent moment", + "rank histogram", + "exact enclosure", + "2979", + "rank7 joint endframes", + "paired extraction", + "joint lowering", + "4088", + "4091", + "4067", + "staged544-mixed-fourth-entry", + "staged-fourth-mixture-ceiling", + "r437-two-window-carrier" + ], + "existing_checks": [ + { + "scope": "PR #4088: mixed local/owner/phase ledger and moment reported.", + "evidence": "research/staged544-mixed-fourth-entry/README.md", + "remaining": "Independent binding of changed construction to paid histogram, distinct from full admission.", + "status": "author-checked", + "source_commit": "72d7efc181edb21806a878f174847c96dd410ab0", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "PR #4091: sharp pooled ceiling proposed for fixed THREE/FOUR words and inventory.", + "evidence": "research/staged-fourth-mixture-ceiling/README.md", + "remaining": "Validate dual, exhaustive pattern scope and exact source equivalence before applying to newer fills.", + "status": "author-checked", + "source_commit": "a29bdd8f56cd18b805917775cdec09ff935eddc4", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "PR #4067: paid (2,2) to (1,3) exchange and fixed OQ4 host reported.", + "evidence": "research/r437-two-window-carrier/README.md", + "remaining": "Independent full cost delta; remaining caller/precision/stock obligations are separate.", + "status": "author-checked", + "source_commit": "a95f11c52e3ca1e4c5171b3c20d8f4b0815aee2e", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + } + ] + }, + "scope_status": "unchanged", + "actual_scope": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?", + "scope_change_reason": null, + "bottleneck_link": "Completes independent coverage of every constructive component of #4554 (with #4556-#4562).", + "budget_minutes": 30, + "stop_rule": "Stop at PHASE-TABLES.md sections 2 and 4.", + "outcome": "completed", + "evidence": "Exact induction steps for h<=200 (1,200 cases); 60 random layered programs within eps(V+1)(C0+3)^(h+1) (worst 0.0094) and N^b C^h; FFT DAG instantiation (observed < FFT-ERROR chain < section-2 bound); bilinear extension and perturbed shear denominators (2,000 cases); 5 statements bound.", + "author": "Ryanson Jonathan, with Claude Code (Claude Opus 5.5) assistance", + "ai_assistance": "Substantial Claude Code (Claude Opus 5.5) assistance; no independent human review.", + "open_questions": [ + { + "question": "A faster native cyclic transform satisfying NATIVE-OUTER-ASSEMBLY.md's contract.", + "why_it_matters": "#4554's stated missing hypothesis for any Word RAM kappa.", + "first_step": "Bind the #956 rational kernel lowering (#4562) to a cyclic embedding with paid preparation.", + "blocked_by": "Newton/Toeplitz preparation bounds." + } + ] +} diff --git a/research/pr4554-lowering-check/VERIFY.json b/research/pr4554-lowering-check/VERIFY.json new file mode 100644 index 0000000000..1b454b9501 --- /dev/null +++ b/research/pr4554-lowering-check/VERIFY.json @@ -0,0 +1,18 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [], + "commands": [ + { + "argv": [ + "{python}", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 120, + "checks": "Exact induction steps; random rounded linear programs vs per-node error and magnitude bounds; FFT DAG instantiation vs FFT-ERROR chain; bilinear extension and perturbed shear denominators; binding to vendored spec (sha256)." + } + ] +} diff --git a/research/pr4554-lowering-check/verify.py b/research/pr4554-lowering-check/verify.py new file mode 100644 index 0000000000..357c937a89 --- /dev/null +++ b/research/pr4554-lowering-check/verify.py @@ -0,0 +1,212 @@ +#!/usr/bin/env python3 +"""Independent check of PR #4554 (Dugongue) PHASE-TABLES.md section 2 (finite-precision lowering of linear +straight-line programs) and section 4 (perturbed shear denominators), from the vendored text only. + + 1. Induction steps, exactly in rationals for every depth h <= 200 and C0 in a grid: with C = max(2, C0), + 2 N^b C^(h-1) <= N^b C^h and C0 N^b C^(h-1) <= N^b C^h; with X = C0 + 3 and the closed form + B_h = eps (V+1) X^(h+1): 2 B_(h-1) + eps <= B_h (add) and (C0+1) B_(h-1) + eps V + eps <= B_h (scale). + 2. Random layered programs (ADD, SUB, SCALE by prepared complex scalars |c| <= C0 whose prepared values carry an + error <= eps; inputs |x| <= N^b; depth <= a log2(N+2)) run in exact rationals and in a rounded copy (prepared + scalars perturbed, every data gate rounded toward zero to 2^-P): at every node the error is <= eps(V+1)(C0+3)^(h+1) + and the exact value is <= N^b C^h, where h is the node's data depth. + 3. Instantiation on radix-2 FFT DAGs (depth 2 log2 N, C0 = 1, inputs <= beta): the theorem's bound and the sharper + FFT-ERROR.md chain (checked in #4557) both dominate the observed rounded-twiddle FFT error. + 4. The bilinear product extension |a_hat b_hat - a b| <= |a| E_b + |b| E_a + E_a E_b (+ one rounding), and section 4: + approximate shear denominators within 1/4 keep |t1| >= 3/4 and |t2| >= 1/2. +Ryanson Jonathan with Claude Code (Claude Opus 5.5) assistance. Apache-2.0. +""" +from fractions import Fraction as Fr +import cmath +import hashlib +import json +import math +from pathlib import Path +import random +import sys + +HERE = Path(__file__).resolve().parent +SPEC = HERE / "PR4554-PHASE-TABLES.md" # #4554 head 4f46c6a9, research/native-word-fft/PHASE-TABLES.md, blob 2dab8f3e +SPEC_SHA256 = "9a5844176b67f81ad8468b77d0049935d05517a80e80586edfa8401f400917f8" + + +def need(condition, message): + if not condition: + raise AssertionError(message) + + +def check_induction(): + cases = 0 + for C0 in (Fr(1, 2), Fr(1), Fr(3, 2), Fr(2), Fr(7, 3), Fr(10)): + C = max(Fr(2), C0) + X = C0 + 3 + for h in range(1, 201): + need(2 * C ** (h - 1) <= C ** h and C0 * C ** (h - 1) <= C ** h, "magnitude induction") + # With eps and V scaled out (eps = 1, V >= 0 arbitrary): check for V = 0 and a large V. + for V in (Fr(0), Fr(10 ** 6)): + Bh, Bp = (V + 1) * X ** (h + 1), (V + 1) * X ** h + need(2 * Bp + 1 <= Bh, "add-error induction") + need((C0 + 1) * Bp + V + 1 <= Bh, "scale-error induction") + cases += 1 + return cases + + +def tz(v, P): + S = 1 << P + n = v * S + q = n.numerator // n.denominator if n >= 0 else -((-n.numerator) // n.denominator) + return Fr(q, S) + + +def cabs(z): + return math.hypot(float(z[0]), float(z[1])) + + +def cadd(x, y, sign=1): + return (x[0] + sign * y[0], x[1] + sign * y[1]) + + +def cmul(x, y): + return (x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]) + + +def check_random_programs(rng): + runs, worst = 0, 0.0 + for _ in range(60): + N = rng.choice((4, 8, 16, 32)) + b = rng.choice((0, 1)) + C0 = rng.choice((Fr(1, 2), Fr(1), Fr(2))) + a = rng.choice((1, 2)) + depth = int(a * math.log2(N + 2)) + P = rng.choice((16, 24, 32)) + eps = Fr(1, 1 << (P - 1)) # covers the complex rounding of both parts (sqrt(2)/2^P < 2^(1-P)) + width = rng.randint(2, 5) + nb = Fr(N) ** b + layer = [(Fr(rng.randint(-int(nb) * 8, int(nb) * 8), 8), Fr(rng.randint(-int(nb) * 8, int(nb) * 8), 8)) + for _ in range(width)] + layer = [(max(min(x, nb / 2), -nb / 2), max(min(y, nb / 2), -nb / 2)) for x, y in layer] + exact = list(layer) + approx = list(layer) + h_of = [0] * width + C = max(Fr(2), C0) + V = nb * C ** depth + for d in range(1, depth + 1): + new_e, new_a, new_h = [], [], [] + for _ in range(width): + op = rng.choice(("ADD", "SUB", "SCALE")) + i, j = rng.randrange(width), rng.randrange(width) + if op == "SCALE": + ang = rng.uniform(0, 2 * math.pi) + r = float(C0) * rng.random() + c = (Fr(r * math.cos(ang)).limit_denominator(1 << 20), Fr(r * math.sin(ang)).limit_denominator(1 << 20)) + if cabs(c) > float(C0): + c = (c[0] / 2, c[1] / 2) + ch = (c[0] + Fr(rng.randint(-1000, 1000), 1000) * eps / 2, c[1] + Fr(rng.randint(-1000, 1000), 1000) * eps / 2) + ve, va = cmul(c, exact[i]), cmul(ch, approx[i]) + else: + sign = 1 if op == "ADD" else -1 + ve, va = cadd(exact[i], exact[j], sign), cadd(approx[i], approx[j], sign) + va = (tz(va[0], P), tz(va[1], P)) + new_e.append(ve) + new_a.append(va) + new_h.append(d) + exact, approx, h_of = new_e, new_a, new_h + for k in range(width): + bound = float(eps * (V + 1)) * float(C0 + 3) ** (h_of[k] + 1) + err = math.hypot(float(approx[k][0] - exact[k][0]), float(approx[k][1] - exact[k][1])) + need(err <= bound, "error <= eps (V+1)(C0+3)^(h+1)") + need(cabs(exact[k]) <= float(nb * C ** h_of[k]) * (1 + 1e-12), "magnitude <= N^b C^h") + worst = max(worst, err / bound) + runs += 1 + return runs, worst + + +def check_fft(rng): + rows = [] + for logN in (3, 5, 7, 9): + N = 1 << logN + beta = 15 + P = 30 + eps = 2.0 ** (1 - P) + x = [complex(rng.randint(0, beta), 0) for _ in range(N)] + + def fft(a, perturb): + a = list(a) + j = 0 + for i in range(1, N): + bit = N >> 1 + while j & bit: + j ^= bit + bit >>= 1 + j |= bit + if i < j: + a[i], a[j] = a[j], a[i] + m = 2 + while m <= N: + for st in range(0, N, m): + for k in range(m // 2): + w = cmath.exp(-2j * math.pi * k / m) + if perturb: + w += complex(rng.uniform(-1, 1), rng.uniform(-1, 1)) * eps / 2 + q = 2.0 ** P + t = w * a[st + k + m // 2] + t = complex(math.trunc(t.real * q) / q, math.trunc(t.imag * q) / q) + else: + t = w * a[st + k + m // 2] + u = a[st + k] + a[st + k], a[st + k + m // 2] = u + t, u - t + m *= 2 + return a + + err = max(abs(p - q) for p, q in zip(fft(x, True), fft(x, False))) + depth = 2 * logN + V = beta * 2 ** depth # N^b C^h with C = 2 for unit twiddles, inputs <= beta + theorem = eps * (V + 1) * 4 ** (depth + 1) + chain = 36 * N ** 3 * beta * 2.0 ** (-P) * 4 # FFT-ERROR (1) at S = 2^(P-2) (coordinate error eps/2) + need(err <= chain <= theorem, "observed <= FFT-ERROR chain <= section-2 theorem") + rows.append({"N": N, "observed": err, "fft_error_chain": chain, "section2_theorem": theorem}) + return rows + + +def check_bilinear_and_shear(rng): + cases = 0 + for _ in range(2000): + a = complex(rng.uniform(-50, 50), rng.uniform(-50, 50)) + b = complex(rng.uniform(-50, 50), rng.uniform(-50, 50)) + Ea, Eb = rng.uniform(0, 1), rng.uniform(0, 1) + ah = a + Ea * cmath.exp(2j * math.pi * rng.random()) + bh = b + Eb * cmath.exp(2j * math.pi * rng.random()) + need(abs(ah * bh - a * b) <= abs(a) * Eb + abs(b) * Ea + Ea * Eb + 1e-9, "bilinear extension") + mu = complex(rng.uniform(-20, 20), rng.uniform(-20, 20)) + t1 = 1 + abs(mu) ** 2 + t2 = mu - t1 + for _ in range(3): + e1 = 0.25 * rng.random() * cmath.exp(2j * math.pi * rng.random()) + e2 = 0.25 * rng.random() * cmath.exp(2j * math.pi * rng.random()) + need(abs(t1 + e1) >= 0.75 - 1e-12 and abs(t2 + e2) >= 0.5 - 1e-12, "perturbed shear denominators") + cases += 1 + return cases + + +def check_binding(): + raw = SPEC.read_bytes() + need(hashlib.sha256(raw).hexdigest() == SPEC_SHA256, "vendored PHASE-TABLES.md") + text = raw.decode("utf-8") + for s in ("E_new<=E_left+E_right+epsilon.", "E_new <= (C0+1)*E_input + epsilon*V + epsilon.", + "E_at_depth_h <= epsilon*(V+1)*(C0+3)^(h+1).", "|a|*E_b+|b|*E_a+E_a*E_b", + "absolute values at least 3/4 and 1/2."): + need(s in text, f"spec states: {s}") + return 5 + + +def main(): + rng = random.Random(45546) + runs, worst = check_random_programs(rng) + report = {"status": "PASS_PR4554_LOWERING_CHECK", "induction_cases": check_induction(), + "random_programs": runs, "worst_error_over_bound": worst, "fft_instances": check_fft(rng), + "bilinear_and_shear_cases": check_bilinear_and_shear(rng), "spec_statements_bound": check_binding(), + "new_kappa_claimed": False} + print(json.dumps(report, indent=1)) + + +if __name__ == "__main__": + sys.exit(main())