From 4f46c6a9732f63cba411698c743f6b6b681273c9 Mon Sep 17 00:00:00 2001 From: Dugongue <49054532+Dugongue@users.noreply.github.com> Date: Wed, 21 Oct 2026 23:09:13 -0700 Subject: [PATCH] Develop native Word RAM Gaussian adapters with chirp blocks and reusable band inverses --- research/native-word-fft/.gitattributes | 1 + research/native-word-fft/CONTENT-HASHES.json | 38 +++ .../EXP-IMPLEMENTATION-REVIEW.md | 153 ++++++++++++ research/native-word-fft/FFT-ERROR.md | 159 +++++++++++++ research/native-word-fft/FINAL-REVIEW.md | 37 +++ research/native-word-fft/GAUSSIAN-CHIRP.md | 219 +++++++++++++++++ .../native-word-fft/INDEPENDENT-REVIEW.md | 138 +++++++++++ .../native-word-fft/KERNEL-PRECEDENTS.json | 50 ++++ research/native-word-fft/LICENSE | 201 ++++++++++++++++ .../native-word-fft/NATIVE-BANDED-INVERSE.md | 116 +++++++++ research/native-word-fft/NATIVE-EXP-TABLE.md | 31 +++ .../native-word-fft/NATIVE-OUTER-ASSEMBLY.md | 59 +++++ research/native-word-fft/NATIVE-RECOUNT.md | 87 +++++++ research/native-word-fft/NATIVE-STABILITY.md | 118 +++++++++ research/native-word-fft/NOTICE.md | 37 +++ research/native-word-fft/OVERLAP.md | 43 ++++ research/native-word-fft/PHASE-TABLES.md | 221 +++++++++++++++++ research/native-word-fft/PROOF.md | 113 +++++++++ research/native-word-fft/README.md | 101 ++++++++ .../native-word-fft/RECOVERY-PRECEDENTS.json | 22 ++ research/native-word-fft/ROOTS-RECOVERY.md | 223 ++++++++++++++++++ research/native-word-fft/SOURCES.json | 18 ++ research/native-word-fft/TASK.json | 70 ++++++ research/native-word-fft/VALIDATION.md | 50 ++++ research/native-word-fft/VERIFY.json | 20 ++ research/native-word-fft/band_solve.py | 110 +++++++++ research/native-word-fft/chirp_check.py | 51 ++++ research/native-word-fft/exp_table.py | 80 +++++++ research/native-word-fft/gaussian_block.py | 62 +++++ research/native-word-fft/gaussian_inverse.py | 92 ++++++++ .../native-word-fft/native_convolution.py | 73 ++++++ research/native-word-fft/native_fft.py | 165 +++++++++++++ research/native-word-fft/periodic_solve.py | 92 ++++++++ research/native-word-fft/phase_tables.py | 121 ++++++++++ research/native-word-fft/recover.py | 130 ++++++++++ research/native-word-fft/root_seed.py | 107 +++++++++ research/native-word-fft/verify.py | 54 +++++ 37 files changed, 3462 insertions(+) create mode 100644 research/native-word-fft/.gitattributes create mode 100644 research/native-word-fft/CONTENT-HASHES.json create mode 100644 research/native-word-fft/EXP-IMPLEMENTATION-REVIEW.md create mode 100644 research/native-word-fft/FFT-ERROR.md create mode 100644 research/native-word-fft/FINAL-REVIEW.md create mode 100644 research/native-word-fft/GAUSSIAN-CHIRP.md create mode 100644 research/native-word-fft/INDEPENDENT-REVIEW.md create mode 100644 research/native-word-fft/KERNEL-PRECEDENTS.json create mode 100644 research/native-word-fft/LICENSE create mode 100644 research/native-word-fft/NATIVE-BANDED-INVERSE.md create mode 100644 research/native-word-fft/NATIVE-EXP-TABLE.md create mode 100644 research/native-word-fft/NATIVE-OUTER-ASSEMBLY.md create mode 100644 research/native-word-fft/NATIVE-RECOUNT.md create mode 100644 research/native-word-fft/NATIVE-STABILITY.md create mode 100644 research/native-word-fft/NOTICE.md create mode 100644 research/native-word-fft/OVERLAP.md create mode 100644 research/native-word-fft/PHASE-TABLES.md create mode 100644 research/native-word-fft/PROOF.md create mode 100644 research/native-word-fft/README.md create mode 100644 research/native-word-fft/RECOVERY-PRECEDENTS.json create mode 100644 research/native-word-fft/ROOTS-RECOVERY.md create mode 100644 research/native-word-fft/SOURCES.json create mode 100644 research/native-word-fft/TASK.json create mode 100644 research/native-word-fft/VALIDATION.md create mode 100644 research/native-word-fft/VERIFY.json create mode 100644 research/native-word-fft/band_solve.py create mode 100644 research/native-word-fft/chirp_check.py create mode 100644 research/native-word-fft/exp_table.py create mode 100644 research/native-word-fft/gaussian_block.py create mode 100644 research/native-word-fft/gaussian_inverse.py create mode 100644 research/native-word-fft/native_convolution.py create mode 100644 research/native-word-fft/native_fft.py create mode 100644 research/native-word-fft/periodic_solve.py create mode 100644 research/native-word-fft/phase_tables.py create mode 100644 research/native-word-fft/recover.py create mode 100644 research/native-word-fft/root_seed.py create mode 100644 research/native-word-fft/verify.py diff --git a/research/native-word-fft/.gitattributes b/research/native-word-fft/.gitattributes new file mode 100644 index 0000000000..625449502b --- /dev/null +++ b/research/native-word-fft/.gitattributes @@ -0,0 +1 @@ +* -text diff --git a/research/native-word-fft/CONTENT-HASHES.json b/research/native-word-fft/CONTENT-HASHES.json new file mode 100644 index 0000000000..1a72cacbf5 --- /dev/null +++ b/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" +} diff --git a/research/native-word-fft/EXP-IMPLEMENTATION-REVIEW.md b/research/native-word-fft/EXP-IMPLEMENTATION-REVIEW.md new file mode 100644 index 0000000000..053a348d50 --- /dev/null +++ b/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=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 E0, +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. diff --git a/research/native-word-fft/FFT-ERROR.md b/research/native-word-fft/FFT-ERROR.md new file mode 100644 index 0000000000..53691b10de --- /dev/null +++ b/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. diff --git a/research/native-word-fft/FINAL-REVIEW.md b/research/native-word-fft/FINAL-REVIEW.md new file mode 100644 index 0000000000..ae5eb6e825 --- /dev/null +++ b/research/native-word-fft/FINAL-REVIEW.md @@ -0,0 +1,37 @@ +# Independent final arithmetic review + +Reviewed the package at `root/repo/research/native-word-fft` without editing its sources. The reviewer owns only this arithmetic review. No correctness defect was found in the inspected Gaussian block, actual Gaussian coefficient, periodic Woodbury, or rounded diagonal-dominance arguments, under their stated precision, near-identity and common-word-width hypotheses. + +## Concrete checks and mathematical scope + +A fresh optimized-mode replay of `gaussian_block.py` passed 12 blocks and 45 outputs at absolute tolerance 2^-60. The executed family is a=1/16, u=v=1, q in {1,2,4,8}, and integer shifts c in {-1,0,1}; the universal block decomposition is a separate mathematical result. The reversed source and convolution index x+q-1 implement the diagonal-Hankel-diagonal identity correctly. With transform length 4q, source length q and kernel length 2q-1, there is no cyclic aliasing. General Convolver callers must themselves supply a prepared transform of the same length and sufficient zero padding; the Gaussian caller does so. + +A fresh replay of `gaussian_inverse.py` passed the 3 systems (s,t)=(13,14),(17,19),(25,27), alpha=2, h=4, over 55 output rows at tolerance 2^-100. The maximum finite-reference residual was approximately 4.604745e-51. Its reference sums only periodic copies -3 through 3 using Decimal; this is finite validation, not an independent proof about the infinite sum. + +The actual coefficient numerator is correct: shifted=t(i+d)-s*q_i represents s*(r*d+beta_i), while beta[(i+d) mod s] supplies the periodic column label. Keeping i+d unwrapped in shifted and wrapping only the label is necessary and is implemented. Machin pi uses table.S=2^P; querying the resulting argument with fraction_bits=P has the correct scale. Negative exponential cutoff is explicit. Its omitted error must be charged separately, at most (2h+1)*exp(-C*ell) per row in this finite-band implementation; the fixed test parameters have ample margin. Arbitrary caller precision choices must enforce this budget, as well as the infinite Gaussian tail and pi/exponent-formation errors. + +The periodic reduction uses literal A=A0+UV with V selecting the first and last h columns. Stored X=A0^-1 U and S=I+VX have the correct scales; the correction Xz is divided by 2^p once. Boundary columns are distinct when s>2h. All 2h ordinary band solves, the dense small factorization and final correction are present. The executable contract is real coefficients with ||A-I||_infinity<=1/8, narrower than the universal proof's complex 1/4 interface. In particular, its exact small Schur matrix differs from I by at most 1/7, leaving room for the real BandSystem's 1/4 acceptance check when the advertised guard precision is used. + +## Universal rounded band proof + +The exact row diagonal-dominance margin does not decrease under Schur elimination; its active row sum does not increase. Rounded updates perturb each active row in at most h elimination steps with at most 2(h+1) entries per step. The stated 6(h+1)^2*u allowance closes the row-sum/pivot induction. Expanding the actually stored factor updates gives a local residual bound without multiplying successive errors exponentially. The removed pivot-column residual includes multiplier rounding times its pivot. The h=0 case is directly diagonal; for h>=1 the displayed residual constant dominates (H+2)(2H+1). + +The factor norm identities then bound both triangular inverses. Leading-prefix forward bounds follow from principal factor pairs. Backward suffix inverse bounds also follow by taking the corresponding principal submatrix of the full upper-triangular inverse. Exact fixed-point dot accumulation followed by one coordinate rounding gives the displayed solve residual equations and conservative final error constant. A fixed number of logarithmic-width limbs covers products, complex division and dot accumulation when the specified integral and fractional guards are allocated. + +For the actual periodic Gaussian, the tail inequality follows from |r*d+beta_i|>=|d|-1/2 and |beta_j|<=1/2. The wraparound Woodbury reduction and its constant inverse-norm bounds are valid under the credited source estimate ||N-I||<=1/8 and sufficiently small coefficient/rounding error. Dense small-s handling and the complex extension are proof-level branches, not executed features of these real examples. + +## Model and attribution limits + +Production code obtains roots through the finite Machin/Taylor implementation and exponentials through the paid finite table. Decimal and Fraction are check-only references. Inspected arithmetic products are bounded-coordinate, bounded-parameter or bounded-index products checked against an explicit word limit; there is no unit-cost multiplication of the entire input integer. Python remains an executable bounded-word model rather than an emitted register-machine program, and its operation counters are accounting instruments rather than a formally verified compiler cost certificate. Global caller blocking, requested error budgets and a common word width remain explicit assembly hypotheses. + +This review establishes no new global multiplication bound, no numerical kappa improvement and no novelty priority. Classical FFT, diagonal dominance and Woodbury arguments, upstream Gaussian resampling and its cited Harvey-Hoeven work retain their predecessor credit. Source identities and localization/contraction results must remain pinned and credited in the integrated proof. + +## Reviewed SHA-256 values + +- native_convolution.py: A62A4CB5A1DA458E31C90EA881A1146B59E9FD7CCBF317C1C9836C19BC7FA2C6 +- gaussian_block.py: EB967DDE22E1989EF68057F7644BFF9519C0426137054E2517766FDDEE9D24DF +- gaussian_inverse.py: 08616D26652C7C656ED8928A2C47A0A150A8CFB35BCC7F121CE5A729908E4A35 +- periodic_solve.py: D3940B0D48BA4BBF81A0F2D05BDAC29E3B998AD746771DCF01B03CAC03680393 +- band_solve.py: F4E4471679CC89BBB01AB0B26E807B0DACDC6188AB95BA924FC57243369407BB +- NATIVE-BANDED-INVERSE.md: 1259B4242EBF6AC7F14BFE72E371AFB8DC530B2755D1FD82228D9DC6F01634BA +- exp_table.py: E617B40D499F423A7861B4483BDACDB7DB28232FADF7468EFE871ADAE557FA35 diff --git a/research/native-word-fft/GAUSSIAN-CHIRP.md b/research/native-word-fft/GAUSSIAN-CHIRP.md new file mode 100644 index 0000000000..4e31f30215 --- /dev/null +++ b/research/native-word-fft/GAUSSIAN-CHIRP.md @@ -0,0 +1,219 @@ +# Stable Gaussian block maps by one paid convolution + +Prepared for Dugongue with substantial OpenAI Codex assistance. Apache-2.0. +This is a constructive native Word RAM theorem. It uses no Turing machine, +transition table, machine simulation or free convolution/root/exp oracle. + +## 1. Exact diagonal-Hankel-diagonal identity + +Let a>0 and u,v,c be real. For integer row and column coordinates x,y, define + + W[x,y] = exp(-a*(u*x-v*y+c)^2). + +On a rectangle 0<=x=m+2*l-2. For a power-of-two q, +N0=4q always suffices. Perform the ordinary radix-two forward FFT of reversed +E*Z, multiply by the prepared transform of H, and perform the inverse FFT +with its 1/N0 normalization. Select indices x+l-1 and multiply by D[x]. +The fixed kernel's transform is also paid; it is not a free arbitrary-kernel +convolution. Without sharing H, this is three FFTs per block; with sharing, +it is two per block plus one shared FFT. + +The FFT is the direct iterative bounded-word algorithm from native_fft.py and +FFT-ERROR.md: signed dyadic arithmetic, a generated root seed, explicit +bit-reversal/butterflies, word-sized products, and exact-grid addition. The +root seed and required root powers are generated once for N0 at guard precision. +There is no TM multiplier call. The pure arithmetic algorithm and its error +proof apply to arbitrary complex bounded inputs; integer recovery is not used +inside this Gaussian block. + +For clarity, a conservative sufficient error budget is + + block error <= 2^30*q^8*M^4*(X+1)^2 * 2^-P, (B) + +when all factor/input, root and local rounding errors are at most 2^-P +at their stated certified interfaces. One may use a larger fixed constant. +Choose + + P >= p0 + 30 + 8*ceil(log2(q+2)) + + 8*ceil(C)*p + 2*ceil(log2(X+1)) + 4. (C) + +The integer coefficient 8*ceil(C) safely bounds 4*C*log2(exp(1)); no transcendental +constant is needed to choose precision. Choose this common precision once from +the global parameter and input bounds, rather than recomputing logarithms at +every block. Rounding P up and adding a fixed guard for the supplied exp +implementation makes (B) strictly below 2^-p0. Exponent arguments are computed +with enough additional guard bits that their exp(Cp)-Lipschitz images meet the +factor-error interface; that costs only O(p+log(t+s)) extra bits. The common +physical word width covers this larger internal precision as well. Thus all +working precisions remain O(log n). + +**Proof of (B).** Approximate D,E,H within epsilon=2^-P and use rounded +complex products with error <=2*epsilon. Factor magnitudes are at most M, +and their approximants at most 2M. The preweighted input has magnitude at +most 2M*(X+1) after a fixed rounding guard. Replacing all three exact factors +by their approximants changes each D*E*H product by at most 7M^2*epsilon; +rounding pre/post products and input contributes further constant multiples of +q*M^3*(X+1)*epsilon. The independent three-transform error proof supplies +at most 216*N0^6*beta0^2*epsilon for the convolution, where +beta0<=4M*(X+1), using correspondingly certified twiddles. Multiplication by +D contributes at most 2M times that error, plus its own factor/rounding terms. +Using N0=4q and absorbing the fixed constants and lower q powers yields (B). +The same estimates bound all intermediates by a fixed polynomial in q,M,X. +Since their logarithms are O(log n), both raw word products and their sums fit +a common O(log n) width, including the root-seed implementation's larger guard +width. This proves the constructive block bound. ∎ + +If this map's global output sums at most K such rectangles per row, replace +p0 by p0+ceil(log2(K+1))+2. For fixed K this is a constant extra guard. +If the input and arguments themselves are dyadic approximations, their errors +must be included in the certified factor/input interfaces in (B); they are not +assumed exact simply because they are held in words. A polynomial Lipschitz +bound exp(Cp) times the exponent-argument error supplies that guard. + +## 4. Rectangular blocking and the complete fast banded Gaussian map + +Let 0<=x0 have bounded ratio in both +directions. Fix a radius R and retain a whole q-by-q block whenever some pair +in it satisfies |u*x-v*y+c|<=R. Sum its full Gaussian rectangle; do not apply +a per-entry band mask that would destroy the Hankel identity. + +For each row block, the possible column coordinates form an interval of length +at most [u*(q-1)+2R]/v. It intersects at most + + K = ceil((u*(q-1)+2R)/(v*q)) + 3 + +column blocks. If q is within fixed constant factors of R/max(u,v), K is a +fixed constant. There are O(t/q+1) row blocks, hence O(t/q+1) retained +rectangles. For an active rectangle, + + |c0| <= R+(u+v)*(q-1), + +so its conditioning follows from + + a*(R+2*(u+v)*q)^2 <= C*p. (D) + +Choosing a*R^2=Theta(p) and q a power of two within a fixed factor of +R/max(u,v) ensures (D). Require q>=1; when the resulting q exceeds the array +length, use one appropriately smaller block and the direct small case. + +**Theorem 3.** Under these explicit geometric hypotheses, the retained map +costs + + O(t log(q+2) + s + n^eta polylog n) + +native Word RAM operations, including shared root/exp setup, generated chirps, +all block FFTs, gathers/scatters, summation and IO. When q is a positive power +of log n, the dominant map work is O(t loglog n+s), rather than the naive +O(t*q) neighbor enumeration. Every term is paid. + +**Proof.** Multiply O(t/q+1) rectangles by the O(q log(q+2)) block price, add +O(q log(q+2)) shared fixed-kernel FFT setup, O(t+s) data IO and gathering, +and the n^eta polylog n exp preparation. Use q<=t in the nonsmall branch. +Geometry bounds K uniformly, so summing blocks adds O(t) word work. ∎ + +For the full untruncated Gaussian map, every omitted weight is at most +exp(-a*R^2). Its row error is at most + + s*X*exp(-a*R^2). + +Thus a*R^2 >= (p0+ceil(log2(s*X+1))+2)*log(2) makes omitted tails less than +2^-p0/2. This is compatible with (D) when s,X are polynomially bounded in n +and p0=Theta(log n). Add one precision bit to the block summation budget and +obtain a constructive approximation to the full map. No Gaussian tail is +ignored merely because its entries look numerically small. + +## 5. Scope and actionable use + +This is a positive sufficient replacement for a concrete forward Gaussian +resampling map once its actual grid slopes, smoothing parameter, input bound, +row/column ranges and requested error instantiate (A) and (D). It removes the +naive q factor from the local convolution part with a wholly native arithmetic +algorithm. It does not by itself certify every Gaussian row in the multiplication +source or establish a new global kappa. The root must bind the block layout and +precision to the chosen complete native caller before claiming such a bound. + +The exp tables and root tables are numerical function tables generated from +bounded arithmetic. They are not TM transition tables or simulations. The +Gaussian chirp identity, ordinary FFT convolution, and radix blocking have +classical precedents; no worldwide novelty is asserted. Our result is the +explicit stable finite-word interface and complete cost/error assembly. + +A power-of-two block size satisfying the rationally expressed normalized-span +bound can be found by a doubling loop and squared comparisons. The program +need not receive an exact real square-root primitive. A chosen radius may also +be specified by a rational squared bound and a proved tail inequality. The +actual source instantiation must supply these bounds rather than relying on +an uncomputed real-valued parameter selection. diff --git a/research/native-word-fft/INDEPENDENT-REVIEW.md b/research/native-word-fft/INDEPENDENT-REVIEW.md new file mode 100644 index 0000000000..3a620766f7 --- /dev/null +++ b/research/native-word-fft/INDEPENDENT-REVIEW.md @@ -0,0 +1,138 @@ +# Independent arithmetic review + +Read-only review by the arithmetic worker for the current native Word RAM +campaign. Prepared with OpenAI Codex assistance. This is an informal mathematical +review of current source bytes, not a Lean or human peer-review certificate. + +## Root seed: accepted within its stated model + +Inspected `recovery/root_seed.py` and `recovery/ROOTS-RECOVERY.md` directly. +The construction uses bounded signed magnitudes, fixed-point integer arithmetic, +Machin's identity and a finite Taylor loop. Decimal appears only in the separate +validation routine; it is not called by `root_seed`. + +The following all-parameter points were checked independently: + +1. With P=8(p+16), stopping the alternating arctangent sum when the next + q^(2j+1) exceeds S=2^P gives omitted tail <1/S. There are at most P+1 + retained terms. Floors add at most one scaled unit per retained term. + The Machin combination therefore has error at most 20(P+2)/S. +2. The integer angle floor(2*A/D) differs from 2*pi/D by at most + [40(P+2)+1]/S. It lies in [0,8] after scaling, including D=1. +3. The term recurrence uses the old real and imaginary coordinates in both + numerators before assigning the new pair; no sequential-update error occurs. + Its coordinate error recurrence is E_k <= (8/k)E_(k-1)+1/S. Unfolding it + and using j!/k! <= 1/(k-j)! gives E_k < exp(8)/S <4096/S. +4. For K=128(p+1), the omitted Taylor tail is at most + 2*8^(K+1)/(K+1)!, bounded by 2^(-K/2-3). The last K/2+1 factorial + factors are at least 128. The inequality 2^(-K/2-3)<=32/2^P holds for + every p>=1; p=1 is the limiting equality in this last comparison. +5. The printed internal error numerator simplifies exactly to + 524608*p+529521. For p>=1 this is <2^21*p. Since p<=2^p, + its product with 2^(p+1) is <2^(2p+22)<2^(8p+128). Thus the final + nearest-p-bit rounding proves strict coordinate error <2^-p for every p, + not just for the finite p fixtures. +6. Machin denominators have at most P+O(log P) bits and its next power is + at most 239^2*S. Taylor values are bounded by a fixed multiple of S; + the stated 2P+20 physical word width covers the 2P+17 product envelope, + denominators, accumulators and rounding. D must also fit, as stated. +7. Both arctangent loops have O(P) iterations and the Taylor loop has K=O(p) + iterations, each using a constant number of bounded-word MUL/DIV and + ordinary operations. Therefore O(p) setup is justified in the declared + unit-cost logarithmic-word model. It is not a claim about bit operations. + +The same physical word width must be used by the whole caller. A private wider +seed word model cannot be hidden inside a narrower FFT model. With p=Theta(log n) +and D polynomial in n, the explicit seed width is still Theta(log n), though +its constant is intentionally large. + +## Native FFT: accepted within its stated operational-model scope + +Read the actual `root/native_fft.py`, its current root seed, and the integrated +`FFT-ERROR.md`. The executable now includes the separate 65,535-bit extremal +fixture. This review does not relabel its operational model as an emitted +RamStd program or formal compiler certificate. + +### Correctness and numerical sufficiency + +The bit-reversal permutation and iterative radix-two loops use ordinary +unnormalized forward transforms with negative exponent and inverse with positive +exponent. Twiddles are prepared from the positive primitive root and conjugated +for forward calls. Root indices k*(N/step) lie below N/2, so the prepared table +contains every requested value. Three transforms and one pointwise multiplication +compute the zero-padded ordinary convolution when N>=2L. + +I independently checked the integrated FFT-ERROR.md recurrences. With its +coordinate-to-complex error conversion, the one-transform error bound +36*N^3*beta/S, pointwise error 112*N^5*beta^2/S, and normalized final error +216*N^6*beta^2/S are conservative valid envelopes. The use of approximate +rather than exact second-input amplitude in the twiddle perturbation term is +correct. Its closed larger envelope and the substituted F=24*ell+64 prove a +strict error below 1/4 for every n>=1. The seed/table error supplies the stated +twiddle hypothesis; rounded complex arithmetic supplies the local error. + +For an additional independent route, use complex error tau<=16*N/S for the +successively generated twiddles and local error <=2/S. A transform with ideal +input magnitude X and input error E has output error at most + + N^2*E + N^3*tau*X + 2*N^2/S. + +This follows from per-stage error <=3E+tau*N*X+2/S and 3^log2(N)<=N^2. +It yields forward error <=18*N^4*beta/S, pointwise error +<=56*N^5*beta^2/S, and normalized final error <=76*N^6*beta^2/S. +The forward-error assumption <=N*beta used here is implied by the chosen F. +Since N<4n<2^(ell+2) and beta<=2^ell, even this separate final bound is +less than 2^(8ell+19-F)<1/4. + +Consequently the runtime separation checks are backed by a universal error +proof; checking proximity to some integer alone would not have been enough. +They do not abort on valid inputs. The nearest integer is the exact convolution +coefficient, the high zero convolution coefficients really vanish, and the +subsequent carry loop and output high-padding checks are justified. + +### Word bounds and operation model + +Every MUL in `multiply` is on bounded word data: four fixed-point coordinate +products in complex multiplication, limb/scale products, or polynomially bounded +index/range arithmetic. The independent unbounded product x*y appears only in +`check`, never in the algorithm. Padding by scale can equivalently be a left +shift. Powers of two used for normalization and carry can be shifts/masks; +ordinary bounded-word DIV is permitted by the declared model in any case. + +A crude actual inverse-prefix amplitude bound is 6*N^6*beta^2. Transform +products and combined real/imaginary sums thus have magnitude width at most +2F+6*ceil(log2 N)+2b+8. Pointwise forward products fit that same bound. +This is covered by the single global width 16*(F+16)+20, which also covers +the root seed's larger working precision. Integer counters, indices, polynomial +coefficient-range bounds, array lengths and fixed-table pointers fit too. + +`Words.limit`, overflow diagnostics, Python containers and the independent +reference are model/checking machinery, not an assertion that the emitted +machine stores an unbounded integer or obtains transcendental constants. +The universal bounds justify removing overflow checks in a future machine +lowering. Literal `charged_word_operations` are displayed operational-model +charges, not an exact count from an emitted instruction stream. + +### Direct cost and scope + +With b=max(1,floor(ell/8)), L=ceil(n/b) and N the least power of two at least +2L, we have L=Theta(n/log(n+2)), N=Theta(L), and log N=O(log n). +Bit reversal uses O(N) aggregate mask flips. Three transforms cost O(N log N) +word work; root preparation O(F), root-table recurrence O(N), padding, +pointwise products, recovery, carry and packed IO cost O(N+L). Thus the whole +direct native algorithm costs O(n), with every preparation and IO stage paid. +It uses neither TM simulation nor a TM multiplication antecedent. + +This is a complete positive native arithmetic baseline, not a new best-known +asymptotic exponent: the cited log-RAM literature already contains O(n) +algorithms. The separate stable-native-provider and phase-table lemmas do not +by themselves establish a stronger all-length cyclic transform or new kappa. +An actual stronger provider still needs its own quantitative full-conditioning, +preparation, depth, and caller-boundary admission. + +## Reviewed byte identities + +- `root_seed.py` (recovery): `9b323f84ade7decfa6152dd4385e90ec43b56087927ac23ad4ef73ae4d5d3042` +- `native_fft.py` (root): `83bff7591e4ffa36fe3713f9d639de6ed1edc715d3698bc6d2d0dc3b90146cc5` +- `root_seed.py` (root): `9b323f84ade7decfa6152dd4385e90ec43b56087927ac23ad4ef73ae4d5d3042` +- `FFT-ERROR.md` (native-word-fft): `ca426930df5c449ca8ce432360ba3ac80f12377bcd9a058ad16387d8ad87bef5` diff --git a/research/native-word-fft/KERNEL-PRECEDENTS.json b/research/native-word-fft/KERNEL-PRECEDENTS.json new file mode 100644 index 0000000000..fdf91b89c0 --- /dev/null +++ b/research/native-word-fft/KERNEL-PRECEDENTS.json @@ -0,0 +1,50 @@ +{ + "sources": [ + { + "pr": 956, + "source_kind": "previously fetched proof text", + "commit": "fa4df77c9dfdbac93368d95dcb361e69b3b65464", + "sha256": "33949793c1044204fec332174071714bb085f32054062fcfb3d56fb84dd0a8d9", + "scope": "informal source inspection only", + "path": "research/pr390-heterogeneous-provider-closure/PROOF.md", + "url": "https://github.com/CrocSwap/integer-mult-bounds/pull/956" + }, + { + "pr": 647, + "source_kind": "previously fetched proof text", + "commit": "3b3fda9c602fbe4b45d02d1ead0b757b3c94f52c", + "sha256": "cdc496acb844de3b3bcf03f811d1bd345e476b40f0b1cb7c866066612770e834", + "scope": "historical proof source inspection; full public head rechecked by API this session, no new replay" + } + ], + "new_proof_sha256": "260e78bf8f88667a637b804d013d1e7c6dfdf1e23ae35ffdc64a4366c7321221", + "claim": "native fixed-rational kernel lowering; not cyclic multiplication or a new kappa", + "proof_files": [ + { + "filename": "NATIVE-BANDED-INVERSE.md", + "sha256": "1259b4242ebf6ac7f14bfe72e371afb8dc530b2755d1fd82228d9dc6f01634ba", + "bytes": 10455 + }, + { + "filename": "NATIVE-EXP-TABLE.md", + "sha256": "c9be6814a81e0bfe6b041f1a259a32369f7c030ca22eaa5c71c15911c888c6ff", + "bytes": 5376 + }, + { + "filename": "NATIVE-OUTER-ASSEMBLY.md", + "sha256": "f1e194c01fc1ba4e33c698dbc287aeca8ff0efe822b195c50c95b619242cf9dc", + "bytes": 5613 + }, + { + "filename": "NATIVE-RECOUNT.md", + "sha256": "fc14d31113ef3abd597e654e0aec35ef7972d014268e4b8b7624f202f9ff79db", + "bytes": 10884 + }, + { + "filename": "NATIVE-STABILITY.md", + "sha256": "260e78bf8f88667a637b804d013d1e7c6dfdf1e23ae35ffdc64a4366c7321221", + "bytes": 11321 + } + ], + "scope": "Direct native Gaussian rational kernel, exponential tables, actual periodic Gaussian inverse, and explicit sufficient native Fourier outer assembly; not a new best kappa." +} diff --git a/research/native-word-fft/LICENSE b/research/native-word-fft/LICENSE new file mode 100644 index 0000000000..261eeb9e9f --- /dev/null +++ b/research/native-word-fft/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/research/native-word-fft/NATIVE-BANDED-INVERSE.md b/research/native-word-fft/NATIVE-BANDED-INVERSE.md new file mode 100644 index 0000000000..8be8f9922d --- /dev/null +++ b/research/native-word-fft/NATIVE-BANDED-INVERSE.md @@ -0,0 +1,116 @@ +# A native bounded-word inverse for a near-identity band matrix + +Let A be an s by s complex matrix with bandwidth h: a_ij=0 for |i-j|>h. Suppose A=I+E and ||E||_infinity<=rho<=1/2. All norms below use complex absolute values. The input entries are explicitly stored bounded-word fixed-point numbers; this is not an assumed inverse oracle. Put H=h+1. Gaussian elimination without pivoting, followed by forward and backward substitution, computes an approximation to A^-1 b in O(s H^2) bounded-word operations, O(s H) space and an explicit absolute error at most + + 1000 H^4 u (||b||_infinity+1), + +provided local complex rounding errors are at most u<=1/(1000 H^2), and enough integral bits are reserved as below. Constants are deliberately loose. To obtain output error <=2^-P it suffices to choose + + u <= 2^-P / [1000 H^4 (||b||_infinity+1)]. + +Thus P+log2(||b||+1)+4 log2 H+O(1) fractional bits suffice. If P, log s, log(||b||+1)=O(log n) and h<=n^C for fixed C, all operations use a fixed number of O(log n)-bit words. + +## Exact elimination has bounded growth + +Initially every row has diagonal-dominance margin + + |a_ii|-sum_(j!=i)|a_ij| >= 1-rho >=1/2, + +and absolute row sum at most 1+rho<=3/2. Consider elimination of column k in a remaining row i. Write l=a_ik/a_kk. For the Schur row after deleting column k, the triangle inequality gives + + new margin >= old margin + |a_ik|[1-offsum(row k)/|a_kk|] + >= old margin. + +Similarly its absolute row sum is at most + + old row sum - |a_ik| + |l| offsum(row k) + <= old row sum. + +Consequently all exact pivots are at least 1/2, every surviving row sum is at most 3/2, and every multiplier has magnitude at most 3. These assertions hold for every leading principal subproblem too. The inverse infinity norm is at most 2: if |x_i|=||x||_infinity, then |(Ax)_i|>= (|a_ii|-sumoff)|x_i| >= ||x||/2. + +No pivoting is needed. The lower and upper bandwidths remain h: a row only interacts with at most h preceding pivots and at most h following columns. This proves the claimed operation and space counts directly. + +## Rounded elimination retains these bounds + +Round a multiplier within absolute error u, and each Schur entry update within absolute error u after its exact bounded-word product and subtraction. As long as remaining row sums are at most 2, replacing an exact multiplier by the rounded multiplier perturbs an updated entry by at most 2u, and its final rounding adds at most u. There are at most 2H active entries per updated row and at most H elimination steps affecting that row. The total adverse change in either its row sum or its diagonal-dominance margin is therefore at most 6 H^2 u. The assumption u<=1/(1000 H^2) makes this less than 1/100. By induction the row sums stay below 2, the margins stay above 1/4, pivots have magnitude at least 1/4 and multipliers are bounded by 8. This closes the estimates used in the induction. + +Denote the computed factors by Lhat,Uhat. The factor residual is small without exponentially multiplying rounding errors. For one fixed row i, expand its successive updates literally: + + current a_ij = original a_ij - sum_(k already eliminated) lhat_ik uhat_kj + + sum(local update roundings). + +Each entry has at most H local update roundings. At the removed pivot column, rounding its multiplier contributes at most |uhat_kk|u<=2u. Thus each entry of Lhat Uhat-A has magnitude at most (H+2)u, with at most 2H+1 nonzero entries per row. In particular + + ||Lhat Uhat-A||_infinity <= 6 H^2 u. + +Let Ahat=Lhat Uhat. Its row dominance is at least 1/2-6H^2u>1/4, so ||Ahat^-1||_infinity<=4. Also ||Uhat||_infinity<=2 and ||Lhat||_infinity<=1+8h<=9H. Hence + + ||Lhat^-1||_infinity <= ||Uhat|| ||Ahat^-1|| <=8, + ||Uhat^-1||_infinity <= ||Ahat^-1|| ||Lhat|| <=36H. + +The same bounds hold for every computed leading principal pair. They control all forward-substitution prefixes, so no exponential intermediate-growth assumption is being hidden. + +## Rounded solves and total error + +For each row, accumulate its fixed-point dot product exactly with a bounded-word double-precision accumulator and round only the solved coordinate. Forward substitution then produces yhat with + + Lhat yhat = b+e_f, ||e_f||<=u. + +Thus ||yhat||<=8(||b||+u). Back substitution produces xhat with + + Uhat xhat = yhat+e_b, ||e_b||<=2u, + +because rounding a solved coordinate by at most u contributes at most |uhat_ii|u<=2u to that equation. Hence + + Ahat xhat = b+e_f+Lhat e_b, + ||e_f+Lhat e_b|| <=19H u, + ||xhat|| <=4(||b||+19H u). + +Compare with the true solution x=A^-1 b: + + ||xhat-x|| <=2[6H^2u*4(||b||+19H u)+19H u]. + +For u<=1/(1000H^2), this is bounded by the stated conservative 1000H^4u(||b||+1). All partial dot products have at most H terms, factors bounded by 8, y bounded as above and x bounded by the same argument on upper triangular suffixes or their solved equations. An additional O(log H+log(||b||+1)+1) integral bits therefore covers every exact sum and product; the exact rational complex division uses a pivot modulus at least 1/4. Fixed extra limbs hold products, accumulators and divisions. Every scalar update consequently costs O(1) bounded-word instructions. + +## Gaussian caller and credit + +If the actual Gaussian inverse map N has a rigorously truncated bandwidth h=O(sqrt(p)/alpha), truncation norm error eta and ||N-I||<=rho<1/2, apply this theorem to the stored truncated matrix A. Perturbing its right-hand solve changes the ideal inverse by at most ||N^-1||*eta*||A^-1|| times the right-hand norm, so tail and coefficient errors can be invoiced explicitly. Matrix formation, proof of the bandwidth and row-sum tail bound remain the caller's responsibility. This theorem does not assert those facts merely from the name 'Gaussian'. + +The elimination and diagonal-dominance arguments are classical numerical linear algebra. This document supplies a direct native sufficient algorithm and explicit word/precision/cost accounting; it does not claim a new best multiplication exponent or worldwide novelty. The Gaussian resampling caller is credited to the upstream construction and its cited Harvey-Hoeven sources. + + +## The actual periodic Gaussian inverse, with a constructive reduction + +Here the band hypothesis can be proved from the actual upstream formula, rather than postulated. Let s=2, r=t/s>1, q_i=floor(ri+1/2) and beta_i=ri-q_i in [-1/2,1/2]. The source defines N=C T D. Its entries are + + N_ij = sum_(ell congruent j mod s) + exp[-pi alpha^2 (r(ell-i)+beta_i)^2 + pi alpha^2 beta_j^2]. + +For d=|ell-i|>=1, each positive summand is at most + + exp[-pi alpha^2(d^2-d)], + +because |r(ell-i)+beta_i|>=d-1/2 and |beta_j|<=1/2. Truncate the sum to |ell-i|<=h. The omitted absolute row sum is at most + + 2 sum_(d>=h+1) exp[-pi alpha^2 d(d-1)] + <=4 exp[-pi alpha^2 h(h+1)]. + +The final inequality follows by bounding the successive-term ratio by exp[-2pi alpha^2(h+1)]<1/2. Choosing + + h = ceil(sqrt((P ln 2+ln 16)/(pi alpha^2)))+1 + +therefore makes the row-sum tail less than 2^(-P)/4. This is the required h=O(1+sqrt(P)/alpha), with no unknown localization assumption. For s>2h, each retained entry has at most one integer representative in the window, so formation takes O(s h) explicitly generated coefficients. In the actual assembly s is superpolynomial in the precision, hence this condition eventually holds. Smaller s may instead use dense O(s^3) elimination; when s<=2h this is O(s h^2). + +Suppose the actual source's proved estimate ||N-I||<=2^(-alpha^2 theta) is at most 1/8, as it eventually is in the intended regime. All entries are nonnegative real numbers and diagonal retained entries include 1. Removing the positive tail can only decrease the row norm of N-I. Form the retained periodic matrix A to coefficient row-sum error at most 1/32, by allocating O(log h) extra fractional bits per entry. Then ||A-I||<=1/4 and ||A^-1||<=4/3. Its difference from N is the explicitly controlled tail and coefficient-formation error; allocate both to the requested P-bit accuracy instead of fixing the error at 1/32 when final accuracy is required. + +To handle the wraparound without assuming it is an ordinary band matrix, let A0 be A with only its wraparound corner blocks removed. Those off-diagonal entries are nonnegative. With coefficient errors included conservatively, ||A0-I||<=1/4. The removed corner matrix can be written U V, where V selects the first and last h coordinates, and U has those 2h columns equal to the removed corner columns. Thus rank(UV)<=2h, ||V||<=1 and ||U||<=1/4. All these are literal stored matrices. + +Factor the ordinary band matrix A0 using the theorem above. Solve its system for b and for the 2h columns of U. Denote the results by y=A0^-1 b and X=A0^-1 U. The small matrix + + S=I+V X + +has ||S-I||<=1/3, because ||A0^-1||<=4/3. It is itself strictly row diagonally dominant, with inverse norm at most 3/2. Factor this explicit 2h by 2h matrix with the same constructive elimination, solve S z=V y, and return x=y-Xz. Direct multiplication verifies A x=b. Thus the algorithm is not hiding a small inverse oracle. + +The work is O(s h^2): ordinary band factorization costs O(s h^2), its 2h right-hand solves also cost O(s h^2), and the final dense factorization costs O(h^3)<=O(s h^2). Formation of Xz and every small-matrix operation is included. The storage is O(s h). All exact inverse/operator norms above are fixed constants; numerical approximations to X can use per-entry error budget epsilon/(2h), and all finite sums use O(log s+log h) extra guard bits. The explicit band error theorem then gives any requested final accuracy with P+O(log s+log h+log(||b||+1)) bits. This is O(log n) in the intended regime. Bounded-word coefficient preparation is provided by NATIVE-EXP-TABLE.md, including the charged setup. + +Finally, ||N^-1|| is bounded by 8/7 from the original 1/8 estimate. The identity N^-1-A^-1=N^-1(A-N)A^-1 controls truncation and formation errors directly by at most (8/7)(4/3)||A-N||. This completes a positive direct native algorithm for the actual Gaussian inverse, not merely a conditional band-matrix interface. The exact source localization, contraction and resampling identities remain predecessor facts; the native band/Woodbury evaluation and its accounting are the assembly given here. diff --git a/research/native-word-fft/NATIVE-EXP-TABLE.md b/research/native-word-fft/NATIVE-EXP-TABLE.md new file mode 100644 index 0000000000..ee6e734c82 --- /dev/null +++ b/research/native-word-fft/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/native-word-fft/NATIVE-OUTER-ASSEMBLY.md b/research/native-word-fft/NATIVE-OUTER-ASSEMBLY.md new file mode 100644 index 0000000000..b9c5496896 --- /dev/null +++ b/research/native-word-fft/NATIVE-OUTER-ASSEMBLY.md @@ -0,0 +1,59 @@ +# A direct native Fourier-to-integer assembly + +This is a sufficient assembly theorem with an explicit caller and an explicit missing transform hypothesis. Its multiplication construction uses no Turing-machine algorithm. The independently proved closed Gaussian kernel in NATIVE-STABILITY.md does not, by itself, meet the cyclic Fourier hypothesis below. + +## Input and arithmetic + +Fix a standard random-access machine with word width w=c_w ceil(log2(n+2)), for a sufficiently large fixed c_w, and add/subtract, Boolean operations, shifts, comparison, indexed memory and word multiplication/division. The input consists of two n-bit nonnegative integers in packed words. Output is their ordinary packed product. Input conversion and output packing operate on O(n/w) words, with a constant number of shifts and masks per digit. A product of two fixed-size multiword scalars costs O(1) instructions: split into half-word digits and use a fixed number of word products and carry updates. No full-input unit-cost product is used. + +Choose b=floor(log2(n+2)/8), after a fixed small-input cutoff, and let B=ceil(n/b). Write both inputs in radix 2^b with B digits. Let L be the smallest power of two at least 2B. Thus 2B<=L<4B, and L=Theta(n/log n). All digits and all exact convolution coefficients have O(log n) bits. + +## The explicit transform hypothesis + +Assume an actual uniform native cyclic Fourier program on L complex entries provides the following contract. For arbitrary complex inputs with real and imaginary parts of magnitude at most L^4 2^(4b), it computes either the unnormalized forward transform or the normalized inverse transform, with absolute output error at most delta, where + + delta = 1/(128 L 2^b), + +using O(log n)-bit fixed-point words and total work F(L,n). Its full root construction, address generation, table construction, precision operations and scratch clearing are included in F. A program proof must establish this hypothesis for the full arbitrary-input domain; a kernel identity or nonzero denominator check does not establish it. + +## Constructive assembly and proof + +Pad the digit arrays by zeros to length L. Apply the forward program to each, multiply corresponding approximate outputs using bounded-word complex arithmetic, and apply the normalized inverse program. Round every final real output to the nearest integer. Compute each pointwise product with absolute error at most delta. This is feasible using O(log n) fractional and integral bits. + +The exact forward outputs have modulus at most L 2^b. The difference between an exact product AB and the computed product is at most + + 2 L 2^b delta + delta^2 + delta < 1/16. + +The normalized inverse Fourier operator has infinity-to-infinity norm at most 1: each row has L entries of modulus 1/L. Thus its exact application to the product errors contributes less than 1/16 per output. The inverse program adds at most delta<1/16. The final real error is therefore less than 1/8, so nearest-integer rounding recovers every exact convolution coefficient. The imaginary part is discarded. The forward approximants, product approximants and inverse inputs all lie within the stated magnitude domain, including fixed rounding margins. + +There is no cyclic aliasing because L>=2B. Each coefficient is at most B(2^b-1)^2. A left-to-right radix carry pass uses the recurrence + + s_j = c_j + carry_j, + digit_j = s_j mod 2^b, + carry_(j+1) = floor(s_j/2^b). + +The carry is at most B 2^b by induction, up to a fixed harmless endpoint constant. Each s_j has O(log n) bits, and shifting/masking and O(1) word additions implement each step. Remaining carries terminate in O(1) additional digits because their size is O(log n). The carry and output packing cost O(L), including zero padding. + +Hence the complete native multiplication algorithm has cost + + O(F(L,n) + L), L=Theta(n/log n). + +The input-dependent word width is determined from n; the program and every fixed constant are independent of n. Fixed finitely many small inputs can use grade-school digit multiplication. + +## Native saving corollary and its exact scope + +If a concrete transform satisfies + + F(L,n) = O(L (log L)^(1-kappa) (loglog(L+4))^C) + +for fixed kappa>0 and fixed C, the multiplication cost is + + O(n/(log n)^kappa * (loglog(n+4))^C + n/log n). + +For every fixed 0=2 and alpha= P ln 2 + ln(32(X+1)), + a(2R-1) >= ln 2. + +These conditions can be selected by rationally certified squared comparisons and a doubling loop; no exact square-root instruction is needed. For any real center z, splitting the shifted lattice by its two sides and comparing successive terms gives + + sum_(|j-z|>R) exp[-a(j-z)^2] + <= 4 exp[-a(R-1)^2] / (1-exp[-a(2R-1)]). + +Consequently the periodic omitted tail, including the 1/alpha prefactor, is below 2^(-P-1), after the displayed fixed guard. This bound covers the infinite lattice; a finite s-times-largest-entry bound alone would not establish it. + +Take a power-of-two block width q=Theta(R). Partition output k and the extended source indices j into q-blocks. Retain a whole source/output rectangle whenever some pair lies within the radius R; each output block meets O(1) source blocks because both slopes are constant. The centered exponent on every retained rectangle is O(p). The exact diagonal-Hankel-diagonal identity in GAUSSIAN-CHIRP.md reduces each rectangle to one paid linear convolution of length O(q), with O(p) magnitude/accuracy guards and O(q log q) native word work. It is essential to evaluate the whole rectangle, rather than reintroduce an entrywise truncation mask. + +For the actual asymptotic line lengths, s is superpolynomial in p and hence s>2R. Centers lie in [0,s), so only the three source-period copies in [-s,2s) are needed for retained blocks. Access u_j by j mod s, using bounded-word arithmetic; this adds a constant number of boundary blocks. No scan of infinitely many copies is performed. The tail estimate above accounts for all omitted copies. + +There are O(t/q) rectangles. With O(log p) coefficient/precision guards and the ordinary paid exp/root preparation, the complete forward map costs + + O(t log p + n^eta polylog n) + +for any fixed small eta>0. Its working words have a sufficiently large fixed constant times log n bits. Prefactor multiplication, root generation, FFT kernel transforms, gathering, summation and periodic boundaries are included. This replaces the literal window-product count O(t alpha sqrt(p)), not merely its tape movement. + +## Actual inverse: O(s(1+p/alpha^2)) + +The source defines N=C T D, with q_i=floor((t/s)i+1/2), beta_i=(t/s)i-q_i and |beta_i|<=1/2. Its actual entries and the explicit tail proof are given in NATIVE-BANDED-INVERSE.md. Truncating to periodic distance h=O(1+sqrt(p)/alpha) gives row-sum error as small as 2^-P, because the tail is bounded by + + 4 exp[-pi alpha^2 h(h+1)]. + +The predecessor bound ||N-I||<2^(-alpha^2 theta) is eventually at most 1/8. Form the positive truncated coefficients using the paid exp table. The ordinary-band part is strictly row diagonally dominant; its factors have pivots uniformly bounded away from zero and no unbounded element growth. The periodic corner entries are a literal rank-at-most-2h update. The Woodbury Schur matrix has norm-distance at most 1/3 from identity, so its small inverse is itself computed by the same proved elimination, rather than supplied by an oracle. + +Band factorization, 2h right-hand solves, the small Schur factorization, matrix formation and final recovery together cost + + O(s(1+h)^2) = O(s(1+p/alpha^2)). + +The storage is O(s(1+h)). Accuracy P+O(log s+log h) bits suffices by the explicit rounding/residual proof. This remains O(log n). The external diagonal D and selector C cost O(t) extra work with the paid coefficient table. Shared exp/root preparation adds n^eta polylog n, once for the complete family. + +## Sum over the actual tensor axes + +Let T be the padded array size, d=Theta(p^epsilon), and alpha^2=Theta(d sqrt(p)), as in the source. There are T/t_i disjoint forward lines along axis i, and comparable source/target axis lengths. The sum of the forward costs is O(T d log p). The sum of the inverse costs is + + O(T[d + d p/alpha^2]) = O(T[d+sqrt(p)]). + +Including diagonal maps, selectors and setup, the complete Gaussian axis adapters therefore cost + + O(T[p^epsilon log p + sqrt(p)] + n^eta polylog n). + +A stronger actual reuse bound is available because all parallel fibers of one axis share the same s_i,t_i,alpha and hence the same matrix N. Factor the ordinary band part and prepare X=A0^-1 U and the small Schur factors once per axis. A new right-hand side then uses only two band triangular solves, a small Schur solve and Xz recovery, for O(s_i(1+h)) work. The Schur O(h^2) solve is covered because s_i>h. Across all fibers, the inverse contribution is O(T d(1+h)); the one-time factor work is O(d r(1+h)^2), where r=max(t_i)=2^{Theta(p/d)}=n^{o(1)}. With h=O(1+p^(1/4)/sqrt(d)), the complete stronger adapter bound is + + O(T[p^epsilon log p + p^(1/4+epsilon/2)] + + n^eta polylog n + n^{o(1)}). + +Numerical factor errors do not accumulate across independent right-hand sides; each satisfies the same explicit operator/residual error bound. All stored factors and Woodbury columns are paid for in O(d r(1+h)) preparation space. This is direct factor reuse of an explicitly stored input-independent map, not an oracle or an uncharged matrix inverse. + +This is a strict asymptotic native improvement over direct dense-window evaluation. The direct forward loops alone cost O(T d alpha sqrt(p))=O(T p^(3/4+3epsilon/2)); the new forward contribution is O(T p^epsilon log p). The inverse uses direct band elimination instead of repeated Neumann window applications. These are changes to concrete arithmetic algorithms, not a division of a TM-time estimate by word width. + +## Positive O(1)-amortized selected-bit metadata lemma + +Suppose a simultaneous layer selects d of the L address bits, with d=(log n)^epsilon for a fixed epsilon<1 and L=O(log n). For each of its O(log n) levels, prepare an array deposit[j], 0<=j<2^d, mapping a compressed selected-bit index to its physical-bit submask. Explicit construction costs O(2^d poly(d,log n)) per level, and total preparation is n^o(1). Any fixed rational-kernel template's permutation and quarter-phase tables on these d bits can be prepared with the same type of paid bound. + +Let M be the complementary-address mask. Its submasks can be enumerated in increasing order using + + next = (current-M) & M + +in modular word arithmetic, stopping on return to zero. For each complementary base and each j, base OR deposit[j] is its physical address. Every input point belongs to exactly one such pair, and all T gathers/scatters take O(T) word instructions. No PEXT, PDEP, popcount oracle or d-step extraction is required per point. + +For a separately applied tensor twiddle diagonal, scan coordinates in natural mixed-radix odometer order. Maintain the sum of the per-axis exponent contributions; only changed axes are updated. Since every radix is at least 2, the total number of axis changes is O(T+d), so the metadata work is O(T+d), with bounded-word arithmetic and the paid root-power queries. This lemma applies when the phase is the actual separated sum of per-axis contributions, as in the common tensor level; it does not assert the same cost for an arbitrary nonseparable phase function. + +## Verification scope and next assembly + +The real ordinary-band bounded-word reference implementation band_solve.py was independently rerun: 20 systems, 208 rows, maximum raw arithmetic width 193 bits under a 256-bit budget; singular and overflow controls passed. Fraction is confined to its external exact residual oracle. The implementation uses exact product accumulation in substitution and fixed-grid elimination. Its operation counter is an instrumented model, not an emitted RamStd instruction certificate. The periodic real Woodbury implementation periodic_solve.py was also independently reviewed and rerun: 3 periodic systems, 9 right-hand sides reusing the same prepared factors, 123 rows, maximum raw width 192 bits. Its preparation constructs the ordinary factors, the 2h inverse-update columns and small Schur factors exactly once; each new right-hand side performs O(s h) reuse work with no refactor. These periodic fixtures evaluate exact residuals with an external Fraction oracle. The actual Gaussian N is real, so a complex input is handled by separate real and imaginary right-hand solves, doubling only the constant. The general complex-matrix elimination proof remains informal; it is not an executable path covered by these real fixtures. + +The rational closed Gaussian/Walsh kernel from #956 also has the positive native precision lowering proved in NATIVE-STABILITY.md. Choosing epsilon=1/8 suggests a conservative full native saved-exponent assembly once its tensor caller, arbitrary twiddles, record products, cyclic embedding, precision and recovery are all connected and replayed. That complete assembler is not emitted here. In particular these component results are not a new best kappa and are not a completed proof of a faster full cyclic transform. + +Credit: upstream Gaussian resampling and inverse estimates are predecessor results with their Harvey-Hoeven references; the rational kernel is #956 with its stated lineage; chirp convolution, FFT, diagonal dominance and Woodbury have classical precedents. The contribution here is the explicit native finite-word construction, quantitative conditioning, paid preparation and actual Gaussian caller recount. Prepared for Dugongue with substantial OpenAI Codex assistance. + +Additional root integration: periodic_solve.py now replays real Woodbury systems with reusable factors. gaussian_inverse.py constructs the actual positive Gaussian N entries using a paid Machin pi series and signed exp tables, then solves three periodic systems (55 rows), with reference residual below 2^-100. gaussian_block.py evaluates 12 concrete Gaussian rectangles (45 outputs), checked below 2^-60. These are finite operational-model fixtures, not a full cyclic-transform replay; finite Decimal sums are independent validation only. diff --git a/research/native-word-fft/NATIVE-STABILITY.md b/research/native-word-fft/NATIVE-STABILITY.md new file mode 100644 index 0000000000..86f984363c --- /dev/null +++ b/research/native-word-fft/NATIVE-STABILITY.md @@ -0,0 +1,118 @@ +# Direct bounded-word lowering of the closed rational Gaussian kernel + +This is a native arithmetic construction. It neither calls a Turing machine nor imports a TM simulation or multiplication routine. Its purpose is to turn an existing exact-scalar Gaussian/Walsh provider into a bounded-word algorithm with an explicit accuracy guarantee. + +## A literal sufficient class + +Consider a fixed finite provider template with: +- arbitrary independent live arrays and explicitly cleared owned scratch; +- fixed rational complex scalar ADD/SCALE operations, COPY, CLEAR, and exact monomial coordinate/quarter-phase operations; +- recursive children that are exact unitary Gaussian kernels, possibly on disjoint fibers and conjugated by monomial operators; +- a complete exact endpoint equal to the desired unitary kernel on every live input, with owned scratch erased. + +Read every coefficient from the literal finite template. Replace each rational scalar operation by its signed integer numerator/denominator calculation followed by rounding to a fixed binary grid. Multiplication/division here is on a fixed number of O(log n)-bit words, not on the full input. Nonzero fixed denominators are explicit integers. Metadata, address calculations and phase choices are integer computations. + +This class includes the rational Gaussian kernel in PR #956, whose nonchild coefficients are fixed with denominators supported on 2 and 3 and whose monomial phases are fourth roots of unity. It does not automatically include that PR's separate all-length cyclic Fourier compiler, which prepares additional input-length-dependent scalars through Newton/Toeplitz constructions. + +## Prefix norms are derived, not assumed + +For one scalar gate x<-x+a*y, the Euclidean operator norm is at most 1+|a|. For x<-a*x it is at most max(1,|a|), including the untouched coordinates. COPY has norm at most sqrt(2), CLEAR has norm 1, and coordinate/quarter-phase monomials have norm 1. These statements hold on all arbitrary inputs, not only on a chosen test vector. + +Replace |a| by the explicit rational upper bound |Re(a)|+|Im(a)|. Multiply these elementary bounds over the finite template; enlarge the result to an integer K>=2. This is an effective finite calculation from actual coefficient bytes. It bounds every prefix and suffix of that template when recursive children are treated as their exact unitary operators. + +A coefficientwise pass acts as a direct sum of the same fixed bank operation, so its norm bound does not grow with array length. Disjoint completed fiber calls form a direct sum of unitary maps. Completed owned-scratch macros have the extended map U_live direct-sum 0_owned after their initial clearing, hence norm at most 1. A primitive prefix inside one such macro is covered by the same fixed template bound. + +Classical base butterflies are grouped into completed unitary two-point kernels. Thus an arbitrarily long sequence of completed classical rounds does not multiply their prefix bounds: a partial primitive prefix lies inside only one current butterfly, while all completed rounds are unitary. + +Let D(e) be the maximum recursive call depth at rank e. An expanded prefix or suffix is inside at most one unfinished child at each depth. All its other completed children are unitary. Induction therefore gives the actual bound +$$ +\|\text{expanded prefix}\|_2,\;\|\text{expanded suffix}\|_2 +\le K^{D(e)+1}. +$$ +The exponent is depth, not total work. It does not hide an unproved stability oracle. + +## Explicit roundoff theorem + +Suppose the exact scalar count is at most +$$ +M(e)\le A\,2^e(e+1)^\theta +$$ +for fixed A>=1 and 0=0. Every stored signed fixed-point numerator uses at most I+P+O(1) bits. Fixed rational numerators/denominators add only a template-dependent constant. Choose the physical word-width coefficient large enough for the exact intermediate numerator products, sums and divisions, or use a fixed number of limbs. This yields O(1) word instructions per scalar update. + +When e,b,B=O(log n), D(e)=O(log(e+2)), and A,K,Z are fixed before n, I+P=O(log n). Thus the finite-word program has the same asymptotic arithmetic operation count as the exact-scalar provider, with a computable constant and no TM or bit-cost-machine call. The cutoff can be enlarged until fixed template constants fit; below it the program uses direct unitary butterflies. + +## Instantiation of the actual closed provider + +PR #956 is pinned at fa4df77c9dfdbac93368d95dcb361e69b3b65464. Its optimized source uses m=73, largest child rank 48, productive stock W=58,899, and a complete exact profile. Its own all-input recurrence is +$$ +t(e)\le\sum_r(H_r/W)t(rf(e)) ++A\,[a+m+2\lceil\log_2(e+2)\rceil+1], +$$ +where t(e)=T(e)/2^e, f(e)<=e/73 and the explicit integer-root certificate funds theta=1023/1024. + +The actual contraction gives +$$ +D(e)\le 1+ +\left\lceil \frac{\log(e+1)}{\log(73/48)}\right\rceil. +$$ +All source bank/owner constants are fixed before e. The elementary template bill supplies a finite upper bound for the number and coefficient sizes of its nonchild passes, from which K above is computable. In particular this does not replace a varying exact scalar of unknown condition by a word operation. + +Therefore this pinned rational native Gaussian/Walsh provider admits a bounded-word approximation program of work +$$ +O\bigl(2^e(e+1)^{1023/1024}\bigr) +$$ +and arbitrary requested error 2^-B with the explicit precision formula above. For the unnormalized Walsh conversion H_e=(1-i)^e D_e C_e D_e, first request kernel error 2^{-(B+ceil(e/2)+1)} and add ceil(e/2)+1 integral bits for the final scaling. Since |(1-i)^e|=2^{e/2}, the final error is at most 2^-B. These additions remain O(e). Quarter-phase D_e is exact. The power (1-i)^e is an exactly representable Gaussian integer and can be prepared by O(e) bounded-word updates, absorbed by the displayed work bound. No integer-multiplication theorem follows merely from this kernel statement. + +The same argument applies to a different *actually constructed* fixed rational provider and its proved recurrence; it does not assert a better profile or select unknown scalar coefficients. In particular, a critical-root compiler word must first have an actual admitted algebraic finite-word construction. Removing tape tariffs alone is not that construction. + +## Boundary for cyclic transforms and multiplication + +A separate cyclic transform requires explicit algorithms for roots, scalar preparation, conditioning, numerical errors, and convolution/recovery. Nonzero prepared denominators alone do not imply polynomially bounded inverse magnitudes. The exact-scalar PR #956 cyclic section does not supply the quantitative bounds needed here merely by proving nonzeroness. + +Root/phase tables and ordinary bounded-word FFT base multiplication can instead be constructed directly, but their complete caller and accuracy costs must be provided. No base multiplication TM is permitted in such an assembly. The nine/54-row TM outer ledger cannot simply have all costs divided by word width: each retained operation must be rebuilt as array/word arithmetic and re-invoiced. + +This is a positive finite-word native kernel theorem and a concrete assembly interface. It does not establish a new integer-multiplication exponent or worldwide priority for standard roundoff analysis. + + +## A quantitative bound for the actual shear preparation + +One denominator family in the actual cyclic compiler is already quantitatively safe. Write mu=a+ib with a,b real, and put t1=1+mu*conjugate(mu), t2=mu-t1. Then +$$ +t_1=1+a^2+b^2\ge1,\qquad +\operatorname{Re}(t_2)=-[(a-1/2)^2+b^2+3/4]\le-3/4. +$$ +Consequently |1/t1|<=1 and |1/t2|<=4/3, uniformly in mu. This is stronger than nonzeroness and is a direct sufficient stability bound for these two divisions. If mu is approximated within delta<=1 and |mu|<=U, then the induced t1 perturbation is at most (2U+1)delta; t2 adds another delta. Provided the t2 perturbation is at most 3/8, its reciprocal perturbation is at most (32/9) times that perturbation, by the identity 1/z-1/z'=(z'-z)/(zz'). Thus this shear preparation is polynomially conditioned whenever U is polynomially bounded. The Newton c-inverses and complete cyclic circuit still require their own quantitative analysis. Credit: ram_cost identified the completion-of-squares argument during this session. + +## Polyloglog overhead still retains a strict saving + +For completeness, a direct finite-word compiler need not retain a constant number of words to preserve an exponent saving. Suppose an actually constructed native transform on N data words has exact arithmetic work O(N (log N)^{1-kappa}) with fixed kappa>0, and an explicit numerical analysis and implementation bound every scalar update, preparation and packing overhead together by (loglog(N+4))^C for a fixed C. Then the finite-word cost is +$$ +O\bigl(N(\log N)^{1-\kappa}(\log\log(N+4))^C\bigr) +=O\bigl(N(\log N)^{1-\kappa'}\bigr) +$$ +for every fixed 0=0. This is a positive absorption theorem; it does not assert that the unresolved Newton/Toeplitz preparation meets its explicit hypothesis. It is useful when a separately derived precision budget is O(log N (loglog N)^c): elementary schoolbook limb multiplication and long division on (loglog N)^c limbs have polynomially many limb operations, and therefore satisfy this overhead condition. No unbounded exact-scalar instruction or TM multiplication routine is invoked. + +## Credit and evidence + +The closed exact provider, arbitrary-input identity, ownership, recurrence, finite-profile certificate and explicit pass constants are predecessor results from #956 and its #390/Sussman/P193, native, frame/bank, bridge and normalization lineage. Prepared for Dugongue with substantial OpenAI Codex assistance. The roundoff proof above is informal; no new Lean or GPU process was run. Source inspection is not independent acceptance of the predecessor's complete theorem. diff --git a/research/native-word-fft/NOTICE.md b/research/native-word-fft/NOTICE.md new file mode 100644 index 0000000000..a4497310e9 --- /dev/null +++ b/research/native-word-fft/NOTICE.md @@ -0,0 +1,37 @@ +# Attribution and license + +New exposition and executable models: Dugongue, with substantial OpenAI Codex +assistance from the root worker and three parallel Codex workers. Apache-2.0; +see LICENSE. Independent agent review means an additional AI review of the +stated source bytes, not human peer review or formal kernel certification. + +The actual Gaussian resampling identity, the definition N = C T D, and its +near-identity estimate are predecessor facts from OpenAI/math's manuscript, +retained in this repository at upstream/build/sections/07-resampling.tex, with +its citations to David Harvey and Joris van der Hoeven's integer-multiplication +work. The original OpenAI source/license attribution in the repository remains +in force. No retained upstream file is modified by this package. + +Classical algorithmic antecedents: Cooley–Tukey FFT; Bluestein's chirp filtering +(DOI 10.1109/TAU.1970.1162132); Rabiner–Schafer–Rader's chirp z-transform +(DOI 10.1002/j.1538-7305.1969.tb04268.x); Woodbury low-rank inversion; +row-diagonal-dominance elimination; Machin's formula; and standard digit +convolution/carry recovery. Fürer, arXiv:1402.1811, already establishes a linear +logarithmic-word RAM multiplication baseline. None is presented as our discovery. + +The rational native provider comes from artuskg's PR #956 at +fa4df77c9dfdbac93368d95dcb361e69b3b65464. Retain its credit to Dugongue's PR #390, +its retained contributors, Jacob Sussman's bridge/frame/banking/RAM/Fourier work +and P193 supplier; PRs #363/#373/#375, #377/#381, #397, #401 and #430; and +OpenAI/math's preceding exact-width, sector, CRT/chirp and root-order sources. +PR #647 (artuskg) is a separate exact-scalar construction and remains credited. +The finite-word lowering here does not assert its unproved cyclic conditioning. + +Comparison only: hexapode's pending PR #4535, and the selected PR #4359 lineage. +Their exact heads and scope are in OVERLAP.md. Dugongue's PRs #4538/#4541 supplied +previous Word RAM compiler/transfer packages; this package does not call their +TM simulator. Existing public attribution and source notices remain authoritative. + +Inspected Problem 1.8 limb/inversion proofs, Gaussian expansion notes and a MUB6 +Machin/rounding note are listed by source hash in RECOVERY-PRECEDENTS.json. +Those are source-inspection precedents, not a new formal proof or novelty audit. diff --git a/research/native-word-fft/OVERLAP.md b/research/native-word-fft/OVERLAP.md new file mode 100644 index 0000000000..e8cb2ed4f4 --- /dev/null +++ b/research/native-word-fft/OVERLAP.md @@ -0,0 +1,43 @@ +# Frozen overlap and frontier checkpoint + +Canonical main: f626db5c3a3d781d8f7695038776c621651f75e1. +Fetched again during this session; guidance, priorities and selected record were +unchanged. Selected conditional κ remains 23/10000 = 0.0023, sourced through +PR #4359 and the lineage in certificates/selected-result.json. This package +changes no selected certificate or frontier record. + +Pending comparison: PR #4535, hexapode, +a6d7237f697185e857090e174c23c1cb22c6563c, advertises κ = 0.003020165. +We do not claim to surpass it, and changing model does not establish a better +TM bound. The current package is positive native arithmetic development. + +Previously inspected native precedents: +- PR #956, artuskg, fa4df77c9dfdbac93368d95dcb361e69b3b65464: + research/pr390-heterogeneous-provider-closure/PROOF.md and its NOTICE.md. + Exact-scalar/index RAM kernel and inherited all-length Fourier result; + scalar heights there are unrestricted. Our stability proof is scoped to the + actual fixed-rational unitary kernel, not all of its cyclic coefficient DAG. +- PR #647, artuskg, 3b3fda9c602fbe4b45d02d1ead0b757b3c94f52c: + distinct exact-scalar E8/Fourier predecessor; finite-word conditioning is not + inferred from its nonzero denominators. +- PR #4541, Dugongue, 46ca24a88b5194725a82c02003aac8aea0c6f53a: + packed TM-to-Word-RAM transfer, conditional O(n/log^κ n). The present native + algorithms replace TM calls; they do not reuse that simulation as arithmetic. +- PR #4538: previous explicit Word RAM bridge, not a native fast multiplier. + +Public updated PR titles/heads were checked again through #4551. Canonical main remained unchanged. PR #4549 supplies a noncommuting-controller support proof; #4550 is a halo-chain disposition; #4551 audits cancellation. They are title-level triage, not mathematical premises of our native algorithms. PR #4536 independently audits #4535 and does not advertise a new bound. The #4548 publication-guard head changed to ae19b6868d5981fbacdf6950a226d969e603b78d. PRs #4542–#4546 +cover physical-transport/headroom/cancellation audits, not these native Gaussian +arithmetic constructions. PRs #4547/#4548 concern portable replay/publication. +Titles are triage evidence; this is not a full mathematical review of their diffs. + +The repository overlap helper was run on bit-new-register-model. Its sparse local +coverage omitted 5,855 HEAD note paths, and remote Discussion refresh hit an API +limit. A local cache-directory write was also unavailable. Direct public PR API +reads supplied the listed heads. No absence-of-overlap or historical-priority +claim follows from that incomplete search. + +The federated DAG search and inspected limb, inverse, Gaussian and Machin notes +are recorded as source-only precedents. We did not search every registered proof +or freshly kernel-check them. The single canonical Lean tree/live DAG was left +unchanged. The resulting scope is explicitly changed from the old register card +because the user assigned native Word RAM replacement. diff --git a/research/native-word-fft/PHASE-TABLES.md b/research/native-word-fft/PHASE-TABLES.md new file mode 100644 index 0000000000..2dab8f3e6d --- /dev/null +++ b/research/native-word-fft/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/native-word-fft/PROOF.md b/research/native-word-fft/PROOF.md new file mode 100644 index 0000000000..c7b4777b05 --- /dev/null +++ b/research/native-word-fft/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/native-word-fft/README.md b/research/native-word-fft/README.md new file mode 100644 index 0000000000..3c77ae69e3 --- /dev/null +++ b/research/native-word-fft/README.md @@ -0,0 +1,101 @@ +# Native Word RAM Gaussian adapters + +This package removes TM subroutines from the arithmetic components developed +here. Its main contribution is two concrete native algorithms for the actual +Gaussian resampling maps, with paid preparation, precision and recovery bounds. +It does **not** establish a new best κ or a complete faster cyclic transform. + +1. **Chirp blocks:** a centered Gaussian rectangle is a diagonal–Hankel–diagonal + product. One explicitly paid convolution evaluates it in O(q log q) Word RAM + operations. Whole-rectangle blocking gives the actual forward map in + O(t log p) operations, plus shared scalar preparation, instead of visiting + Θ(α√p) weights for every output. +2. **Reusable banded inverse:** the actual periodic near-identity matrix has an + explicit Gaussian tail. Truncate to h = O(1 + √p/α), factor its ordinary band, + and handle the corners with a rank-at-most-2h Woodbury correction. Uniform + pivot and Schur bounds give O(s h²) paid preparation and O(s h) per new + right-hand side. Nothing assumes an inverse oracle. + +For T entries and d = Θ(p^ε) tensor axes with α² = Θ(d√p), reuse across the +actual parallel fibers yields Gaussian adapter work + + O(T[p^ε log p + p^(1/4+ε/2)] + n^η polylog n + n^o(1)), + +where p = Θ(log n) and η > 0 is a fixed small setup exponent. See +[NATIVE-RECOUNT.md](NATIVE-RECOUNT.md) for the full hypotheses and cost proof. +This improves the adapter's arithmetic count; it does not improve the complete +multiplication exponent without the remaining transform integration. + +The package also supplies a direct exact O(n) packed Word RAM multiplier with +an all-input rounding proof, generated roots, signed exponential tables, paid +constant-cost phase queries, and a finite-precision rational Gaussian kernel. +The O(n) multiplication baseline is classical, not a new bound: see +[Fürer, 2014, Theorem 4.1(a)](https://arxiv.org/abs/1402.1811). + +## Model and files + +All native algorithms use one sufficiently large fixed constant times log n +bits per word, with ordinary bounded-word MUL/DIV, shifts, masks and indexed +memory. Full n-bit operand multiplication is confined to independent test +oracles. Every setup cost, including roots, coefficients and fixed-kernel FFTs, +is paid. Large word constants simplify the proof; no practical speedup is claimed. + +| Result | Argument | Executable evidence | +|---|---|---| +| Gaussian block and forward map | GAUSSIAN-CHIRP.md; NATIVE-RECOUNT.md | chirp_check.py; gaussian_block.py; native_convolution.py | +| Actual periodic Gaussian inverse | NATIVE-BANDED-INVERSE.md | band_solve.py; periodic_solve.py with factor reuse; gaussian_inverse.py constructs actual N entries | +| Direct exact O(n) native multiplier | PROOF.md; FFT-ERROR.md | native_fft.py, including a 65,535-bit case | +| Generated roots and exact carries/CRT | ROOTS-RECOVERY.md | root_seed.py; recover.py | +| Paid exp and phase tables | NATIVE-EXP-TABLE.md; PHASE-TABLES.md | exp_table.py; phase_tables.py | +| Rational Gaussian/Walsh kernel lowering | NATIVE-STABILITY.md | Informal source-specific proof; no new kernel replay claimed | +| Sufficient full-transform caller | NATIVE-OUTER-ASSEMBLY.md | Informal assembly theorem; faster caller not supplied | + +These are informal universal arguments and finite executable controls, not Lean +proofs, emitted RamStd programs, human peer review or a certificate of the full +multiplication theorem. Python integers emulate explicitly bounded magnitudes; +Decimal and Fraction appear only in independent references. The real-grid +Gaussian executable covers a concrete block family; the periodic inverse and +actual Gaussian coefficient path have separate real-valued fixtures. General +complex-matrix elimination remains proof-level; real Gaussian N acts on complex +data by two real solves. Source bytes are authenticated before any replayed module executes. + +## Reproduce + +Python 3.12, standard library only; commands work on Windows and Linux: + + python3.12 scripts/research_package.py --package research/native-word-fft + +Standalone, including an authenticated Git-free export: + + python3.12 -I -O -B research/native-word-fft/verify.py + +VERIFY.json registers the bounded replay. TASK.json preserves the assigned +scope change; the old register-model card is provenance only. SOURCES.json +pins canonical main and the selected record. Additional source and overlap +coverage is documented in OVERLAP.md, KERNEL-PRECEDENTS.json and +RECOVERY-PRECEDENTS.json. Source inspection is not kernel verification. + +## Open next questions + +Connect the native adapters, rational kernel, phases, synthetic polynomial +transforms, exact record products and recovery in one complete cyclic caller. +Freeze every normalization and one common precision/memory bound, then check the +full paid recurrence before reporting any new κ. Expand the actual Gaussian fixtures into the full asymptotic parameter and +normalization regime, and connect the reusable periodic factors to every axis +of that caller. TASK.json gives +concrete first steps and their relevance. + +## Credit + +Prepared by Dugongue with substantial OpenAI Codex assistance, using the root +worker and three parallel workers. The Gaussian identities and inverse estimates +come from OpenAI/math's retained manuscript and Harvey–van der Hoeven. Chirp +convolution, FFT, Woodbury and diagonal dominance are classical precedents. +The rational provider is artuskg's PR #956, building on Dugongue's PR #390 and +Jacob Sussman's bridge/frame/RAM/Fourier lineage; NOTICE.md lists the relevant +PRs. PR #647 is a distinct exact-scalar precedent, not a proved finite-word +caller. PR #4535 supplies the current pending comparison; PRs #4538/#4541 are +our previous Word RAM compiler/TM-transfer work. Their TM simulation is not used +in these native constructions. Problem 1.8 limb/inversion modules and earlier +Gaussian/Machin proof notes informed the bounds; their inspected hashes are +recorded without claiming fresh Lean audits. No worldwide novelty claim is made. diff --git a/research/native-word-fft/RECOVERY-PRECEDENTS.json b/research/native-word-fft/RECOVERY-PRECEDENTS.json new file mode 100644 index 0000000000..62864a07f3 --- /dev/null +++ b/research/native-word-fft/RECOVERY-PRECEDENTS.json @@ -0,0 +1,22 @@ +[ + { + "source": "MpLimbs20261007.lean", + "sha256": "cbbb9b64e17dd9059378e92819c7da492a1643a6505ea35d24e486d8d6798600", + "evidence": "source inspection only; no new kernel audit" + }, + { + "source": "Problem18IMPPrimeInverse20261001.lean", + "sha256": "01f2247cea102e321304ebe22aa6cfa9251edd448525ef131ef20f354c54543c", + "evidence": "source inspection only; no new kernel audit" + }, + { + "source": "PR3877 Gaussian expansion headroom; Dugongue, Mason Hale interval lineage", + "sha256": "2af6ee0d075d7b5dbad5c58314418a84e57d54dd15fc77bf81ea016666aaa95a", + "evidence": "source inspection only; not independently validated" + }, + { + "source": "OpenAI MUB6 preprint, 2026-09-24, rounding section", + "sha256": "b3c819a70e33bace099d89e214471739a84aba20cf062c03fed67ff2053ae005", + "evidence": "source inspection only; not independently validated" + } +] diff --git a/research/native-word-fft/ROOTS-RECOVERY.md b/research/native-word-fft/ROOTS-RECOVERY.md new file mode 100644 index 0000000000..5791526c61 --- /dev/null +++ b/research/native-word-fft/ROOTS-RECOVERY.md @@ -0,0 +1,223 @@ +# Native root seeds and exact product recovery + +Prepared by Dugongue with substantial OpenAI Codex assistance. Apache-2.0. +Every algorithm here acts directly on bounded words, arrays and arithmetic. +There is no Turing-machine primitive, simulator, conversion or multiplication +black box on an entire variable-length input. Ordinary MUL/DIV on word-sized +operands are allowed, with explicit sufficient width below. + +## 1. A generated primitive-root seed with certified precision + +For every integer D>=1 and requested p>=1, root_seed.py constructs Gaussian +dyadic coordinates x/2^p,y/2^p satisfying + + max(|x/2^p - cos(2pi/D)|, |y/2^p - sin(2pi/D)|) < 2^-p. + +Take P=8(p+16), S=2^P, K=128(p+1). It suffices that every word has at least +2P+20 bits and that D and all array/control addresses fit in one word. Signed +integers are a sign bit and an unsigned magnitude; signed arithmetic expands +to bounded comparisons, ordinary ADD/SUB and word MUL/DIV. This is a uniform +finite loop algorithm, not a stored transcendental constant or root oracle. +For p=Theta(log N) and D<=N^C with fixed C, the sufficient word width remains +Theta(log N). Constants are deliberately conservative. + +### Construct pi rather than assuming it + +Use Machin's identity pi=16 arctan(1/5)-4 arctan(1/239). Its identity can be +verified without a numerical pi supplier: tan(2 arctan(1/5))=5/12 and + tan(4 arctan(1/5))=120/119. The tangent subtraction formula then gives + tan(4 arctan(1/5)-arctan(1/239))=1; the angle is in (0,pi/2), so it is pi/4. + +For q=5 and 239, sum the alternating terms + + (-1)^j floor(S / ((2j+1) q^(2j+1))) + +while q^(2j+1)<=S. Generate successive powers by multiplication by the fixed +integer q^2. There are at most P+1 terms. Each floored term differs from its +exact scaled term by less than 1, and the alternating remainder after stopping +is at most 1/S. Consequently the constructed integer A obeys + + |A/S - pi| <= 20(P+2)/S. + +Every denominator is at most S(2P+3); the next q-power is at most 239^2 S. +These are bounded O(P)-bit integers. Compute the integer angle t=floor(2A/D). +Then 0<=t<=8S and + + |t/S - 2pi/D| <= [40(P+2)+1]/S. + +No search for a prime or primitive root is involved. + +### Fixed-point complex Taylor evaluation + +Start the Taylor term and sum at (S,0). For k=1,...,K, replace the current +term (a,b) by + + (trunc(-b*t/(S*k)), trunc(a*t/(S*k))), + +and add it to the sum. Truncation is toward zero, implemented on sign/magnitude +words. Let theta=t/S. Each coordinate's rounding error in term k satisfies +E_k <= (8/k) E_(k-1)+1/S. Expanding this recurrence and using +j!/k!<=1/(k-j)! yields E_k=1. The cosine and sine +functions are 1-Lipschitz on the real line. Combining the angle, rounding and +tail errors gives the explicit coordinate bound + + E_internal <= [40(P+2)+1+4096K+32]/S. + +The displayed numerator times 2^(p+1) is strictly less than S. Finally round +the two fixed-point sums to nearest at p fractional bits. This adds at most +2^(-p-1), proving the stated strict 2^-p error. + +Term and partial-sum magnitudes are below 2^14 S: the exact absolute Taylor +mass is below exp(8), and the accumulated rounding envelope is smaller than S. +Products with t are below 2^(2P+17), denominators and all Machin intermediates +are smaller still. Thus 2P+20 bits is a sufficient native arithmetic bound. +The implementation charges a conservative + + 100(P+1)+40(K+1) + +word-operation envelope. Both P and K are O(p), so setup is O(p) word +operations and O(1) live arithmetic words, besides caller storage. It is not +necessary to form factorials or unbounded exact rational denominators. + +## 2. Rounding and CRT recover exact convolution coefficients + +Let two nonnegative base-beta digit arrays have lengths a,b, beta=2^h>=2. +Every digit is below beta. Their ordinary linear convolution has length +a+b-1 and coefficients + + c_j = sum_(i+k=j) u_i v_k, + 0 <= c_j <= C=min(a,b)(beta-1)^2. + +**Rounded numerical route.** If a computed dyadic real value Y_j/2^P satisfies +|Y_j/2^P-c_j|<1/2 for every j, nearest-integer rounding recovers every c_j +exactly. Range-checking 0<=c_j<=C is useful, but does not itself prove the +error premise. The numerical transform must establish that premise from its +actual precision and prefix-norm bounds; it cannot be inferred from good +sample outputs. The recovery loop costs O(a+b) word operations. Its input +integers fit when 2^P(C+1), plus fixed rounding slack, fits in a word. + +**Exact modular route.** For distinct certified primes p_1,...,p_s, let +Q=product p_i. If Q>C, the canonical residue in [0,Q) is c_j itself. For +signed coefficients |c_j|<=C, use Q>2C and the representative in (-Q/2,Q/2). + +The supplied CRT constructor checks each prime by literal trial division and +constructs inverses of preceding products by modular repeated squaring with +exponent p_i-2. Primality and inverse work is charged setup. Trial division +costs O(sum sqrt(p_i)) word operations; modular inversion costs +O(sum log p_i). This is a constructive positive algorithm, but these setup +costs are not asserted to fit a fast FFT bound. A fast application must provide +its own fully costed modulus construction or choose the numerical route above. + +For each coefficient maintain its reconstructed x in [0,M), initially x=0, +M=1. The precomputed next modulus data give + + t = ((r_i - (x mod p_i)) mod p_i) * inv(M mod p_i) mod p_i, + x <- x + M*t, M <- M*p_i. + +Compute the subtraction as r_i+p_i-(x mod p_i), then reduce; it is nonnegative. +Induction proves all processed congruences and 0<=x> h. + +If carry<=min(a,b)(beta-1), then + + t <= min(a,b) beta(beta-1), + +and the next carry obeys the same bound. Thus this quantity, together with +buffer addresses, is a sufficient word bound. No large input-sized integer +is ever multiplied by the algorithm. + +The induction invariant is + + sum_(j