From 4ed34481357b2331be54da8b2661e38638658438 Mon Sep 17 00:00:00 2001 From: pinterbanget Date: Thu, 22 Oct 2026 13:38:31 +0700 Subject: [PATCH] =?UTF-8?q?Independently=20check=20#4554's=20paid=20phase?= =?UTF-8?q?=20and=20exp=20tables:=20own=20tables=20vs=2016(E+d)/Q=20and=20?= =?UTF-8?q?2^-p,=20worst-direction=20seeds=20(no=20new=20=CE=BA)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- .../pr4554-native-tables-check/.gitattributes | 2 + .../PR4554-NATIVE-EXP-TABLE.md | 31 ++ .../PR4554-PHASE-TABLES.md | 221 ++++++++++++++ research/pr4554-native-tables-check/README.md | 48 +++ .../pr4554-native-tables-check/SOURCES.json | 18 ++ research/pr4554-native-tables-check/TASK.json | 100 ++++++ .../pr4554-native-tables-check/VERIFY.json | 18 ++ research/pr4554-native-tables-check/verify.py | 286 ++++++++++++++++++ 8 files changed, 724 insertions(+) create mode 100644 research/pr4554-native-tables-check/.gitattributes create mode 100644 research/pr4554-native-tables-check/PR4554-NATIVE-EXP-TABLE.md create mode 100644 research/pr4554-native-tables-check/PR4554-PHASE-TABLES.md create mode 100644 research/pr4554-native-tables-check/README.md create mode 100644 research/pr4554-native-tables-check/SOURCES.json create mode 100644 research/pr4554-native-tables-check/TASK.json create mode 100644 research/pr4554-native-tables-check/VERIFY.json create mode 100644 research/pr4554-native-tables-check/verify.py diff --git a/research/pr4554-native-tables-check/.gitattributes b/research/pr4554-native-tables-check/.gitattributes new file mode 100644 index 0000000000..46a4f09e04 --- /dev/null +++ b/research/pr4554-native-tables-check/.gitattributes @@ -0,0 +1,2 @@ +PR4554-PHASE-TABLES.md -text +PR4554-NATIVE-EXP-TABLE.md -text diff --git a/research/pr4554-native-tables-check/PR4554-NATIVE-EXP-TABLE.md b/research/pr4554-native-tables-check/PR4554-NATIVE-EXP-TABLE.md new file mode 100644 index 0000000000..ee6e734c82 --- /dev/null +++ b/research/pr4554-native-tables-check/PR4554-NATIVE-EXP-TABLE.md @@ -0,0 +1,31 @@ +# Paid constant-time native Gaussian weights + +Fix positive constants c, C and 0c+eta+2. Above a fixed finite cutoff, the alternating Taylor polynomial + + q0 = sum_(k=0)^d (-h)^k/k! + +satisfies 0<=q0<=exp(-h) and exp(-h)-q0<=2^-S. The degree is fixed independently of n. Fixed rational polynomial evaluation uses a fixed number of bounded-word operations; compute it with enough guard bits that its total error is below 2^-S, subtract that certified error margin, then round downward to obtain 0<=qhat<=exp(-h) with exp(-h)-qhat<=C_d 2^-S for an explicit template constant C_d. Increase c0 to absorb C_d and every fixed polynomial-evaluation error constant. The initial finite cutoff can use a direct table constructed with higher precision; it depends only on fixed program constants. + +Set Ehat_0=1 and Ehat_(j+1)=floor_S(Ehat_j qhat), where floor_S is rounding downward to the nonnegative binary S-bit grid. Since 0<=qhat<=exp(-h)<=1, every Ehat_j lies in [0,exp(-jh)]. Induction gives + + 0<=exp(-jh)-Ehat_j<=j(C_d+1)2^-S. + +There are J+1=O(n^eta log n) entries, and each recurrence update is bounded-word multiplication, shift and store. Choosing c0 so that (C_d+1)J2^-S<=2^(-p-4) controls the whole prepared table. + +For a query, use binary shifts to compute j=floor(x/h) and u=x-jh in [0,h). Evaluate the same fixed Taylor polynomial at u to error at most 2^(-p-4), multiply by Ehat_j, and round the result with the same guard precision. Then exp(-x)=exp(-jh)exp(-u), all exact factors have magnitude at most 1, and the sum of table, residual and product-rounding errors is less than 2^-p. If x itself approximates a nonnegative intended argument within epsilon, exp(-x) changes by at most epsilon (or exp(epsilon) epsilon if the approximation crosses zero); budget this explicitly when constructing the query. + +All rational coefficients 1/k! have fixed-size denominators. Every stored scalar and exact intermediate has O(log n) bits. J and every address also have O(log n) bits. Setup includes all table entries, and query evaluation includes the residual calculation. There is no multiplication or division on a full n-bit operand, no Turing machine and no root/exponential oracle. + +## Actual Gaussian caller + +The finite-precision resampling formulas in upstream build/sections/07-resampling.tex use coefficients exp(-a z^2) on windows for which a z^2=O(p), after factoring their explicitly stated normalizations. Compute a z^2 by bounded-word arithmetic and invoke the prepared table. Arguments beyond a chosen Cp cutoff can return zero, with exp(-Cp) below the assigned tail budget; the cutoff C must be chosen from that caller's requested error. A fixed approximation to pi at S+O(log p) bits, obtainable by a separately charged ordinary series, suffices for a=pi/alpha^2 or pi alpha^2 after its actual polynomial argument-size bound is included. Polynomially large positive factors can be treated by reversal against a prepared exp(-x) scale, with an additional O(log n) magnitude and precision allowance; this sentence does not replace the actual caller's bound. + +Thus the explicit window-product loops of the S map and E=N-I map can prepare each negative Gaussian weight in O(1) native instructions after the paid O(n^eta log n) table. This justifies a direct arithmetic-loop recount rather than dividing a tape-time estimate by p. It does not itself accelerate the number of window visits, prove a blocked chirp convolution, or establish a new complete multiplication exponent. + +The exponential table construction is a standard range-reduction/table idea specialized here to the actual O(log n)-argument and O(log n)-precision regime. No worldwide novelty claim is made. The Gaussian resampling identities and stability estimates remain credited to the upstream source and its Harvey-Hoeven references. + +## Signed-argument corollary for block chirps + +For 0<=x<=C log2(n+2), exp(x)=1/exp(-x). Request the negative exponential with absolute error epsilon<=2^(-p-3) exp(-2C log2(n+2)); use an explicit power-of-two upper bound for this budget. Since exp(-x)>=exp(-C log2(n+2)), its approximation remains at least half that lower bound. The reciprocal identity bounds the resulting reciprocal error by 2 exp(2C log2(n+2)) epsilon, before the final bounded-word division rounding. Reserve a further 2^(-p-2) budget for that rounding. The added fractional bits and positive-result integral bits are both O(log n). Thus the same paid construction supports exp(x) for the full signed interval |x|<=C log2(n+2), with O(1) bounded-word instructions per query. This supplies the possibly positive quadratic weights in a centered Gaussian block/Hankel factorization. The constant-factor word width is fixed before n. diff --git a/research/pr4554-native-tables-check/PR4554-PHASE-TABLES.md b/research/pr4554-native-tables-check/PR4554-PHASE-TABLES.md new file mode 100644 index 0000000000..2dab8f3e6d --- /dev/null +++ b/research/pr4554-native-tables-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-native-tables-check/README.md b/research/pr4554-native-tables-check/README.md new file mode 100644 index 0000000000..94de43d1de --- /dev/null +++ b/research/pr4554-native-tables-check/README.md @@ -0,0 +1,48 @@ +# Independent check of #4554: paid native phase and exp tables + +**Task:** independent-moment (independent source-bound coverage). **Base:** canonical main `f626db5c`. **No new κ.** + +Dugongue's #4554 (head `4f46c6a9`) prepares two kinds of table: +- root powers `z^e`, with query error `16(E+d)/Q` and O(d) cost per query, in `PHASE-TABLES.md` §1; +- negative Gaussian weights `exp(−x)`, with error `< 2^−p`, O(1) cost per query and `O(n^η log n)` setup, in `NATIVE-EXP-TABLE.md`. + +This package checks both, after my #4556, #4557 and #4558 on the same PR. The texts are vendored as `PR4554-PHASE-TABLES.md` (blob `2dab8f3e`) and `PR4554-NATIVE-EXP-TABLE.md` (blob `ee6e734c`), with their SHA256 pinned. Nothing is imported from #4554's `phase_tables.py` or `exp_table.py`. + +## Results + +All checked claims are confirmed. + +**1. Root-power tables (§1).** The construction follows the spec: `t` rounded squarings per level, `T_j[a] = round(T_j[a−1]·z_j)`, and queries that multiply the d digit-selected entries starting from exact 1, with nearest rounding to `1/Q`. +- **Coverage.** Six `(D, t)` pairs from `(12, 1)` to `(2^20, 4)`. Each uses all four floor/ceil grid seeds around `exp(2πi/D)`, so every coordinate error is `< 1/Q`, and `Q` is the smallest power of two with `Q ≥ 64(E+d)`. +- **Query error.** Over 20,852 queries, every error is at most `16(E+d)/Q` against a decimal reference with AGM π. The worst is 0.076 of the bound. +- **Magnitudes.** Coordinates stay `≤ 2Q`, products `≤ 4Q²` and sums `≤ 8Q²`, all within `2p+8` signed bits. +- **Error-factor count.** A literal unfolding of the computation tree confirms at most `2E+d` error factors. `(2E+d)η ≤ 1/16`, so `(1+η)^(2E+d) − 1 ≤ 16(E+d)/Q`, checked exactly. + +**2. Exp table.** Concrete constants are `(c, C, η) = (1, 2, 1/4)` and odd `d = 13`, which gives `η(d+1) = 3.5 > c+η+2 = 3.25`, with `c0 = 6`. The test sizes are `n = 2^16`, `2^20` and `2^24`, with up to 3,073 table entries. +- **Taylor base.** The exact rational Taylor polynomial, rounded down to S bits, satisfies `0 ≤ q̂ ≤ e^(−h)`, with `C_d ≤ 0.96`. Then `(C_d+1)J·2^(−S) ≤ 2^(−p−4)`. +- **Table entries.** Every entry satisfies `0 ≤ e^(−jh) − Ê_j ≤ j(C_d+1)2^(−S)`; the worst is 0.70 of the bound. +- **Queries.** 204 queries per size, at the edges, at multiples of `h` and at random points, all have error `< 2^−p`. The worst is 1.8·10⁻⁴ of `2^−p`. +- **Cutoff.** For these constants, the tail condition `h^(d+1)/(d+1)! ≤ 2^(−S)` holds for every `n = 2^m − 2` with `2 ≤ m ≤ 4096`. The "fixed finite cutoff" is empty here. +- **Signed corollary.** `exp(+x) = 1/exp(−x)` with the stated budget `ε ≤ 2^(−p−3)e^(−2C log₂(n+2))` and a `2^(−p−2)` reciprocal rounding gives error `< 2^−p` (worst 0.12). + +**3. Binding.** 10 statements are bound verbatim to the vendored texts. + +## Verify + +```text +python3 -B research/pr4554-native-tables-check/verify.py +``` + +This runs in a few seconds with the standard library only. + +## Limitations + +- **Out of scope:** §2 of `PHASE-TABLES.md` (finite-precision lowering of linear programs) and the Gaussian caller's argument bounds. +- **Exp table constants.** The exp table is checked for one concrete `(c, C, η, d, c0)` and three sizes. The universal statement rests on the written recurrence, which these checks corroborate. +- **Reference.** References are decimal, at well above the needed precision; they are not interval-certified. + +## Open next questions + +- An independent check of `PHASE-TABLES.md` §2, the finite-precision lowering of linear straight-line programs, on an actual transform DAG. First step: instantiate it on a radix-2 FFT DAG and compare the `ε·N^A` bound with the exact error chain in my #4557. + +Credit: Dugongue (#4554). Prepared by Ryanson Jonathan with substantial Claude Code (Claude Opus 5.5) assistance. diff --git a/research/pr4554-native-tables-check/SOURCES.json b/research/pr4554-native-tables-check/SOURCES.json new file mode 100644 index 0000000000..a2a5ea8bbd --- /dev/null +++ b/research/pr4554-native-tables-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-native-tables-check/TASK.json b/research/pr4554-native-tables-check/TASK.json new file mode 100644 index 0000000000..9f017e44ef --- /dev/null +++ b/research/pr4554-native-tables-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": "Independent coverage of the paid phase/exp tables that the native Gaussian adapters query in O(1).", + "budget_minutes": 45, + "stop_rule": "Stop at PHASE-TABLES.md section 1 and NATIVE-EXP-TABLE.md.", + "outcome": "completed", + "evidence": "Phase tables: 20,852 queries within 16(E+d)/Q (worst 0.076) over 6 (D,t) and 4 worst-direction seeds, magnitude envelope, exact 2E+d factor count; exp table at (c,C,eta,d,c0)=(1,2,1/4,13,6) for n=2^16,2^20,2^24: qhat bracket, every table entry within j(C_d+1)2^-S, 612 queries < 2^-p, empty cutoff, reciprocal corollary; 10 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": "Independent check of PHASE-TABLES.md section 2 (linear-program lowering) on an actual transform DAG.", + "why_it_matters": "It is the generic precision theorem for native transforms.", + "first_step": "Instantiate on a radix-2 FFT DAG and compare with the exact error chain of #4557.", + "blocked_by": "None." + } + ] +} diff --git a/research/pr4554-native-tables-check/VERIFY.json b/research/pr4554-native-tables-check/VERIFY.json new file mode 100644 index 0000000000..b3651f5a2c --- /dev/null +++ b/research/pr4554-native-tables-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": 300, + "checks": "Own root-power tables vs 16(E+d)/Q with worst-direction seeds and magnitude envelope; own exp(-x) table recurrence, query error < 2^-p, cutoff and reciprocal corollary; binding to vendored texts (sha256)." + } + ] +} diff --git a/research/pr4554-native-tables-check/verify.py b/research/pr4554-native-tables-check/verify.py new file mode 100644 index 0000000000..db1e952f85 --- /dev/null +++ b/research/pr4554-native-tables-check/verify.py @@ -0,0 +1,286 @@ +#!/usr/bin/env python3 +"""Independent check of PR #4554 (Dugongue) paid native tables: PHASE-TABLES.md section 1 and NATIVE-EXP-TABLE.md. + +Written from the vendored texts only; nothing is imported from #4554's phase_tables.py or exp_table.py. + 1. Root-power tables. z_j by t rounded squarings per level, T_j[a] = round(T_j[a-1] z_j), queries multiply the d + digit-selected entries starting from exact 1; nearest rounding to the 1/Q grid (half away from zero). Seeds are + all four floor/ceil grid neighbours of exp(2 pi i/D), so each coordinate error is < 1/Q. For Q >= 64(E+d), every + query error is compared with 16(E+d)/Q against a high-precision decimal reference (AGM pi), and coordinates, + products and sums are checked against 2Q, 4Q^2, 8Q^2 (2p+8 signed bits). The leaf/node count 2E+d and the + inequality (1+eta)^(2E+d) - 1 <= 16(E+d)/Q are checked exactly. + 2. Exp table, with concrete (c, C, eta) = (1, 2, 1/4), odd d = 13 (eta(d+1) = 3.5 > c+eta+2 = 3.25) and c0 = 6: + h = 2^-floor(eta log2(n+2)), J = ceil(C log2(n+2)/h), S = p + ceil(log2(J+1)) + c0. The base qhat is the exact + rational Taylor polynomial rounded down to S bits; 0 <= qhat <= e^-h and e^-h - qhat <= C_d 2^-S are checked, and + every table entry satisfies 0 <= e^(-jh) - Ehat_j <= j(C_d+1)2^-S. Queries x in [0, C log2(n+2)] (edges, multiples + of h, random) have error < 2^-p. The cutoff above which h^(d+1)/(d+1)! <= 2^-S holds is located exactly. + The reciprocal corollary for exp(+x) is checked with the stated budgets. +Ryanson Jonathan with Claude Code (Claude Opus 5.5) assistance. Apache-2.0. +""" +import decimal +from decimal import Decimal as Dec +from fractions import Fraction as Fr +import hashlib +import json +import math +from pathlib import Path +import random +import sys + +HERE = Path(__file__).resolve().parent +SPECS = {"PR4554-PHASE-TABLES.md": "9a5844176b67f81ad8468b77d0049935d05517a80e80586edfa8401f400917f8", "PR4554-NATIVE-EXP-TABLE.md": "c9be6814a81e0bfe6b041f1a259a32369f7c030ca22eaa5c71c15911c888c6ff"} # #4554 head 4f46c6a9; blobs 2dab8f3e (phase), ee6e734c (exp) + + +def need(condition, message): + if not condition: + raise AssertionError(message) + + +def agm_pi(digits): + decimal.getcontext().prec = digits + 10 + a, b, t, pw = Dec(1), Dec(1) / Dec(2).sqrt(), Dec(1) / Dec(4), Dec(1) + for _ in range(int(math.log2(digits)) + 4): + an = (a + b) / 2 + b = (a * b).sqrt() + t -= pw * (a - an) ** 2 + a = an + pw *= 2 + return (a + b) ** 2 / (4 * t) + + +def dec_cos_sin(x, digits): + decimal.getcontext().prec = digits + 10 + c, s, term, k = Dec(0), Dec(0), Dec(1), 0 + eps = Dec(10) ** (-(digits + 8)) + while abs(term) > eps or k < 4: + if k % 4 == 0: + c += term + elif k % 4 == 1: + s += term + elif k % 4 == 2: + c -= term + else: + s -= term + k += 1 + term = term * x / k + return c, s + + +def dec_exp(x, digits): + decimal.getcontext().prec = digits + 10 + return Dec(x).exp() if not isinstance(x, Fr) else (Dec(x.numerator) / Dec(x.denominator)).exp() + + +# ---------------------------------------------------------------- 1. root-power tables +class Mag: + def __init__(self, Q, p): + self.Q, self.p = Q, p + self.coord = self.prod = self.sum = 0 + + def mul(self, a, b): + Q = self.Q + pr = (a[0] * b[0], a[1] * b[1], a[0] * b[1], a[1] * b[0]) + self.prod = max(self.prod, *(abs(v) for v in pr)) + re, im = pr[0] - pr[1], pr[2] + pr[3] + self.sum = max(self.sum, abs(re), abs(im)) + need(max(abs(re), abs(im)) < 1 << (2 * self.p + 7), "signed sums fit 2p+8 bits") + + def rnd(v): + return (v + Q // 2) >> self.p if v >= 0 else -((-v + Q // 2) >> self.p) + + out = (rnd(re), rnd(im)) + self.coord = max(self.coord, abs(out[0]), abs(out[1])) + return out + + +def phase_case(D, t, extra_bits, seed_choice, ref_digits): + R = 1 << t + d = 1 + while R ** d < D: + d += 1 + E = R ** d + p = (64 * (E + d)).bit_length() + extra_bits + Q = 1 << p + need(Q >= 64 * (E + d), "Q >= 64(E+d)") + pi = agm_pi(ref_digits) + decimal.getcontext().prec = ref_digits + 10 + c, s = dec_cos_sin(2 * pi / D, ref_digits) + fx, fy = int((c * Q).to_integral_value(rounding=decimal.ROUND_FLOOR)), int((s * Q).to_integral_value(rounding=decimal.ROUND_FLOOR)) + z0 = (fx + (seed_choice & 1), fy + (seed_choice >> 1 & 1)) + need(abs(Dec(z0[0]) / Q - c) < Dec(1) / Q and abs(Dec(z0[1]) / Q - s) < Dec(1) / Q, "seed coordinate error < 1/Q") + M = Mag(Q, p) + zj = z0 + tables = [] + for j in range(d): + if j: + for _ in range(t): + zj = M.mul(zj, zj) + row = [(Q, 0)] + for a in range(1, R): + row.append(M.mul(row[-1], zj)) + tables.append(row) + bound = Dec(16 * (E + d)) / Q + worst = Dec(0) + es = range(D) if D <= 4096 else sorted(set([0, 1, D - 1] + [random.Random(D + t).randrange(D) for _ in range(300)])) + for e in es: + acc, rest = (Q, 0), e + for j in range(d): + acc = M.mul(acc, tables[j][rest & (R - 1)]) + rest >>= t + ce, se = dec_cos_sin(2 * pi * e / D, ref_digits) if e else (Dec(1), Dec(0)) + err = ((Dec(acc[0]) / Q - ce) ** 2 + (Dec(acc[1]) / Q - se) ** 2).sqrt() + need(err <= bound, "query error <= 16(E+d)/Q") + worst = max(worst, err / bound) + need(M.coord <= 2 * Q and M.prod <= 4 * Q * Q and M.sum <= 8 * Q * Q, "magnitude envelope") + # The multiplicative argument: at most 2E+d factors of (1+eta), eta = 2/Q. + eta = Fr(2, Q) + x = (2 * E + d) * eta + need(x <= Fr(1, 16) and 2 * x <= Fr(16 * (E + d), Q), "(1+eta)^(2E+d)-1 <= 2(2E+d)eta <= 16(E+d)/Q") + return len(es), float(worst) + + +def check_phase(): + queries, worst = 0, 0.0 + for D, t in ((12, 1), (97, 2), (1000, 3), (4096, 4), (65537, 4), (1 << 20, 4)): + for choice in range(4): + q, w = phase_case(D, t, 0, choice, 40) + queries += q + worst = max(worst, w) + # Leaf/node count by literal unfolding of the computation tree. + for t in (1, 2, 3): + R = 1 << t + for d in (1, 2, 3): + leaves = {} + nodes = {} + for j in range(d): + zl, zn = R ** j, R ** j - 1 # t squarings per level: R^j leaves, R^j - 1 products + for a in range(R): + leaves[(j, a)] = a * zl + nodes[(j, a)] = a * (zn + 1) # a products by z_j (the first one by the exact 1) + E = R ** d + for e in range(E): + digits = [(e >> (t * j)) & (R - 1) for j in range(d)] + lv = sum(leaves[(j, a)] for j, a in enumerate(digits)) + nd = sum(nodes[(j, a)] for j, a in enumerate(digits)) + d + need(lv == e and nd <= e + d and lv + nd <= 2 * E + d, "at most 2E+d error factors") + return queries, worst + + +# ---------------------------------------------------------------- 2. exp table +def floor_S(v, S): + return Fr((v * (1 << S)).__floor__(), 1 << S) + + +def exp_case(n, c=1, C=2, eta=Fr(1, 4), d=13, c0=6, queries=200): + need(d % 2 == 1 and eta * (d + 1) > c + eta + 2, "odd d with eta(d+1) > c+eta+2") + lg = (n + 2).bit_length() - 1 # floor(log2(n+2)); p <= c log2(n+2) + p = c * lg + h = Fr(1, 1 << math.floor(eta * math.log2(n + 2))) + J = math.ceil(C * math.log2(n + 2) / h) + S = p + (J + 1 - 1).bit_length() + c0 # ceil(log2(J+1)) = bitlength(J) for J >= 1 + need((1 << (S - p - c0)) >= J + 1, "ceil(log2(J+1)) bits") + tail = h ** (d + 1) / math.factorial(d + 1) + tail_ok = tail <= Fr(1, 1 << S) + q0 = sum((-h) ** k / math.factorial(k) for k in range(d + 1)) + qhat = floor_S(q0, S) + digits = S // 3 + 30 + e_h = dec_exp(-h, digits) + decimal.getcontext().prec = digits + 10 + need(Dec(qhat.numerator) / qhat.denominator <= e_h and qhat >= 0, "0 <= qhat <= e^-h") + Cd = (e_h - Dec(qhat.numerator) / qhat.denominator) * (1 << S) + if tail_ok: + need(Cd <= 2, "C_d <= 2 when the tail is below 2^-S") + need(3 * J <= 1 << (S - p - 4), "(C_d+1) J 2^-S <= 2^-(p+4)") + Ehat = Fr(1) + ref = Dec(1) + worst_table = Dec(0) + table = [Ehat] + for j in range(1, J + 1): + Ehat = floor_S(Ehat * qhat, S) + ref = ref * e_h + gap = ref - Dec(Ehat.numerator) / Ehat.denominator + need(gap >= -Dec(10) ** (-digits + 5) and gap <= j * (Cd + 1) / (1 << S) + Dec(10) ** (-digits + 5), + "0 <= e^(-jh) - Ehat_j <= j(C_d+1)2^-S") + worst_table = max(worst_table, gap * (1 << S) / (j * (Cd + 1))) + table.append(Ehat) + rng = random.Random(n) + xmax = Fr(C * math.log2(n + 2)).limit_denominator(1 << 30) + frac_bits = p + 8 + xs = [Fr(0), h, 5 * h, Fr(int(xmax * (1 << frac_bits)), 1 << frac_bits)] + xs += [Fr(rng.randrange(int(xmax * (1 << frac_bits))), 1 << frac_bits) for _ in range(queries)] + worst_q = Dec(0) + for x in xs: + j = int(x / h) # shifts: h is a power of two + u = x - j * h + need(0 <= u < h and j <= J, "range reduction") + pu = floor_S(sum((-u) ** k / math.factorial(k) for k in range(d + 1)), S) + val = floor_S(table[j] * pu, S) + err = abs(Dec(val.numerator) / val.denominator - dec_exp(-x, digits)) + need(err < Dec(1) / (1 << p), "query error < 2^-p") + worst_q = max(worst_q, err * (1 << p)) + return {"n": n, "p": p, "J": J, "S": S, "tail_below_2^-S": tail_ok, "C_d": float(Cd), + "table_worst_over_bound": float(worst_table), "query_worst_over_2^-p": float(worst_q)} + + +def cutoff(c=1, C=2, eta=Fr(1, 4), d=13, c0=6): + """Smallest k such that the Taylor tail is below 2^-S for every n = 2^m - 2, k <= m <= 4096 (1: no cutoff needed).""" + last_bad = 0 + for m in range(2, 4097): + lg = m + h_exp = math.floor(eta * m) + J = math.ceil(C * m * (1 << h_exp)) + S = c * lg + J.bit_length() + c0 + # tail = 2^(-h_exp(d+1)) / (d+1)! <= 2^-S + if Fr(1, 1 << (h_exp * (d + 1))) / math.factorial(d + 1) > Fr(1, 1 << S): + last_bad = m + return last_bad + 1 + + +def check_reciprocal(n=1 << 12, C=2): + """exp(+x) = 1/exp(-x): request eps <= 2^(-p-3) e^(-2 C log2(n+2)), round the reciprocal within 2^(-p-2).""" + p = (n + 2).bit_length() - 1 + L = C * math.log2(n + 2) + k = math.ceil(2 * L * math.log2(math.e)) + p + 3 # 2^-k <= 2^(-p-3) e^(-2L) + digits = k // 3 + 40 + rng = random.Random(7) + worst = Dec(0) + for _ in range(100): + x = Fr(rng.randrange(int(L * 2 ** 20)), 2 ** 20) + decimal.getcontext().prec = digits + 10 + true_neg = dec_exp(-x, digits) + approx = true_neg + Dec(rng.uniform(-1, 1)) / Dec(2) ** k # any approximation within eps + need(approx >= Dec(-L).exp() / 2, "approximation stays above half the lower bound") + recip = (Dec(1) / approx * Dec(2) ** (p + 2)).to_integral_value() / Dec(2) ** (p + 2) + err = abs(recip - dec_exp(x, digits)) + need(err < Dec(1) / Dec(2) ** p, "exp(+x) error < 2^-p") + worst = max(worst, err * Dec(2) ** p) + return float(worst) + + +def check_binding(): + bound = 0 + for name, sha in SPECS.items(): + need(hashlib.sha256((HERE / name).read_bytes()).hexdigest() == sha, f"vendored {name}") + bound += 1 + ph = (HERE / "PR4554-PHASE-TABLES.md").read_text(encoding="utf-8") + for s in ("Q>=64(E+d)", "absolute error at most 16(E+d)/Q", "(1+eta)^(2E+d) - 1", "w>=2p+8 bits"): + need(s in ph, f"PHASE-TABLES.md states: {s}") + bound += 1 + ex = (HERE / "PR4554-NATIVE-EXP-TABLE.md").read_text(encoding="utf-8") + for s in ("eta(d+1)>c+eta+2", "0<=exp(-jh)-Ehat_j<=j(C_d+1)2^-S.", "S=p+ceil(log2(J+1))+c0", + "epsilon<=2^(-p-3) exp(-2C log2(n+2))"): + need(s in ex, f"NATIVE-EXP-TABLE.md states: {s}") + bound += 1 + return bound + + +def main(): + q, w = check_phase() + report = {"status": "PASS_PR4554_TABLES_CHECK", "phase_queries": q, "phase_worst_error_over_bound": w, + "exp_cases": [exp_case(n) for n in (1 << 16, 1 << 20, 1 << 24)], + "exp_tail_cutoff_log2n": cutoff(), "exp_reciprocal_worst_over_2^-p": check_reciprocal(), + "spec_statements_bound": check_binding(), "new_kappa_claimed": False} + print(json.dumps(report, indent=1)) + + +if __name__ == "__main__": + sys.exit(main())