diff --git a/research/pr4554-gaussian-chirp-check/.gitattributes b/research/pr4554-gaussian-chirp-check/.gitattributes new file mode 100644 index 0000000000..b263cf196f --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/.gitattributes @@ -0,0 +1 @@ +PR4554-GAUSSIAN-CHIRP.md -text diff --git a/research/pr4554-gaussian-chirp-check/PR4554-GAUSSIAN-CHIRP.md b/research/pr4554-gaussian-chirp-check/PR4554-GAUSSIAN-CHIRP.md new file mode 100644 index 0000000000..4e31f30215 --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/PR4554-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/pr4554-gaussian-chirp-check/README.md b/research/pr4554-gaussian-chirp-check/README.md new file mode 100644 index 0000000000..931fa25ced --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/README.md @@ -0,0 +1,56 @@ +# Independent check of #4554: stable Gaussian block maps by one paid convolution + +**Task:** independent-moment (independent source-bound coverage). **Base:** canonical main `f626db5c`. **No new κ.** + +Dugongue's #4554 (head `4f46c6a9`) evaluates Gaussian resampling blocks `W[x,y] = exp(−a(ux − vy + c)²)` with a diagonal–Hankel–diagonal factorisation, so each block costs one linear convolution. The supporting results are: +- conditioning after centring each block (A); +- a precision rule (C) that implies the block error bound (B); +- a blocking geometry (Theorem 3). + +This package checks `GAUSSIAN-CHIRP.md`, vendored as `PR4554-GAUSSIAN-CHIRP.md` (blob `4e31f302`, SHA256-pinned). Nothing is imported from #4554's `chirp_check.py`, `gaussian_block.py` or `native_convolution.py`. It follows my #4556–#4559 on the same PR. + +## Results + +All checked claims are confirmed. + +1. **Theorem 1.** + - `log D[x] + log E[y] + log H[x+y] = −a(ux − vy + c)²` holds as a polynomial identity in ℚ[a,u,v,c,x,y]. + - `Σ_y W[x,y]Z[y] = D[x]·(H ∗ reverse(E·Z))[x+l−1]` holds, with an explicit **linear** convolution, at 60 digits on 40 random rectangles (`m, l ≤ 7`). +2. **Conditioning (§2).** With `Cp := a(|u|q+|v|q+|c0|)²`, every exponent of D, E and H is at most `Cp` in absolute value. This is checked exactly in rationals on 300 random parameter sets, with signed `u`, `v`, `c0` and `q ≤ 24`. +3. **Blocking geometry (§4).** Brute force over 120 random rational geometries (3,522 retained blocks): + - each row block keeps at most `K = ⌈(u(q−1)+2R)/(vq)⌉ + 3` column blocks; + - every retained block has `|c0| ≤ R+(u+v)(q−1)`; + - (D) implies (A). + + The tail rule `aR² ≥ (p0+⌈log₂(sX+1)⌉+2)ln 2` keeps the omitted row mass below `2^−p0/2` (116 cases). +4. **Precision rule.** (C) implies `(B) ≤ 2^(−p0−4)` in all 576 parameter cases, using only `log₂e < 2`. +5. **End-to-end block.** I wrote my own fixed-point evaluator at the `P` that (C) prescribes: + - chirps and twiddles on the `2^−P` grid; + - nearest-rounded complex products; + - linear convolution through a radix-2 FFT of size `N0 = 4q`, then division by `N0`; + - selection at `x+l−1` and multiplication by `D`. + + Against a direct 120-digit sum, all 6 runs (`q ∈ {8, 16}`, with the conditioning scale `p` up to 52) have error far below `2^−p0` and within (B). Observed errors are about `2^−331` at `P = 512` and `2^−219` at `P = 234`. + + **Observation.** (C) is valid but conservative. The observed precision loss is about 180 bits at `p = 52`, against the roughly 484 bits that (C) reserves above `p0`. The `8⌈C⌉p` term, sized for `M⁴`, dominates the rule. +6. **Binding.** 7 statements are bound verbatim to the vendored text. + +## Verify + +```text +python3 -B research/pr4554-gaussian-chirp-check/verify.py +``` + +This runs in a few seconds with the standard library only. + +## Limitations + +- **Exp table.** The exp table that generates `D`, `E` and `H` is checked separately in #4559. Here the chirps come from a decimal reference rounded to the grid, which matches the `2^−P` interface of (B). +- **Cost.** Theorem 3's operation count is not re-counted. Only its geometric inputs, `K`, `|c0|` and (D), are checked. +- **Finite checks.** These are finite and random checks of universal statements whose proofs are written. + +## Open next questions + +- A sharper precision rule than (C). The observed loss suggests that something like `2Cp·log₂e + O(log q)` guard bits suffice. First step: redo the (B) derivation, tracking `M` powers separately for the forward and inverse transforms. + +Credit: Dugongue (#4554); the Harvey–van der Hoeven Gaussian resampling sources; classical chirp/Bluestein convolution. Prepared by Ryanson Jonathan with substantial Claude Code (Claude Opus 5.5) assistance. diff --git a/research/pr4554-gaussian-chirp-check/SOURCES.json b/research/pr4554-gaussian-chirp-check/SOURCES.json new file mode 100644 index 0000000000..a2a5ea8bbd --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/SOURCES.json @@ -0,0 +1,18 @@ +{ + "schema_version": 1, + "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "ed33caef5c5827693e9572d252a66f99cbb73d40", + "sha256": "5732707cc1e460b090b53f2463162749ad860de878a18f7408938671d295a9d9" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "64bfc5af6709ed9811171e0a2108f05e1219cb65", + "sha256": "352a6184713b8843d96fd2c2071eb7aeaa655906b7c51d45ad444912666c8a96" + } + ] +} diff --git a/research/pr4554-gaussian-chirp-check/TASK.json b/research/pr4554-gaussian-chirp-check/TASK.json new file mode 100644 index 0000000000..fc14318967 --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/TASK.json @@ -0,0 +1,100 @@ +{ + "schema_version": 1, + "task_id": "independent-moment", + "priority_revision": 7, + "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1", + "original_card": { + "id": "independent-moment", + "question": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?", + "starting_points": [ + "docs/research/research-family-status.md", + "certificates/selected-result.json", + "research/staged544-mixed-fourth-entry/README.md", + "research/staged-fourth-mixture-ceiling/README.md", + "research/r437-two-window-carrier/README.md" + ], + "methods": [ + "Independently reconstruct one changed paid ledger or pooled dual with exact source and model binding.", + "Check a genuinely uncovered moment implication with a different enclosure method and semantic controls." + ], + "acceptance_gate": "Demonstrate independence, source binding, positive cases and semantic failure controls. Name the existing checks and the missing coverage. More digits or another equivalent evaluation of an unchanged root is not by itself new verification coverage.", + "first_step": { + "command": [ + "python3", + "-B", + "scripts/research_baseline.py" + ], + "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY", + "inspect": "Choose one uncovered implication in #4088 literal mixed fill, #4091 pooled dual/model binding or #4067 changed frame/host bill. #4091 pins an earlier #4088 head: validate word, capacity and stock equivalence before applying that ceiling to the later fill.", + "scope": "Source-binding orientation only; no scientific replay, theorem acceptance or supplier admission. Budget any heavy replay separately." + }, + "first_experiment": "Independently reconstruct one changed ledger or pooled certificate with a positive case and a semantic failure control. Compare exact sources and name what existing checks miss. Keep this separate from supplier-boundary-audit (ordinary/outer dispatch) and uniform-transfer-audit (all-size proof). Stop after the chosen implication; more digits of an unchanged root are not additional mechanism coverage.", + "bottleneck_link": "Independent pricing prevents unsupported frontier claims and can reveal genuine accounting slack.", + "allowed_variations": [ + "Reconstruct either a source ledger or a numeric enclosure with a different method." + ], + "scope_change_triggers": [ + "Importing the primary function and calling its output independent verification.", + "Calling an aggregate histogram or source archive an independently accepted physical supplier." + ], + "overlap_terms": [ + "independent moment", + "rank histogram", + "exact enclosure", + "2979", + "rank7 joint endframes", + "paired extraction", + "joint lowering", + "4088", + "4091", + "4067", + "staged544-mixed-fourth-entry", + "staged-fourth-mixture-ceiling", + "r437-two-window-carrier" + ], + "existing_checks": [ + { + "scope": "PR #4088: mixed local/owner/phase ledger and moment reported.", + "evidence": "research/staged544-mixed-fourth-entry/README.md", + "remaining": "Independent binding of changed construction to paid histogram, distinct from full admission.", + "status": "author-checked", + "source_commit": "72d7efc181edb21806a878f174847c96dd410ab0", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "PR #4091: sharp pooled ceiling proposed for fixed THREE/FOUR words and inventory.", + "evidence": "research/staged-fourth-mixture-ceiling/README.md", + "remaining": "Validate dual, exhaustive pattern scope and exact source equivalence before applying to newer fills.", + "status": "author-checked", + "source_commit": "a29bdd8f56cd18b805917775cdec09ff935eddc4", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "PR #4067: paid (2,2) to (1,3) exchange and fixed OQ4 host reported.", + "evidence": "research/r437-two-window-carrier/README.md", + "remaining": "Independent full cost delta; remaining caller/precision/stock obligations are separate.", + "status": "author-checked", + "source_commit": "a95f11c52e3ca1e4c5171b3c20d8f4b0815aee2e", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + } + ] + }, + "scope_status": "unchanged", + "actual_scope": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?", + "scope_change_reason": null, + "bottleneck_link": "Independent coverage of the Gaussian block-map construction that #4554 offers as the native replacement for the forward resampling map.", + "budget_minutes": 45, + "stop_rule": "Stop at GAUSSIAN-CHIRP.md; Theorem 3 operation count not re-counted.", + "outcome": "completed", + "evidence": "Exponent identity as a polynomial identity; linear-convolution identity at 60 digits (40 rectangles); 300 exact conditioning cases; 3,522 retained blocks for K, |c0| and (D)=>(A); 116 tail cases; 576 precision-rule cases; own fixed-point block evaluator at the (C) precision: 6 runs with error far below 2^-p0 and within (B); 7 statements bound.", + "author": "Ryanson Jonathan, with Claude Code (Claude Opus 5.5) assistance", + "ai_assistance": "Substantial Claude Code (Claude Opus 5.5) assistance; no independent human review.", + "open_questions": [ + { + "question": "Is a sharper precision rule than (C) valid?", + "why_it_matters": "(C) reserves ~2.7x the observed precision loss.", + "first_step": "Redo (B) tracking M powers per transform.", + "blocked_by": "None." + } + ] +} diff --git a/research/pr4554-gaussian-chirp-check/VERIFY.json b/research/pr4554-gaussian-chirp-check/VERIFY.json new file mode 100644 index 0000000000..650ea02264 --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/VERIFY.json @@ -0,0 +1,18 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [], + "commands": [ + { + "argv": [ + "{python}", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 300, + "checks": "Diagonal-Hankel identity (polynomial + numeric linear convolution), exact conditioning, blocking geometry brute force, tail and precision rules, own fixed-point block evaluator at the (C) precision vs 120-digit direct sums; binding to vendored spec (sha256)." + } + ] +} diff --git a/research/pr4554-gaussian-chirp-check/verify.py b/research/pr4554-gaussian-chirp-check/verify.py new file mode 100644 index 0000000000..8ab494bb03 --- /dev/null +++ b/research/pr4554-gaussian-chirp-check/verify.py @@ -0,0 +1,322 @@ +#!/usr/bin/env python3 +"""Independent check of PR #4554 (Dugongue) GAUSSIAN-CHIRP.md, written from the vendored text only. + +Nothing is imported from #4554's chirp_check.py, gaussian_block.py or native_convolution.py. + 1. Theorem 1: log D[x] + log E[y] + log H[x+y] = -a(ux - vy + c)^2 as a polynomial identity in Q[a,u,v,c,x,y], and + sum_y W[x,y] Z[y] = D[x] (H * reverse(E Z))[x+l-1] with an explicit linear (not cyclic) convolution, numerically + at 60 digits on random rectangles. + 2. Section 2: under (A) with C p := a(|u|q + |v|q + |c0|)^2, every exponent of D (x < q), E (y < q) and + H (z < 2q - 1) has absolute value at most C p, exactly in rationals. + 3. Section 4 geometry, by brute force over random rational slopes, offsets, radii and block sizes: the retained + column blocks of each row block number at most K = ceil((u(q-1)+2R)/(vq)) + 3, every retained block has + |c0| <= R + (u+v)(q-1), and (D) implies (A). The tail rule a R^2 >= (p0 + ceil(log2(sX+1)) + 2) ln 2 makes the + omitted row mass below 2^-p0/2. + 4. The precision rule (C) implies (B) < 2^-p0 (integer comparisons, using only log2 e < 2). + 5. End to end: my own fixed-point block evaluator at the precision P prescribed by (C) (chirps and twiddles on the + 2^-P grid, nearest-rounded complex products, linear convolution via a radix-2 FFT of size N0 = 4q, division by + N0, selection at x+l-1, multiplication by D) against a direct 60-digit sum: error < 2^-p0 and within (B). +Ryanson Jonathan with Claude Code (Claude Opus 5.5) assistance. Apache-2.0. +""" +import decimal +from decimal import Decimal as Dec +from fractions import Fraction as Fr +import hashlib +import json +import math +from pathlib import Path +import random +import sys + +HERE = Path(__file__).resolve().parent +SPEC = HERE / "PR4554-GAUSSIAN-CHIRP.md" # #4554 head 4f46c6a9, research/native-word-fft/GAUSSIAN-CHIRP.md, blob 4e31f302 +SPEC_SHA256 = "d8f2c45a985fc51d9523c3261ec39fc622afe737bf26a868e7f601dc7525b075" +VARS = ("a", "u", "v", "c", "x", "y") + + +def need(condition, message): + if not condition: + raise AssertionError(message) + + +class Poly: + def __init__(self, t=None): + self.t = {k: v for k, v in (t or {}).items() if v} + + @staticmethod + def var(n): + k = [0] * len(VARS) + k[VARS.index(n)] = 1 + return Poly({tuple(k): Fr(1)}) + + @staticmethod + def const(c): + return Poly({(0,) * len(VARS): Fr(c)}) + + def __add__(self, o): + o = o if isinstance(o, Poly) else Poly.const(o) + t = dict(self.t) + for k, v in o.t.items(): + t[k] = t.get(k, 0) + v + return Poly(t) + + __radd__ = __add__ + + def __neg__(self): + return Poly({k: -v for k, v in self.t.items()}) + + def __sub__(self, o): + return self + (-(o if isinstance(o, Poly) else Poly.const(o))) + + def __mul__(self, o): + o = o if isinstance(o, Poly) else Poly.const(o) + t = {} + for k1, v1 in self.t.items(): + for k2, v2 in o.t.items(): + k = tuple(i + j for i, j in zip(k1, k2)) + t[k] = t.get(k, 0) + v1 * v2 + return Poly(t) + + __rmul__ = __mul__ + + +def check_identity(): + a, u, v, c, x, y = (Poly.var(n) for n in VARS) + logD = -a * (u * u + u * v) * x * x - 2 * a * u * c * x - a * c * c + logE = -a * (v * v + u * v) * y * y + 2 * a * v * c * y + z = x + y + logH = a * u * v * z * z + target = -a * (u * x - v * y + c) * (u * x - v * y + c) + need(not (logD + logE + logH - target).t, "exponent identity") + # Numeric check with an explicit linear convolution. + rng = random.Random(1) + decimal.getcontext().prec = 70 + worst = Dec(0) + for _ in range(40): + A, U, V, Cc = (Dec(rng.uniform(0.01, 0.4)), Dec(rng.uniform(0.3, 2)), Dec(rng.uniform(0.3, 2)), + Dec(rng.uniform(-2, 2))) + m, l = rng.randint(1, 7), rng.randint(1, 7) + Z = [(Dec(rng.uniform(-1, 1)), Dec(rng.uniform(-1, 1))) for _ in range(l)] + D = [(-A * (U * U + U * V) * xx * xx - 2 * A * U * Cc * xx - A * Cc * Cc).exp() for xx in range(m)] + E = [(-A * (V * V + U * V) * yy * yy + 2 * A * V * Cc * yy).exp() for yy in range(l)] + H = [(A * U * V * zz * zz).exp() for zz in range(m + l - 1)] + B = [(E[yy] * Z[yy][0], E[yy] * Z[yy][1]) for yy in range(l)] + rev = B[::-1] + conv = [(Dec(0), Dec(0)) for _ in range(len(H) + l - 1)] + for i, h in enumerate(H): + for j, bb in enumerate(rev): + conv[i + j] = (conv[i + j][0] + h * bb[0], conv[i + j][1] + h * bb[1]) + for xx in range(m): + lhs = [Dec(0), Dec(0)] + for yy in range(l): + w = (-A * (U * xx - V * yy + Cc) ** 2).exp() + lhs[0] += w * Z[yy][0] + lhs[1] += w * Z[yy][1] + rhs = (D[xx] * conv[xx + l - 1][0], D[xx] * conv[xx + l - 1][1]) + worst = max(worst, abs(lhs[0] - rhs[0]) + abs(lhs[1] - rhs[1])) + need(worst < Dec(10) ** -50, "numeric identity with linear convolution at index x+l-1") + return 41 + + +def check_conditioning(rng): + cases = 0 + for _ in range(300): + q = rng.randint(1, 24) + a = Fr(rng.randint(1, 50), rng.randint(50, 400)) + u = Fr(rng.randint(-40, 40), rng.randint(1, 20)) + v = Fr(rng.randint(-40, 40), rng.randint(1, 20)) + c0 = Fr(rng.randint(-60, 60), rng.randint(1, 20)) + Cp = a * (abs(u) * q + abs(v) * q + abs(c0)) ** 2 + for xx in range(q): + need(abs(a * (u * u + u * v) * xx * xx + 2 * a * u * c0 * xx + a * c0 * c0) <= Cp, "D exponent") + need(abs(-a * (v * v + u * v) * xx * xx + 2 * a * v * c0 * xx) <= Cp, "E exponent") + for zz in range(2 * q - 1): + need(abs(a * u * v * zz * zz) <= Cp, "H exponent") + cases += 1 + return cases + + +def check_geometry(rng): + blocks = 0 + for _ in range(120): + u = Fr(rng.randint(5, 40), 10) + v = Fr(rng.randint(5, 40), 10) + c = Fr(rng.randint(-200, 200), 10) + R = Fr(rng.randint(1, 80), 10) + q = 1 << rng.randint(0, 4) + t, s = rng.randint(q, 96), rng.randint(q, 96) + K = math.ceil((u * (q - 1) + 2 * R) / (v * q)) + 3 + for x0 in range(0, t, q): + kept = 0 + for y0 in range(0, s, q): + keep = any(abs(u * xx - v * yy + c) <= R for xx in range(x0, min(t, x0 + q)) + for yy in range(y0, min(s, y0 + q))) + if keep: + kept += 1 + c0 = u * x0 - v * y0 + c + need(abs(c0) <= R + (u + v) * (q - 1), "|c0| <= R + (u+v)(q-1)") + a = Fr(1, 7) + need(a * (R + 2 * (u + v) * q) ** 2 >= a * (u * q + v * q + abs(c0)) ** 2, "(D) implies (A)") + blocks += 1 + need(kept <= K, "at most K retained column blocks") + tails = 0 + for p0 in range(1, 200, 7): + for sX in (1, 7, 1000, 10 ** 9): + k = p0 + math.ceil(math.log2(sX + 1)) + 2 + need(Fr(sX, 2 ** k) < Fr(1, 2 ** (p0 + 1)), "tail below 2^-p0/2") # exp(-aR^2) = 2^-k + tails += 1 + return blocks, tails + + +def precision_rule(p0, q, C, p, X): + l2 = lambda v: math.ceil(math.log2(v)) + return p0 + 30 + 8 * l2(q + 2) + 8 * math.ceil(C) * p + 2 * l2(X + 1) + 4 + + +def check_precision_rule(): + cases = 0 + for p0 in (1, 20, 64, 200): + for q in (1, 2, 16, 1024): + for C in (Fr(1, 2), 1, 3): + for p in (1, 10, 50): + for X in (0, 1, 100, 10 ** 6): + P = precision_rule(p0, q, C, p, X) + # 2^30 q^8 M^4 (X+1)^2 2^-P <= 2^(-p0-4), with M^4 = e^(4Cp) < 2^(8Cp) as log2 e < 2. + lhs_log2_upper = 30 + 8 * math.log2(q) + 8 * C * p + 2 * math.log2(X + 1) - P + need(lhs_log2_upper <= -p0 - 4, "(C) implies (B) <= 2^(-p0-4)") + cases += 1 + need(math.e < 4, "log2 e < 2") + return cases + + +# ---------------------------------------------------------------- end-to-end fixed-point block +def end_to_end(rng): + results = [] + decimal.getcontext().prec = 120 + for a, u, v, c, q, p0 in ((Fr(1, 8), Fr(1), Fr(3, 2), Fr(1, 3), 8, 20), + (Fr(1, 20), Fr(5, 4), Fr(1, 2), Fr(-2), 16, 24), + (Fr(1, 4), Fr(1, 2), Fr(1, 2), Fr(0), 8, 30)): + m = l = q + X = 1 + C = 1 + p = math.ceil(a * (abs(u) * q + abs(v) * q + abs(c)) ** 2) # (A) with C = 1 + P = precision_rule(p0, q, C, p, X) + 8 + S = 1 << P + A, U, V, Cc = (Dec(f.numerator) / Dec(f.denominator) for f in (a, u, v, c)) + + def grid(d): + return int((d * S).to_integral_value(rounding=decimal.ROUND_HALF_EVEN)) + + D = [(-A * (U * U + U * V) * x * x - 2 * A * U * Cc * x - A * Cc * Cc).exp() for x in range(m)] + E = [(-A * (V * V + U * V) * y * y + 2 * A * V * Cc * y).exp() for y in range(l)] + H = [(A * U * V * z * z).exp() for z in range(m + l - 1)] + Dg, Eg, Hg = [grid(d) for d in D], [grid(e) for e in E], [grid(h) for h in H] + N0 = 4 * q + pi = Dec(0) + # 120-digit pi by Machin, for twiddles on the grid (certified-interface error <= 2^-P). + def arctan_inv(k): + total, term, j, x2 = Dec(0), Dec(1) / k, 0, k * k + while abs(term) > Dec(10) ** -115: + total += term / (2 * j + 1) * (-1 if j % 2 else 1) + term /= x2 + j += 1 + return total + pi = 16 * arctan_inv(5) - 4 * arctan_inv(239) + + def cexp(theta): + cs, sn, term, k = Dec(0), Dec(0), Dec(1), 0 + while abs(term) > Dec(10) ** -115 or k < 4: + if k % 4 == 0: + cs += term + elif k % 4 == 1: + sn += term + elif k % 4 == 2: + cs -= term + else: + sn -= term + k += 1 + term = term * theta / k + return cs, sn + + tw = [tuple(grid(w) for w in cexp(2 * pi * k / N0)) for k in range(N0 // 2)] + + def rmul(x, y): + re, im = x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0] + return ((re + (S >> 1)) >> P, (im + (S >> 1)) >> P) + + def fft(arr, inverse): + arr = list(arr) + n = len(arr) + j = 0 + for i in range(1, n): + bit = n >> 1 + while j & bit: + j ^= bit + bit >>= 1 + j |= bit + if i < j: + arr[i], arr[j] = arr[j], arr[i] + size = 2 + while size <= n: + step = n // size + for st in range(0, n, size): + for k in range(size // 2): + w = tw[k * step] + if not inverse: + w = (w[0], -w[1]) + tt = rmul(w, arr[st + k + size // 2]) + uu = arr[st + k] + arr[st + k] = (uu[0] + tt[0], uu[1] + tt[1]) + arr[st + k + size // 2] = (uu[0] - tt[0], uu[1] - tt[1]) + size *= 2 + return arr + + Hhat = fft([(h, 0) for h in Hg] + [(0, 0)] * (N0 - len(Hg)), False) + for trial in range(2): + Z = [(rng.randint(-S, S), rng.randint(-S, S)) for _ in range(l)] + Zf = [(Fr(zr, S), Fr(zi, S)) for zr, zi in Z] + need(all(abs(zr) <= 1 and abs(zi) <= 1 for zr, zi in Zf), "|Z| bounded") + B = [rmul((Eg[y], 0), Z[y]) for y in range(l)] + Bhat = fft(B[::-1] + [(0, 0)] * (N0 - l), False) + prod = [rmul(x, y) for x, y in zip(Hhat, Bhat)] + conv = fft(prod, True) + conv = [((cr + N0 // 2) // N0, (ci + N0 // 2) // N0) for cr, ci in conv] + worst = Dec(0) + for x in range(m): + out = rmul((Dg[x], 0), conv[x + l - 1]) + ref = [Dec(0), Dec(0)] + for y in range(l): + w = (-A * (U * x - V * y + Cc) ** 2).exp() + ref[0] += w * Dec(Z[y][0]) / S + ref[1] += w * Dec(Z[y][1]) / S + err = ((Dec(out[0]) / S - ref[0]) ** 2 + (Dec(out[1]) / S - ref[1]) ** 2).sqrt() + worst = max(worst, err) + M = Dec(C * p).exp() + bound_B = Dec(2) ** 30 * q ** 8 * M ** 4 * (X + 1) ** 2 / Dec(2) ** P + need(worst < Dec(2) ** -p0 and worst <= bound_B, "end-to-end error < 2^-p0 and within (B)") + results.append({"q": q, "p": p, "P": P, "p0": p0, "log2_error": float((worst + Dec(10) ** -100).ln() / Dec(2).ln())}) + return results + + +def check_binding(): + raw = SPEC.read_bytes() + need(hashlib.sha256(raw).hexdigest() == SPEC_SHA256, "vendored GAUSSIAN-CHIRP.md bytes") + text = raw.decode("utf-8") + for s in ("sum_y W[x,y]*Z[y] = D[x] * (H * reverse(E*Z))[x+l-1],", "a*(|u|*q+|v|*q+|c0|)^2 <= C*p.", + "block error <= 2^30*q^8*M^4*(X+1)^2 * 2^-P,", "K = ceil((u*(q-1)+2R)/(v*q)) + 3", + "|c0| <= R+(u+v)*(q-1),", "a*(R+2*(u+v)*q)^2 <= C*p.", "a*R^2 >= (p0+ceil(log2(s*X+1))+2)*log(2)"): + need(s in text, f"spec states: {s}") + return 7 + + +def main(): + rng = random.Random(45543) + blocks, tails = check_geometry(rng) + report = {"status": "PASS_PR4554_GAUSSIAN_CHIRP_CHECK", "identity_checks": check_identity(), + "conditioning_cases": check_conditioning(rng), "retained_blocks_checked": blocks, "tail_cases": tails, + "precision_rule_cases": check_precision_rule(), "end_to_end": end_to_end(rng), + "spec_statements_bound": check_binding(), "new_kappa_claimed": False} + print(json.dumps(report, indent=1)) + + +if __name__ == "__main__": + sys.exit(main())