Native Word RAM Gaussian adapters: chirp blocks and reusable banded inverse

research/native-word-fft/.gitattributes

@@ -0,0 +1 @@
+* -text

research/native-word-fft/CONTENT-HASHES.json

@@ -0,0 +1,38 @@
+{
+  ".gitattributes": "8b9c3028d6467ff3be7c793141927841d70f47f52d0590677438f488ca154992",
+  "band_solve.py": "f4e4471679cc89bbb01ab0b26e807b0dacdc6188ab95ba924fc57243369407bb",
+  "chirp_check.py": "3745bc34e548103434ee6442fe5287a830cc89b0cfb26591027c75cc086ecdad",
+  "EXP-IMPLEMENTATION-REVIEW.md": "6100ef60b1e5f53e63c44958772b4d4027026914e474007199e05052d492716d",
+  "exp_table.py": "e617b40d499f423a7861b4483bdacdb7db28232fadf7468efe871adae557fa35",
+  "FFT-ERROR.md": "ca426930df5c449ca8ce432360ba3ac80f12377bcd9a058ad16387d8ad87bef5",
+  "FINAL-REVIEW.md": "beddfb16d94c383b85b0db061a1eb683f4d68eff5932901f5b6bce78aca0d6e9",
+  "GAUSSIAN-CHIRP.md": "d8f2c45a985fc51d9523c3261ec39fc622afe737bf26a868e7f601dc7525b075",
+  "gaussian_block.py": "b58397fee1d9695400f3f5afd033d5f9aaae9b4318fa29c1ced4f618e20c1363",
+  "gaussian_inverse.py": "08616d26652c7c656ed8928a2c47a0a150a8cfb35bcc7f121ce5a729908e4a35",
+  "INDEPENDENT-REVIEW.md": "d30ddb3a17b5e0327b49bccd2bfd612619597b0113d4b5a00f8866333de34989",
+  "KERNEL-PRECEDENTS.json": "d5af5e07cd1d1d023864d3ed31189653d6d96fca38e2618d8e007e4c16dc2aee",
+  "LICENSE": "c71d239df91726fc519c6eb72d318ec65820627232b2f796219e87dcf35d0ab4",
+  "NATIVE-BANDED-INVERSE.md": "1259b4242ebf6ac7f14bfe72e371afb8dc530b2755d1fd82228d9dc6f01634ba",
+  "NATIVE-EXP-TABLE.md": "c9be6814a81e0bfe6b041f1a259a32369f7c030ca22eaa5c71c15911c888c6ff",
+  "NATIVE-OUTER-ASSEMBLY.md": "f1e194c01fc1ba4e33c698dbc287aeca8ff0efe822b195c50c95b619242cf9dc",
+  "NATIVE-RECOUNT.md": "fc14d31113ef3abd597e654e0aec35ef7972d014268e4b8b7624f202f9ff79db",
+  "NATIVE-STABILITY.md": "260e78bf8f88667a637b804d013d1e7c6dfdf1e23ae35ffdc64a4366c7321221",
+  "native_convolution.py": "a62a4cb5a1da458e31c90ea881a1146b59e9fd7ccbf317c1c9836c19bc7fa2c6",
+  "native_fft.py": "83bff7591e4ffa36fe3713f9d639de6ed1edc715d3698bc6d2d0dc3b90146cc5",
+  "NOTICE.md": "9e0b4a18fdd6b062c500a0123cfd269f4f0b1b02b6d24e39b786cf045489400a",
+  "OVERLAP.md": "a17118cf3d608771d0d38168b0e4712ffa4dd07caf499abbcdcd69425db9de0f",
+  "periodic_solve.py": "d3940b0d48ba4bbf81a0f2d05bdac29e3b998ad746771dcf01b03cac03680393",
+  "PHASE-TABLES.md": "9a5844176b67f81ad8468b77d0049935d05517a80e80586edfa8401f400917f8",
+  "phase_tables.py": "4d70323775a8efc3af82eff7ee567482bfa5d14d8c86289fe21c29c0d37422e1",
+  "PROOF.md": "db7d6ceb495e4b111764bdea1905f66c7e6f0921bd8c1e3048054445ad448c41",
+  "README.md": "7574dbc6ae49d79a96ecfdc3049a17a7e2002e2e02d5248d2967e5ddd89d923d",
+  "recover.py": "e05e9baf6709560c1498b094826a1cc9e4e84c1133226ea36739f2ae37301f52",
+  "RECOVERY-PRECEDENTS.json": "c06c4abc2c3e18776de7bb9fb400983b2ec1df8af8d0c108927ac3b7384d7c57",
+  "root_seed.py": "9b323f84ade7decfa6152dd4385e90ec43b56087927ac23ad4ef73ae4d5d3042",
+  "ROOTS-RECOVERY.md": "5e115c3ded1734dbcb119c889755a13ce9e9357598ec20a3eb1d3cdd727638f4",
+  "SOURCES.json": "c7050d7b52c6c262c6172dc1f01e9a249a4d104ca84c1c2bf780fa1cbb080630",
+  "TASK.json": "9a7e6849260cd0a6f1f49f9b27cbb9510ae49af38f0b6051ff26408bbc1f8dc0",
+  "VALIDATION.md": "427ba13876056552bd6a8ee30a7f22a79aa4fd6543b2c7cfddbbec98a0f5de8c",
+  "VERIFY.json": "8d54792bc63ff42d4a9c913deb3c80d7a8f1fc357216098b16f2a84d83dda304",
+  "verify.py": "86e9d740d1b01a8e46c4f4b6198d33a3ca6c4982d5c823352030573a8350138c"
+}

research/native-word-fft/EXP-IMPLEMENTATION-REVIEW.md

@@ -0,0 +1,153 @@
+# Independent review of the signed native exponential table
+
+Reviewed exp_table.py read-only. This is a mathematical and finite-word review,
+not Lean verification or independent human peer review. Apache-2.0;
+Dugongue with substantial OpenAI Codex assistance.
+
+## Verdict and one scope correction
+
+The approximation scheme is sufficient for the stated strict 2^-p error under
+its declared dyadic-input/range contracts. The construction uses O(log n)-bit
+words, but its table contains J+1=Theta(n^(1/4)log n) words along the growing
+parameter sequence. It does NOT have O(log n) total memory words. Arithmetic
+scratch outside the stored table is O(1) words. Each query executes at most
+833 Taylor iterations, hence constant work for the fixed bounded c,C family.
+
+## Parameters
+
+Let ell=bitlength(n+2)>=2, 1<=c,C<=8, k=c+4C<=40,
+
+    p=c*ell, hbits=floor(ell/4), h=2^-hbits,
+    J=C*ell*2^hbits,
+    d=8(k+64)+1, an odd integer at most 833,
+    jbits=bitlength(J+1), dbits=bitlength(d),
+    P=k*ell+jbits+dbits+32, S=2^P.
+
+An input argument is exactly the signed dyadic integer argument/2^fraction_bits,
+with 0<=fraction_bits<=P and absolute value at most C*ell. These are actual
+input-domain restrictions, not consequences of arbitrary list entries fitting.
+
+## Taylor truncation holds for every n, without runtime factorials
+
+The alternating exp(-x) remainder for 0<=x<=h after odd degree d is between
+zero and h^(d+1)/(d+1)!. We prove this is at most 1/S.
+
+For ell>=8, hbits>=ell/8 and therefore
+
+    hbits*(d+1) >= (k+64)*ell.
+
+A convenient bound is jbits<=2ell+4: C*ell<=8ell<=2^(ell+2) for ell>=2,
+and multiplying by 2^hbits followed by the harmless +1 gives that inequality.
+Also dbits<=10. Hence
+
+    P <= (k+2)*ell+46 < (k+64)*ell.
+
+The h power alone already gives a tail at most 2^-P; ignoring the factorial
+only weakens the estimate.
+
+For 2<=ell<=7, J<=112, so jbits<=7 and P<=7k+49. Put m=d+1, an even integer.
+The last m/2 factors of m! are at least 2, so
+
+    m! >= 2^(m/2) = 2^(4k+257) >= 2^(7k+49) >= 2^P,
+
+because k<=40. Again the omitted h factor only strengthens this bound.
+Thus the check() factorial inequality has an all-input proof; no factorial
+construction or unbounded factorial integer is needed by the native algorithm.
+
+## Floor recurrence and downward seed
+
+For x/S=theta in [0,h], let u_i=theta^i/i!, t_0=S and
+
+    t_i=floor(t_(i-1)*x/(S*i)).
+
+Set epsilon_i=u_i-t_i/S>=0. Then
+
+    epsilon_i <= (theta/i)epsilon_(i-1)+1/S.
+
+In the alternating sum, pair consecutive errors instead of adding their
+absolute values independently. For every odd i>=3,
+
+    epsilon_i-epsilon_(i-1) <= 1/S,
+
+and for every even i,
+
+    epsilon_(i-1)-epsilon_i >= -1/S.
+
+Also epsilon_1<=1/S. Therefore the floored degree-d polynomial differs from
+the exact degree-d alternating polynomial by between -(d-1)/(2S) and
+(d+1)/(2S). Adding the <=1/S Taylor remainder yields a conservative absolute
+polynomial error <=(d+3)/(2S).
+
+Subtracting d+2 integer grid units as done in q makes q/S <= exp(-h), while
+
+    0 <= exp(-h)-q/S <= (2d+4)/S.                         (1)
+
+The clamp at zero is inactive: exp(-h)>=exp(-1)>1/3, whereas S exceeds
+2^32*d by the precision definition, so the right side of (1) is far below
+1/12. In particular q>S/4>0. Terms are nonincreasing for theta<=1, so every
+integer alternating partial sum used by poly lies between 0 and S; its actual
+SUB operations need no negative intermediate.
+
+## Row propagation and residual query
+
+Let R_0=S and R_j=floor(R_(j-1)*q/S). The downward seed and downward products
+give
+
+    0 <= exp(-jh)-R_j/S <= j(2d+5)/S.                    (2)
+
+This follows inductively by one multiplication discrepancy plus its <1/S
+rounding. For a query magnitude a, write a=jh+t, 0<=t<h; at the exact upper
+endpoint a=C*ell, j=J and t=0. Thus the row index is always in [0,J].
+
+The residual poly has error <=(d+3)/(2S), and both row/poly values are in [0,S].
+Multiplying and flooring once more shows that its returned v obeys
+
+    |v/S-exp(-a)| <= (J+1)(2d+8)/S
+                    < 2^(jbits+dbits+2)/S
+                    <= E=2^(jbits+dbits+4)/S.            (3)
+
+The strict middle inequality uses d>=553, so 2d+8<4d. The +4 envelope is
+slightly stronger than the proposed +5 envelope and leaves an additional
+bit of slack for the positive-argument reciprocal below.
+
+## Positive-argument reciprocal and final rounding
+
+Let u=exp(-a), where 0<=a<=C*ell. Since exp(1)<4,
+
+    u > 2^(-2C*ell).
+
+Substituting P in (3) gives E=2^(-p-4C*ell-28), so E<u/2. Consequently v>0,
+v/S>=u/2, and
+
+    |1/(v/S)-1/u| <= 2E/u^2 < 2^(4C*ell+1) E.
+
+The literal floor(S^2/v) adds less than 1/S to the normalized reciprocal.
+Thus the total pre-final-rounding error is strictly less than
+
+    2^(-p-27)+1/S < 2^(-p-26).
+
+For nonpositive arguments no reciprocal is taken and (3) is smaller still.
+Nearest rounding from P fractional bits to p adds at most 2^(-p-1), proving
+strict final error below 2^-p. This covers zero and both signs. It does not
+justify a relative error contract; the established bound is absolute.
+
+## Bounded words, setup and query complexity
+
+All poly/row multiplicands are between 0 and S, so their products are <=S^2.
+A denominator S*i is at most S*d. The reciprocal numerator is exactly S^2,
+which needs 2P+1 bits. The positive reciprocal value is at most
+S*2^(2C*ell+2), well within the chosen 2P+32-bit word because P>=4C*ell.
+The argument shift is at most C*ell*S, the exact endpoint row J exists, and
+P>hbits ensures all shifts are nonnegative. Final-rounding addition has the
+same sufficient width margin. A common global machine word width must cover
+this subroutine's 2P+32 and all other caller subroutines/addresses.
+
+There are J+1 stored rows and constant-size row steps. Computing the initial
+poly and a query residual poly each takes d<=833 iterations. Since
+2^hbits<=2^(1/4)(n+2)^(1/4), setup is O((n+2)^(1/4)log(n+2)), each query is
+O(1), and total stored memory is O((n+2)^(1/4)log(n+2)) words. Word width is
+O(log(n+2)). No Decimal or factorial is called by construction/query; both
+occur only in the independent finite check().
+
+No mathematical implementation correction was required. Documentation should
+retain the stored-table memory bound and the precise dyadic/range contract.

research/native-word-fft/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.