diff --git a/research/pr4554-native-fft-check/.gitattributes b/research/pr4554-native-fft-check/.gitattributes new file mode 100644 index 0000000000..c20c6dc8cf --- /dev/null +++ b/research/pr4554-native-fft-check/.gitattributes @@ -0,0 +1,2 @@ +PR4554-PROOF.md -text +PR4554-FFT-ERROR.md -text diff --git a/research/pr4554-native-fft-check/PR4554-FFT-ERROR.md b/research/pr4554-native-fft-check/PR4554-FFT-ERROR.md new file mode 100644 index 0000000000..53691b10de --- /dev/null +++ b/research/pr4554-native-fft-check/PR4554-FFT-ERROR.md @@ -0,0 +1,159 @@ +# Independent numerical error bound for the native radix-two FFT + +Prepared by Dugongue with substantial OpenAI Codex assistance. Apache-2.0. +This argument is independent of the transform's arithmetic implementation. +It uses no TM computation or simulation. + +## Assumptions and norms + +Let N=2^h>=2 and let each input coefficient have complex modulus at most beta, +beta>=1. The exact forward and inverse butterflies are unnormalized, with +true twiddles of modulus one. The inverse is divided by N once at the end. +Let S=2^F. Suppose every supplied approximate twiddle differs from its true +value by at most rho=16N/S in each real/imaginary coordinate. Suppose each +complex multiplication's local fixed-grid rounding contributes at most 2/S +in each coordinate. Addition/subtraction at the common grid is exact, and +all actual word operations satisfy their no-overflow bounds. + +Use complex modulus for this proof, and the maximum such modulus over an array. +The twiddle error is then at most delta=32N/S; every local multiplication or +final normalization rounding can be bounded by eta=4/S. Assume delta<=1 and +eta<=1; the final parameters below easily imply both. Conjugation of twiddles +for the inverse preserves these error bounds. + +## One forward transform + +If A_j bounds the approximate data modulus after j stages, then + + A_(j+1) <= (2+delta) A_j + eta <= 4 A_j + eta. + +Since A_0<=beta, induction gives + + A_j <= 4^j(beta+eta/3) <= 2*4^j beta, + A_h <= 2N^2 beta. + +Let E_j be the error from the exact partial FFT. Subtract exact and approximate +butterflies, assigning the twiddle discrepancy to the approximate second input: + + E_(j+1) <= 2 E_j + delta A_j + eta, E_0=0. + +The elementary sum + + sum_(j=0)^(h-1) 2^(h-1-j) 4^j = N(N-1)/2 + +therefore gives + + E_h <= delta beta N(N-1) + eta(N-1) + <= 36N^3 beta/S. (1) + +The exact forward data modulus is at most N beta. The same bounds hold for the +other operand's independent forward transform. + +## Pointwise product + +Let A,B be the exact forward arrays and Ahat,Bhat the computed arrays. A rounded +pointwise product has error + + |round(Ahat*Bhat)-A*B| + <= |Ahat-A| |Bhat| + |A| |Bhat-B| + eta + <= (36N^3 beta/S)(2N^2 beta + N beta) + eta + <= 112N^5 beta^2/S. (2) + +The computed product modulus is at most + + (2N^2 beta)^2 + eta <= 5N^4 beta^2. (3) + +No independence, random error cancellation or favorable data signs are used. + +## Inverse transform and normalization + +Applying the same amplitude recurrence to the approximate product input gives +an inverse-stage amplitude bound + + A'_j <= 6*4^j N^4 beta^2. + +The inverse error recurrence now starts with (2), rather than zero. Thus + + E'_h <= N*(112N^5 beta^2/S) + + delta*6N^4 beta^2*N(N-1)/2 + eta(N-1) + <= 212N^7 beta^2/S. + +Division by the exact integer N scales this preexisting error by 1/N. Its own +fixed-grid coordinate rounding contributes at most eta in complex modulus. +Consequently every normalized output coefficient has error at most + + 216N^6 beta^2/S. (4) + +In particular, the more conservative requested envelope is valid: + + error <= 2^28 N^8 (F+1)^2 beta^4/S. (5) + +Indeed N>=1, beta>=1 and F+1>=1, while 216<2^28. The factors (F+1)^2 are not +needed for this independent estimate; retaining them safely accommodates the +stated caller ledger. They must not stand in for an unproved twiddle or local +rounding premise. + +## Parameters imply strict recovery error below 1/4 + +Let n>=1, + + ell = bitlength(n+2), + b = max(1,floor(ell/8)), + beta = 2^b, + N = least power of two >= 2ceil(n/b), + F = 24ell+64, S=2^F. + +Then ell>=2, b<=ell, and + + N < 4ceil(n/b) <= 4n < 2^(ell+2), + beta^4 <= 2^(4ell). + +The first inequality is the ordinary bound on rounding a positive integer up +to a power of two. Also + + F+1=24ell+65 <= 2^(ell+7). + +For ell=2 the inequality is immediate; increasing ell by one adds only 24 to +the left while doubling the right, and the previous left side is already at +least 24. Substitute these bounds in (5): + + 2^28 N^8(F+1)^2 beta^4/S + < 2^[28+8(ell+2)+2(ell+7)+4ell-(24ell+64)] + = 2^(-10ell-6) + <= 2^-26 < 1/4. + +The same parameter inequalities imply delta=32N/S<1 and eta=4/S<1, discharging +the earlier stability conditions. The empty-input case n=0 should be returned +as zero directly, rather than introducing an unnecessary FFT convention. + +For zero-padded nonnegative digit arrays, N>=2ceil(n/b) prevents cyclic wrap +in the linear convolution. The exact normalized inverse coefficients are real +integers. Equation (5) therefore places each computed real coordinate strictly +within 1/4 of the right integer, and nearest-integer recovery is exact. The +linear carry theorem in ROOTS-RECOVERY.md then produces the exact product. + +## Word-size sanity bound + +The approximate inverse data modulus before normalization is at most +6N^6 beta^2. On the F-bit grid its signed integer coordinates have magnitude +at most 6S N^6 beta^2. Twiddle coordinates have magnitude at most 2S. A raw +real multiplication used in a complex multiply is consequently bounded by +12S^2 N^6 beta^2, and a constant additional factor covers combining two such +products. Thus a sufficient magnitude width is + + 2F + 6ceil(log2 N) + 2b + 8, + +plus caller addresses and a sign representation. This is O(log n) for the +chosen F. The independently generated root seed uses an even larger fixed +constant times F for its internal precision; one common global word width +must cover both stages. Allowing each subroutine an unaccounted private wider +word would not establish a single Word RAM implementation. + +## Scope + +This establishes numerical sufficiency for an actually implemented iterative +FFT whose twiddles, local rounding and word bounds meet the hypotheses. It does +not verify an arbitrary FFT routine by observing good samples, or prove that +an unrelated power-saving transform inherits these polynomial prefix norms. +The native implementation and its exact operation ledger remain separately +reviewable. No improved TM exponent or formal Lean proof is claimed. diff --git a/research/pr4554-native-fft-check/PR4554-PROOF.md b/research/pr4554-native-fft-check/PR4554-PROOF.md new file mode 100644 index 0000000000..c7b4777b05 --- /dev/null +++ b/research/pr4554-native-fft-check/PR4554-PROOF.md @@ -0,0 +1,113 @@ +# Native bounded-word integer multiplication + +Prepared by Dugongue with substantial OpenAI Codex assistance. Apache-2.0. + +## The theorem and its model + +For every n ≥ 1 and two nonnegative integers of at most n bits, the algorithm in +native_fft.py computes their exact packed product on a uniform Word RAM in O(n) +word operations and O(n/log(n+2)) words of memory. The program uses add/subtract, +comparison, shifts, masks, indexed memory and multiplication/division of bounded +words. Its word width is a sufficiently large fixed constant times log(n+2). +The input/output consists of packed radix digits; converting conventional packed +binary words to/from these digits takes O(n/log(n+2)) operations. + +This is an informal algorithmic theorem accompanied by a portable executable +bounded-word model, not an emitted machine program or a Lean theorem. Linear +Word RAM FFT multiplication is classical: see Fürer, 2014, Theorem 4.1(a), +https://arxiv.org/abs/1402.1811. Our contribution is a self-contained, paid, +source-bound native implementation and the explicit precision, error and assembly +interfaces needed to remove the prior TM multiplication subroutine. O(n) is not +claimed as a new record, or as a higher κ. + +## Fully specified construction + +Set ell = bitlength(n+2), b = max(1,floor(ell/8)), β = 2^b, +L = ceil(n/b), N the least power of two at least 2L, F = 24ell+64 and S = 2^F. +Use a single global magnitude-word width w = 16(F+16)+20; signed values have a +separate sign bit. These are fixed formulas, with w = Θ(log(n+2)). The large +constant simplifies this proof; it is not a claim of practical speed. + +The packed ABI supplies L radix-β digits per operand, with unused high input +bits zero. Pad each to N. Construct exp(2πi/N) to coordinate error < 1/S using +root_seed.py: Machin's formula and bounded Taylor arithmetic, not a root oracle. +ROOTS-RECOVERY.md proves its universal precision, word bound and O(F) cost. +Generate the first N/2 twiddles by successive rounded complex products. + +Perform two ordinary iterative radix-two forward FFTs, multiply their entries, +perform the unnormalized inverse FFT, then divide by N. Each real scalar product +is taken exactly within the common word and the complex product is rounded back +to the S grid. Forward twiddles are conjugates of the positive-angle table; +inverse twiddles use that table directly. Additions are exact at the common grid. +Round the real inverse coefficients to integers and carry in radix β. + +## Discharging numerical premises + +Let δ = 1/S. A seed's complex error is < 2δ, and one rounded complex product +has error < 2δ. If a_j is the error in the j-th stored power, then + + a_(j+1) ≤ (1+2δ)a_j + 4δ, a_0 = 0. + +Since Nδ < 1/16, summing this geometric recurrence gives a_j ≤ 8jδ ≤ 16Nδ. +Every twiddle therefore has coordinate error at most 16N/S, including the +conjugated table. Local fixed-grid product errors are less than 2/S per +coordinate; floor division instead of nearest rounding still satisfies this +bound. The independent proof in FFT-ERROR.md now applies. It proves final +complex error at most 216N^6β²/S, below the looser envelope + + 2^28 N^8(F+1)^2β^4/S < 2^(-10ell-6) ≤ 2^-26 < 1/4. + +The separate integer division by N has its own grid-rounding error included in +that estimate. Thus every rounded real coefficient is the correct integer; +the imaginary coordinate is within 1/4 of zero. The runtime separation tests +always succeed for valid input. They are checks of a proved property, rather +than a substitute for a closeness-to-the-right-integer proof. + +N ≥ 2L prevents cyclic aliasing. Each coefficient lies between 0 and +L(β−1)². Carry propagation uses s_j = c_j + carry_j, output digit s_j mod β +and carry_(j+1) = floor(s_j/β). Inductively carry_j ≤ L(β−1), so every carry +operation has O(log n) bits. After 2L digits the carry is zero. Product size +< 2^(2n) proves that trimmed output digits and unused high bits are zero. +This includes n = 1 and all finite small cases without a separate FFT convention. + +## One common bounded-word model + +FFT-ERROR.md bounds inverse intermediate complex magnitudes by 6N^6β². Raw +scalar product magnitudes, including sums of two such products, need at most +2F + 6ceil(log2 N) + 2b + 8 bits. The root routine needs at most +16(F+16)+20 bits, which dominates this FFT bound for every ell ≥ 2 since +ceil(log2 N) < ell+2 and b ≤ ell. Address, loop and carry words fit too. +The runtime rejects any stored arithmetic magnitude reaching 2^w. +Its bound-checking interpreter uses arbitrary-precision Python integers only +to emulate those bounded words. No runtime operation multiplies the n-bit +operands as integers. Python's independent x*y appears only in the test oracle. +The physical word constant is chosen before execution from the bit budget n; +subroutines do not silently switch to larger private words. + +## Complete cost + +Root setup costs O(F) operations. Its internal precision 8(F+16) and its +128(F+1) Taylor iterations remain O(log n). Twiddle setup costs O(N) rounded +complex word products. Bit reversal costs O(N): the number of inner bit flips +summed over the whole index loop is less than 2N. There are log2 N stages, +N/2 butterflies per stage, with a constant number of scalar operations, address +updates, loads and stores per butterfly. Three transforms cost O(N log N). +Pointwise products, validation of padding, rounding, carries, allocation, +packing and final output cost O(N). The counted arithmetic ledger is exposed +by the executable; it is not a literal emitted-instruction trace. + +For ell ≥ 16, b ≥ ell/16, so N < 4ceil(n/b) ≤ 64n/ell + 4 and log2 N ≤ ell+2. +Consequently N log N = O(n). The finite remaining lengths are absorbed in a +fixed constant. All setup costs are paid; no input-independent preprocessing +is excluded. Storage is O(N) words, plus O(1) root-generation words. + +## Use in the faster route + +This construction is a native implementation for every record multiplication +in a future Fourier/resampling compiler. NATIVE-OUTER-ASSEMBLY.md proves the +quantitative caller and full multiplication cost if a faster cyclic transform +is actually supplied. NATIVE-STABILITY.md proves a bounded-word rational +Gaussian kernel and safe shear divisions; PHASE-TABLES.md proves paid native +phase queries. Those positive results do not silently supply the missing cyclic +conditioning bound. No existing TM instruction table, TM simulation or TM +multiplier is used by this construction. diff --git a/research/pr4554-native-fft-check/README.md b/research/pr4554-native-fft-check/README.md new file mode 100644 index 0000000000..d0d141ffdb --- /dev/null +++ b/research/pr4554-native-fft-check/README.md @@ -0,0 +1,55 @@ +# Independent check of #4554: native FFT multiplier error chain and parameters + +**Task:** independent-moment (independent source-bound coverage). **Base:** canonical main `f626db5c`. **No new κ.** + +Dugongue's #4554 (head `4f46c6a9`) gives a native bounded-word FFT multiplier. `PROOF.md` sets its parameters, and `FFT-ERROR.md` proves a recovery error of at most `216N⁶β²/S < 2^(−10ℓ−6) < 1/4`. This package checks both. They are vendored as `PR4554-PROOF.md` (blob `c7b4777b`) and `PR4554-FFT-ERROR.md` (blob `53691b10`), with their SHA256 pinned. Nothing is imported from #4554's `native_fft.py`, `root_seed.py` or `recover.py`. It follows my #4556 on the same PR. + +## Results + +All checked claims are confirmed. + +**A. Error chain, in exact rational arithmetic.** The recurrences are run with equality, as the worst case, for `N = 2^h` (`h ≤ 20`) and `β = 2^b` (`b ≤ 8`), at `S = 2^F` with `δ = 32N/S` and `η = 4/S`. That is 80 cases. +- **Forward amplitudes.** They stay below `4^j(β+η/3) ≤ 2·4^jβ`. +- **Forward errors.** They stay below `δβN(N−1) + η(N−1) ≤ 36N³β/S`, which is (1). +- **Pointwise products.** The error is at most `112N⁵β²/S` (2), and the computed modulus is at most `5N⁴β²` (3). +- **Inverse transform.** Amplitudes stay at most `6·4^jN⁴β²`, and the final error before normalization is at most `212N⁷β²/S`. +- **Normalized output.** The error is at most `216N⁶β²/S` (4), which is at most the envelope (5). The worst recurrence value is 0.196 of (4), as expected for a bound proved with loose amplitudes. +- **Twiddles.** The recurrence `a_(j+1) = (1+2δ)a_j + 4δ` stays at most `8jδ` (3,000 steps). Analytically, `2((1+x)^j − 1) ≤ 2(e^(jx) − 1) ≤ 4jx` when `jx ≤ 1`. +- **Parameters, exact integers, at the largest `n` for each `ℓ = 2…200`.** + - `N < 4L ≤ 4n < 2^(ℓ+2)` and `F+1 ≤ 2^(ℓ+7)`; + - `2^28 N^8 (F+1)^2 β^4 · 2^(10ℓ+6) < S`; + - `δ, η < 1`; + - the raw-product budget is within `2F + 6⌈log₂N⌉ + 2b + 8 ≤ 16(F+16)+20`. + +**B. A separately written fixed-grid multiplier** following the `PROOF.md` construction. Its steps are: +- **Root.** π by Machin, cos and sin by Taylor at `F+64` bits, then floored to the `2^−F` grid. +- **Transforms.** Twiddles come from successive floor-rounded complex products, and forward twiddles are conjugated. Two forward FFTs are followed by rounded pointwise products, an unnormalized inverse and floor division by `N`. +- **Recovery.** Each coefficient is rounded to the nearest integer, and carries propagate in radix `β`. + +It ran 367 products: all pairs for `n ≤ 4`, plus random and all-ones operands up to `n = 12,000` bits. All are exact, checked against Python's multiplication, which is used as the oracle only. +- **Observed error.** The worst observed coefficient error is **5.7·10⁻⁷ of (4)**, always `< 1/4`. +- **Magnitudes.** Inverse magnitudes stay below `6N⁶β²·S`. Carries stay at most `L(β−1)`, and the final carry is zero. +- **Word width.** Every stored magnitude is below `2^w` with `w = 16(F+16)+20`. The FFT itself uses at most 12% of that width, since the width is set by the root routine, as the PR says. + +**C. Binding and replay.** 12 statements are bound verbatim to the vendored texts. #4554's registered package replays at its head (`outcome=completed`, 16.6 s; see #4556). + +## Verify + +```text +python3 -B research/pr4554-native-fft-check/verify.py +``` + +This runs in about 16 s with the standard library only. + +## Limitations + +- **Root precision.** `ROOTS-RECOVERY.md`'s universal root-precision proof is not re-derived here. My root is computed independently and is only required to meet the same `< 2/S` seed premise. +- **Operation ledger.** The operation-count ledger and `O(n)` cost are not re-counted. +- **Other components.** The Gaussian adapters, phase tables and outer assembly are out of scope; #4556 covers the banded inverse. +- **Finite checks.** Part B is a finite family. The universal bound rests on the written recurrences that part A checks. + +## Open next questions + +- An independent check of `ROOTS-RECOVERY.md`, covering root-seed precision with the Machin and Taylor truncation and the `O(F)` cost. First step: certified truncation bounds for both series at precision `8(F+16)`. + +Credit: Dugongue (#4554); Fürer (2014) for linear Word RAM FFT multiplication. Prepared by Ryanson Jonathan with substantial Claude Code (Claude Opus 5.5) assistance. diff --git a/research/pr4554-native-fft-check/SOURCES.json b/research/pr4554-native-fft-check/SOURCES.json new file mode 100644 index 0000000000..a2a5ea8bbd --- /dev/null +++ b/research/pr4554-native-fft-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-fft-check/TASK.json b/research/pr4554-native-fft-check/TASK.json new file mode 100644 index 0000000000..ff0e838203 --- /dev/null +++ b/research/pr4554-native-fft-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 native multiplier that removes the TM multiplication subroutine from the Word RAM route.", + "budget_minutes": 45, + "stop_rule": "Stop at PROOF.md + FFT-ERROR.md error chain and parameters.", + "outcome": "completed", + "evidence": "Exact rational worst-case recurrences for (1)-(5) over 80 (N,beta) cases; twiddle recurrence; exact parameter inequalities for ell=2..200 incl. word budget; separately written fixed-grid FFT multiplier exact on 367 products up to 12,000 bits with observed error 5.7e-7 of the proved bound; 12 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 #4554 ROOTS-RECOVERY.md (root-seed precision and cost).", + "why_it_matters": "The FFT bound assumes the <2/S seed premise.", + "first_step": "Certified Machin/Taylor truncation bounds at precision 8(F+16).", + "blocked_by": "None." + } + ] +} diff --git a/research/pr4554-native-fft-check/VERIFY.json b/research/pr4554-native-fft-check/VERIFY.json new file mode 100644 index 0000000000..99b837b966 --- /dev/null +++ b/research/pr4554-native-fft-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": "Exact error-chain recurrences (1)-(5), twiddle recurrence, parameter inequalities ell<=200; separate fixed-grid FFT multiplier exact on 367 products with observed error vs proved bound; binding to vendored PROOF.md/FFT-ERROR.md (sha256)." + } + ] +} diff --git a/research/pr4554-native-fft-check/verify.py b/research/pr4554-native-fft-check/verify.py new file mode 100644 index 0000000000..4bd32f48c7 --- /dev/null +++ b/research/pr4554-native-fft-check/verify.py @@ -0,0 +1,267 @@ +#!/usr/bin/env python3 +"""Independent check of PR #4554 (Dugongue) native FFT multiplier: PROOF.md and FFT-ERROR.md, from vendored text only. + +Nothing is imported from #4554's native_fft.py, root_seed.py or recover.py. + A. The error chain, run as worst-case recurrences in exact rational arithmetic for N = 2^h (h <= 20) and + beta = 2^b (b <= 8), at the stated S = 2^F: forward amplitudes A_(j+1) = (2+delta)A_j + eta stay below + 4^j(beta+eta/3) <= 2*4^j beta; errors E_(j+1) = 2E_j + delta*A_j + eta stay below the closed form + delta*beta*N(N-1) + eta(N-1) <= 36N^3 beta/S (1); the pointwise error is <= 112N^5 beta^2/S (2) with product + modulus <= 5N^4 beta^2 (3); inverse amplitudes stay <= 6*4^j N^4 beta^2; E'_h <= 212N^7 beta^2/S; the normalized + error is <= 216N^6 beta^2/S (4) <= 2^28 N^8 (F+1)^2 beta^4/S (5). The twiddle recurrence + a_(j+1) = (1+2d)a_j + 4d gives a_j <= 8jd when 2jd <= 1. The parameter step (5) < 2^(-10 ell - 6) <= 2^-26 and + delta, eta < 1 are checked exactly at the largest n of every ell in 2..200, with the word budget + 2F + 6ceil(log2 N) + 2b + 8 <= 16(F+16)+20. + B. A separately written fixed-grid FFT multiplier following the PROOF.md construction: own root of unity (pi by + Machin, cos/sin by Taylor, floor to the 2^-F grid), twiddles by successive floor-rounded complex products, + conjugated forward twiddles, two forward transforms, rounded pointwise products, unnormalized inverse, floor + division by N, nearest-integer recovery and radix-beta carries. Products are exact against Python's + multiplication (oracle only) on random, all-ones and small exhaustive operands; the observed coefficient error + and intermediate magnitudes are compared with (4), 1/4 and the 6N^6 beta^2 amplitude bound, and every stored + value with the global word width. +Ryanson Jonathan with Claude Code (Claude Opus 5.5) assistance. Apache-2.0. +""" +from fractions import Fraction as Fr +import hashlib +import json +from pathlib import Path +import random +import sys + +HERE = Path(__file__).resolve().parent +SPECS = {"PR4554-PROOF.md": "db7d6ceb495e4b111764bdea1905f66c7e6f0921bd8c1e3048054445ad448c41", "PR4554-FFT-ERROR.md": "ca426930df5c449ca8ce432360ba3ac80f12377bcd9a058ad16387d8ad87bef5"} # #4554 head 4f46c6a9, blobs c7b4777b (PROOF), 53691b10 (FFT-ERROR) + + +def need(condition, message): + if not condition: + raise AssertionError(message) + + +def params(n): + ell = (n + 2).bit_length() + b = max(1, ell // 8) + L = -(-n // b) + N = 1 + while N < 2 * L: + N *= 2 + F = 24 * ell + 64 + return ell, b, L, N, F + + +# ---------------------------------------------------------------- A: the error chain +def chain(h, b, F): + N, beta, S = 1 << h, Fr(1 << b), Fr(1 << F) + delta, eta = 32 * N / S, 4 / S + A, E = beta, Fr(0) + for j in range(h): + need(A <= 4 ** j * (beta + eta / 3) <= 2 * 4 ** j * beta, "forward amplitude") + A, E = (2 + delta) * A + eta, 2 * E + delta * A + eta + need(A <= 2 * N * N * beta, "A_h <= 2N^2 beta") + need(E <= delta * beta * N * (N - 1) + eta * (N - 1) <= 36 * N ** 3 * beta / S, "(1)") + P = E * (2 * N * N * beta + N * beta) + eta + need(P <= 112 * N ** 5 * beta ** 2 / S, "(2)") + need((2 * N * N * beta) ** 2 + eta <= 5 * N ** 4 * beta ** 2, "(3)") + A2, E2 = (2 * N * N * beta) ** 2 + eta, P + for j in range(h): + need(A2 <= 6 * 4 ** j * N ** 4 * beta ** 2, "inverse amplitude") + A2, E2 = (2 + delta) * A2 + eta, 2 * E2 + delta * A2 + eta + need(A2 <= 6 * N ** 6 * beta ** 2, "inverse data modulus") + need(E2 <= 212 * N ** 7 * beta ** 2 / S, "E'_h") + final = E2 / N + eta + need(final <= 216 * N ** 6 * beta ** 2 / S <= 2 ** 28 * N ** 8 * (F + 1) ** 2 * beta ** 4 / S, "(4) and (5)") + return final * S / (216 * N ** 6 * beta ** 2) + + +def check_chain(): + cases, worst = 0, Fr(0) + for h in range(1, 21): + for b in (1, 2, 4, 8): + worst = max(worst, chain(h, b, 24 * max(2, h) + 64)) + cases += 1 + tw = 0 + S = Fr(1 << 20) + d = 1 / S + a = Fr(0) + for j in range(1, 3001): # 2jd <= 1 throughout + a = (1 + 2 * d) * a + 4 * d + need(a <= 8 * j * d, "twiddle recurrence") + tw += 1 + params_checked = 0 + for ell in range(2, 201): + n = (1 << ell) - 3 # largest n with this ell; N is nondecreasing in n + e2, b, L, N, F = params(n) + need(e2 == ell, "ell") + S = 1 << F + need(N < 4 * L <= 4 * n < (1 << (ell + 2)) and F + 1 <= 1 << (ell + 7), "parameter bounds") + need(2 ** 28 * N ** 8 * (F + 1) ** 2 * 2 ** (4 * b) * 2 ** (10 * ell + 6) < S, "(5) < 2^(-10 ell - 6)") + need(32 * N < S and 4 < S, "delta, eta < 1") + logN = N.bit_length() - 1 + need(2 * 12 * N ** 6 * 2 ** (2 * b) * S * S < 1 << (2 * F + 6 * logN + 2 * b + 8), "raw product budget") + need(2 * F + 6 * logN + 2 * b + 8 <= 16 * (F + 16) + 20, "root width dominates") + params_checked += 1 + return cases, float(worst), tw, params_checked + + +# ---------------------------------------------------------------- B: an independent fixed-grid multiplier +def fixed_pi(bits): + one = 1 << bits + + def arctan_inv(x): + total, term, k, sign = 0, one // x, 1, 1 + x2 = x * x + while term: + total += sign * (term // k) + term //= x2 + k += 2 + sign = -sign + return total + + return 4 * (4 * arctan_inv(5) - arctan_inv(239)) + + +def root(N, F): + """exp(2 pi i/N) floored to the 2^-F grid, from a Taylor series at F+64 bits.""" + G = F + 64 + one = 1 << G + x = 2 * fixed_pi(G) // N + c, s, term, k = 0, 0, one, 0 + while term: + 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 // (one * k) + return c >> 64, s >> 64 + + +class Word: + def __init__(self, w): + self.limit = 1 << w + self.peak = 0 + + def __call__(self, v): + a = abs(v) + need(a < self.limit, "stored magnitude below 2^w") + self.peak = max(self.peak, a) + return v + + +def cmul(a, b, S, W): + re = W(W(a[0] * b[0]) - W(a[1] * b[1])) + im = W(W(a[0] * b[1]) + W(a[1] * b[0])) + return (W(re // S), W(im // S)) + + +def fft(data, tw, N, S, W, forward): + a = list(data) + j = 0 + for i in range(1, N): # bit reversal + 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: + step = N // m + for start in range(0, N, m): + for k in range(m // 2): + w = tw[k * step] + if forward: + w = (w[0], -w[1]) + t = cmul(w, a[start + k + m // 2], S, W) + u = a[start + k] + a[start + k] = (W(u[0] + t[0]), W(u[1] + t[1])) + a[start + k + m // 2] = (W(u[0] - t[0]), W(u[1] - t[1])) + m *= 2 + return a + + +def multiply(x, y, n, stats): + ell, b, L, N, F = params(n) + S, beta = 1 << F, 1 << b + W = Word(16 * (F + 16) + 20) + w1 = root(N, F) + tw = [(S, 0)] + for _ in range(N // 2 - 1): + tw.append(cmul(tw[-1], w1, S, W)) + dx = [((x >> (b * i)) & (beta - 1)) * S for i in range(L)] + [0] * (N - L) + dy = [((y >> (b * i)) & (beta - 1)) * S for i in range(L)] + [0] * (N - L) + X = fft([(v, 0) for v in dx], tw, N, S, W, True) + Y = fft([(v, 0) for v in dy], tw, N, S, W, True) + Z = fft([cmul(p, q, S, W) for p, q in zip(X, Y)], tw, N, S, W, False) + inv_peak = max(max(abs(p[0]), abs(p[1])) for p in Z) + coeffs = [(p[0] // N, p[1] // N) for p in Z] + # Oracle only: exact convolution coefficients by Kronecker substitution. + slot = 2 * b + L.bit_length() + 2 + px = sum(((x >> (b * i)) & (beta - 1)) << (slot * i) for i in range(L)) + py = sum(((y >> (b * i)) & (beta - 1)) << (slot * i) for i in range(L)) + pz = px * py + exact = [(pz >> (slot * k)) & ((1 << slot) - 1) for k in range(N)] + err = max(max(abs(Fr(c[0], S) - e), abs(Fr(c[1], S))) for c, e in zip(coeffs, exact)) + need(err < Fr(1, 4), "recovery margin") + bound = Fr(216 * N ** 6 * beta ** 2, S) + need(err <= bound, "observed error within (4)") + need(inv_peak <= 6 * N ** 6 * beta ** 2 * S, "inverse amplitude bound") + digits = [(2 * c[0] + S) // (2 * S) for c in coeffs] # nearest integer + carry, out = 0, 0 + for k in range(2 * L + 1): + s_ = (digits[k] if k < N else 0) + carry + out |= (s_ % beta) << (b * k) + carry = s_ // beta + need(carry <= L * (beta - 1), "carry bound") + need(carry == 0, "final carry zero") + stats["worst_error_over_bound"] = max(stats["worst_error_over_bound"], float(err / bound)) + stats["peak_bits_over_width"] = max(stats["peak_bits_over_width"], W.peak.bit_length() / (16 * (F + 16) + 20)) + return out + + +def check_multiplier(): + rng = random.Random(45541) + stats = {"products": 0, "worst_error_over_bound": 0.0, "peak_bits_over_width": 0.0, "max_n": 0} + cases = [(n, x, y) for n in range(1, 5) for x in range(1 << n) for y in range(1 << n)] + for n in (7, 16, 31, 64, 100, 257, 1000, 4096, 12000): + cases.append((n, (1 << n) - 1, (1 << n) - 1)) + for _ in range(2): + cases.append((n, rng.getrandbits(n), rng.getrandbits(n))) + for n, x, y in cases: + need(multiply(x, y, n, stats) == x * y, "exact product") + stats["products"] += 1 + stats["max_n"] = max(stats["max_n"], n) + return stats + + +def check_binding(): + bound = 0 + for name, sha in SPECS.items(): + raw = (HERE / name).read_bytes() + need(hashlib.sha256(raw).hexdigest() == sha, f"vendored {name}") + bound += 1 + text = (HERE / "PR4554-FFT-ERROR.md").read_text(encoding="utf-8") + for s in ("E_(j+1) <= 2 E_j + delta A_j + eta,", "<= 36N^3 beta/S.", "<= 112N^5 beta^2/S.", "<= 212N^7 beta^2/S.", + "216N^6 beta^2/S.", "= 2^(-10ell-6)", "2F + 6ceil(log2 N) + 2b + 8,"): + need(s in text, f"FFT-ERROR.md states: {s}") + bound += 1 + proof = (HERE / "PR4554-PROOF.md").read_text(encoding="utf-8") + for s in ("a_(j+1) ≤ (1+2δ)a_j + 4δ,", "F = 24ell+64", "w = 16(F+16)+20"): + need(s in proof, f"PROOF.md states: {s}") + bound += 1 + return bound + + +def main(): + cases, worst, tw, pc = check_chain() + report = {"status": "PASS_PR4554_NATIVE_FFT_CHECK", "chain_cases": cases, "chain_worst_final_over_216_bound": worst, + "twiddle_steps": tw, "ell_values_checked": pc, "multiplier": check_multiplier(), + "spec_statements_bound": check_binding(), "new_kappa_claimed": False} + print(json.dumps(report, indent=1)) + + +if __name__ == "__main__": + sys.exit(main())