Independently check #4554's native FFT multiplier: exact error-chain recurrences, parameters to ell=200, separate fixed-grid multiplier (no new κ)

research/pr4554-native-fft-check/.gitattributes

@@ -0,0 +1,2 @@
+PR4554-PROOF.md -text
+PR4554-FFT-ERROR.md -text

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.

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.

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.

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"
+    }
+  ]
+}