diff --git a/.github/workflows/h3-band-v2.yml b/.github/workflows/h3-band-v2.yml new file mode 100644 index 000000000..d00c154e3 --- /dev/null +++ b/.github/workflows/h3-band-v2.yml @@ -0,0 +1,40 @@ +name: h3-band-v2 +on: + pull_request: + branches: [main] + paths: + - 'experiments/h3_band_v2/**' + - '.github/workflows/h3-band-v2.yml' + push: + branches: [main] + paths: + - 'experiments/h3_band_v2/**' + - '.github/workflows/h3-band-v2.yml' + workflow_dispatch: +permissions: + contents: read +concurrency: + group: h3-band-v2-${{ github.ref }} + cancel-in-progress: true +jobs: + replay: + name: h3-band-v2 + runs-on: ubuntu-latest + timeout-minutes: 40 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 + with: + persist-credentials: false + - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 + with: + python-version: '3.11' + - name: Install the pinned mpmath the reviewed transcripts were made with + run: | + python -m pip install --require-hashes --only-binary=:all: -r experiments/h3_band_v2/requirements.txt + - name: Structure lemma, records, replays and every mutation control + run: python -B -m unittest discover -s experiments/h3_band_v2 -p 'test_*.py' -v + - name: The same suite under -O, replays and negative controls included + run: python -B -O -m unittest discover -s experiments/h3_band_v2 -p 'test_*.py' -v + - name: State scientific scope + run: | + echo 'Replayed a draft H3 band certificate successor that is under review. A green run is a run, not a review, and it changes no claim status. Scientific effect NONE.' >> "$GITHUB_STEP_SUMMARY" diff --git a/experiments/h3_band_v2/README.md b/experiments/h3_band_v2/README.md new file mode 100644 index 000000000..1f172c740 --- /dev/null +++ b/experiments/h3_band_v2/README.md @@ -0,0 +1,156 @@ +# H3 band certificate, successor v2 + +**Status: a draft for review. Scientific effect NONE. Review status REVIEW_REQUIRED.** + +The author is Anthropic Claude (session `session_01Fij9F1cYrAUDW7PhWQze1G`). The package needs a reviewer from another model family. A reader from the same provider earns zero independence credit. For that reviewer, [`REVIEW_GUIDE.md`](REVIEW_GUIDE.md) states and derives, with exact line references, the four claims that rest on a written derivation rather than on a run. Nothing here changes any claim, status, HOLD or consumer pin. Those decisions belong to the H3/SIDE24 owners and the consumers. + +This is the numbered successor of two frozen programs, both authored by Kimi: +- `h3_band_ceil.py` (sha256 `b97c5428…`); +- `h3_band_floor.py` (sha256 `a907eeed…`). + +Both live on main under `incoming/claude-bundle-v3-native-copies-20260928/K3_SIDE24_LB/UPPER2D/H3_closure/` as `.py.txt`. They bound E[G_r] for r in (0, 0.05]. G_r is the normalizer of the LPW_CONSTANT v4 / OBL-D1-PROMOTE band, defined in the frozen `H3_BAND_CEIL.md` and `H3_BAND_FLOOR.md`. In brief: +- (q_M, q_S, a_M, a_S, b_M, b_S) is the conditional Gaussian law of the six soft variables, given the corrected pins of the exact side-24 periodized field (`h2/cov_exact.py`); +- D_M = a_M q_M − r b_M² and D_S = a_S q_S − r b_S²; +- G_r = −D_M D_S on the event {q_M < 0, D_M > 0, D_S < 0}. + +This package does not re-derive that reduction. It certifies the expectation as the frozen programs define it. + +The review that motivated this successor is O7 on [main #419](https://github.com/d6g8k5htny-coder/main/issues/419#issuecomment-6245249123). It found that the frozen programs do not certify cell 1, (0, 0.0025], or the published uniform constants. + +## What the replay certifies, subject to review + +Each row is certified for every r in the cell. The decimals are rounded outward exactly: the ceiling upward and the floor downward. They come from `results/ceil.json` and `results/floor.json`. + +| cell | law method | ceiling E[G_r] ≤ | floor E[G_r] ≥ | +|---|---|---|---| +| (0, 1/16384] | R-majorant | 3.65718683240099 | 2.61755715030352 | +| (1/16384, 1/400] | Neumann | 3.66271691969455 | 2.58373305711774 | +| (1/400, 1/200] | Neumann | 3.66435028046303 | 2.58174158877405 | +| (1/200, 1/100] | Neumann | 3.67360595202420 | 2.54243388604045 | +| (1/100, 1/50] | Neumann | 3.69805856644731 | 2.45593199230241 | +| (1/50, 3/100] | Neumann | 3.71138266685765 | 2.44021698557766 | +| (3/100, 1/25] | Neumann | 3.72414747013089 | 2.41538548140069 | +| (1/25, 1/20] | Neumann | 3.73661723530033 | 2.28610115562878 | +| **uniform (0, 1/20]** | | **3.73661723530033** | **2.28610115562878** | + +Compared with the frozen programs: +- **The ceiling over (0, 0.0025]** is now at most 3.66271691969455. That is below the frozen cell-1 value 3.66282864761194, which LPW v4 consumes. LPW v4's r₀ = 1/2278031360 lies in the first cell, whose ceiling is 3.65718683240099. +- **The uniform ceiling** is 3.73661723530033. The frozen value was 3.74767948915996. The rigorous remainder (B2) turns out tighter than the frozen correction, so the ceiling moves down. +- **The uniform floor** is 2.28610115562878. The frozen value was 2.30659559567154. The two-tail fix (B3) raises the δ piece in the last cell to 0.121, so the floor moves down. It stays far above c_Z = 1.6154892676…. + +Whether any consumer re-pins to these values is the consumer's decision. + +## Changes from the frozen programs + +The law pipeline, CB0–CB6, is kept stage by stage, and the two envelopes share it in one program. The upper-envelope algebra (`ceil_hi`) and the lower-envelope algebra (`band_lo`) are unchanged, except for the items below and two checks. `band_lo` no longer carries the frozen floor's check that Edelta and Ebeta exceed −1e−6, and both envelopes now require their second moments to be nonnegative before a square root is taken. Neither affects validity when the enclosures are valid. + +- **B1. Majorants.** `matmax` and `vecmax` now take max(|a|, |b|) of each interval. The frozen versions took the signed upper endpoint. A new guard checks that every majorant dominates every entry it majorizes. For the unmutated code it holds by construction: it is a regression guard, and it is what catches `signed_maxima`. +- **First cell.** True majorants give a majorant growth of R ≈ 2622. So the R-majorant cell is narrowed to (0, 1/16384], where R·r* ≈ 0.160 < 0.2, and the Neumann method covers the rest of (0, 0.0025]. The 0.2 and the condition θ(R) ≤ 1/2 are margins kept from the frozen program: the induction itself needs θ(R) ≤ 1, and convergence needs R·r* < 1 and r*/ρ < 1, which the tail bounds check. The frozen first cell, (0, 0.0025], fails the guard and the majorant series genuinely diverges there (R/400 ≈ 6.6); this is the `wide_first_cell` control. On the first cell the pole coefficients of μ and Σ are dropped on CB2s alone, with no numeric fallback, because 1/t is unbounded there; their vanishing follows from CB2s (c) and (d) by direct algebra. +- **B2. Neumann remainder.** The correction now comes from the exact identities + + G (w − W_P) = (zV − G_P W_P) − (G − G_P) W_P + G (G⁻¹ − H_P) = (I − G_P H_P) − (G − G_P) H_P + + together with the exact expansions of μ = A w and Σ = GY − A G⁻¹ Aᵀ. Here: + - W_P and H_P are the truncated Taylor series of G⁻¹zV and G⁻¹, and G_P is the truncated pin Gram; + - the residual polynomials zV − G_P W_P and I − G_P H_P are computed in full, through order 2(J+2); + - so are the high orders of A_P W_P and of A_P H_P A_Pᵀ, through orders 2(J+2) and 3(J+2); + - every Gram block has a Cauchy tail from order J+3; + - ‖G(t)⁻¹‖ is bounded on each cell by a midpoint-preconditioned Neumann contraction. + + The frozen formula is kept only as the `old_remainder` control. +- **B3.** The floor bounds the two-sided event {|δ| > η} with two Gaussian tails. Every Gaussian tail in either envelope goes through one function, `event_prob_upper`, and the number of tails follows from the event kind. At start-up that function and the tail helper are self-tested on exact point laws, for N(0, 1) and N(1/2, 1), against Gordon's lower bound φ(t) t / (t² + 1). That bound does not use `iv_Phi`. +- **B4. Negative controls.** There are fourteen mutations, each invoked by `--mutation NAME` and each caught: twelve inside the run and two by the replay; see below. A mutated run never reports a pass: when it is not caught inside the run, it ends `H3_BAND_V2 MUTATED RUN …, NOT A CERTIFICATE` and exits 3. The same holds for `structure`: under a mutation the lemma cannot see, it ends `H3_BAND_V2 MUTATED STRUCTURE RUN …, NOT A CERTIFICATE` and exits 3. `--record` is refused under a mutation, and `--record` cannot be combined with `--verify`. +- **Fail-closed guards.** Four conditions that no current cell reaches now stop the run at a named check rather than at an mpmath assertion or a traceback: r < ρ in the Gram tail (`gram_tail`), a positive variance before each square root (`t0_lower`, `trunc_mom_iv`), and a window symmetric about 0 for the two-sided tail in `band_lo`. None changes a certified decimal, and the reviewed transcripts are byte-identical. +- **O2.** The Gordon–Mills branch of `iv_Phi` keeps its endpoints as mpf, not float. +- **CB2s, a new exact lemma, in `fractions.Fraction` only.** It covers the exact Laurent expansions of the pin and soft functionals: + - (a) the pins have no negative powers; + - (b) the soft functionals have no power below z⁻¹; + - (c) Y⁽⁻¹⁾ = E V⁽⁰⁾ for an exact rational matrix E (for the tables of record E is 0/1); + - (d) E zV = 0. + + Hence Ỹ = Y − EV/z is regular. That makes μ = E[Ỹ | V = zV] and Σ = Cov(Ỹ | V) analytic wherever G is invertible: the pole coefficients of μ and Σ vanish exactly. The frozen programs dropped those coefficients after a numeric |c| < 1e−55 check. The numeric checks stay, as consistency checks only: they see orders from −8 up, with a 1e−60 threshold, which could not by itself protect the remainder. Exact vanishing rests on CB2s. A₋₁ = E G₀ is checked at the Gram level. The lemma and the interval series are built from the same table objects, so a change to a table reaches both, and each of the four branches has its own control. +- **Kernel window.** The series now uses `range(0, J+9−n)`; the frozen one used `J+8−n`. Every Gram block is complete through order J+2, and the Cauchy tails start at J+3. The soft-pin block is kept through J+2 as well. +- **Moments.** The Poisson sums over |n| ≤ 3 now carry a certified tail for |n| ≥ 4, built from |He_k(x)| ≤ (|x| + √k)^k. New check CB1d: every consumed a₂ₘ, m ≤ 40, must meet an independent spectral enclosure of a₂ₘ = (1/Z) Σ e^{−k²/2} k^{2m}. It sums the lattice |j| ≤ 116, adds a certified tail to the numerator and divides by the truncated lattice sum Z_T. That still encloses a₂ₘ: every omitted |k| exceeds every retained one, so N_T/Z_T ≤ a₂ₘ ≤ (N_T + tail)/Z_T (scoped CB1 review 5607583647, §5). CB1d is a consistency check, since meeting an enclosure is not inclusion. The kernel sup A_n(ρ) carries a certified lattice tail too, and it uses |k| upper endpoints. +- **Rounding.** The load-bearing upper endpoints come from intervals. That covers: + - CB4 (θ(R), C_W, C_H, the Cauchy ratio 1/ρ, which is rounded upward, and every tail sum, which is now a sum of nonnegative terms rather than a difference). A sanity check requires every computed Laurent coefficient, orders −2 to J+2, to obey its Cauchy bound |c_j| ≤ M ρ^{−j}. The Cauchy bounds the tails use, at higher orders, rest on the derivation; + - the Neumann bounds; + - the Chernoff t₀ lower bounds and the Cauchy–Schwarz pieces; + - the type and cover margins, checked exactly in `Fraction` from the mpf box constants. + + Cell endpoints are enclosed outward, so no gap opens between cells. +- **H2 library.** CB0 executes exactly the bytes whose digest it checked; no import machinery or bytecode cache is involved. +- **Output.** Certified decimals are rounded outward, exactly. The transcript has no timings: they go to stderr. The H2 path is a CLI argument; the frozen programs hard-coded `/mnt/agents/...`. Monte Carlo is opt-in (`--mc`), NON-CERTIFYING, and goes to stderr. numpy is needed only for it. +- **Removed.** The frozen ceiling checked nesting against one global floor value. That check is now a per-cell test across the two records. + +## What this does not establish + +- **Trust in mpmath.** mpmath's `iv` context is trusted as a correct outward-rounding interval arithmetic, as in the frozen programs. This repository's preferred enclosure library, `research/interval/`, is not present on main. +- **The envelope decomposition.** The algebra (type event, covers, p̄ ≥ |p|, the box `B_cert`) is unchanged and was checked in O7. So far it has had only an Anthropic reading. +- **Two mutations are invisible inside a run.** `drop_fail_piece` and `old_remainder` give plausible numbers. They are caught only because a replay must reproduce the reviewed record exactly: the program digest, the transcript digest, and every certified decimal. +- **The upstream reduction.** This package does not re-derive the definition of G_r or the reduction to it. It does not re-pin LPW v4, H5_ZBAND or D1. +- **The B2 remainder has no numeric test.** CB7, the direct point law at r = 1/20 in mp arithmetic, checks the heavy polynomial, not the remainder. Its printed deviation, about 2e−91, is the rounding error of its own dps-100 reference law; against a dps-200 law the Neumann midpoint at r = 1/20 is within about 1e−100. Either figure is far below the rigorous corrections 1.5e−10 and 1.2e−7, so CB7 cannot tell whether the remainder is under-estimated. The remainder rests on the written derivation in `neumann_law`, set out in `REVIEW_GUIDE.md`. A same-family check, with zero independence credit, compared every cell enclosure against a dps-200 direct law at 32 points, and every enclosure contained it. +- **The B3 control covers the shared event bound.** It does not cover a future call site that bypasses `event_prob_upper`. +- **Planar moments.** For m ≤ 4, the torus and planar moments differ by about 1e−122, which is below dps 100. CB1d separates them from m = 16 up. +- **Non-certifying parts.** CB7 is a cross-check and not load-bearing. So is the Monte Carlo. + +## Run + +``` +python3 -m pip install --require-hashes --only-binary=:all: -r experiments/h3_band_v2/requirements.txt +python3 -B experiments/h3_band_v2/h3_band_v2.py structure +python3 -B experiments/h3_band_v2/h3_band_v2.py ceil --verify experiments/h3_band_v2/results/ceil.json +python3 -B experiments/h3_band_v2/h3_band_v2.py floor --verify experiments/h3_band_v2/results/floor.json +python3 -B -m unittest discover -s experiments/h3_band_v2 -p 'test_*.py' -v +``` + +- **Runtime:** each envelope takes about 20 s on CPython 3.11 with mpmath 1.4.1. +- **Byte-exactness:** the reviewed transcripts are byte-exact under that mpmath, and the same under `-O`. CI runs the whole suite, negative controls included, in both modes; under `-O` every subprocess the tests start runs with `-O` too. +- **Output:** a run prints one summary line, `H3_BAND_V2 PASS …`, and exits 0. It exits: + - 1 on any failed check (`H3_BAND_V2 FAIL: …`) or any record mismatch (`H3_BAND_V2 VERIFY_FAIL: …`); + - 2 on a refused combination of options; + - 3 for a mutated run that completed, which is never a certificate. +- **Comparing a stored transcript:** `verify-record --transcript F --results R` checks a transcript against a record without recomputing. It needs only the standard library, and it does not compare the verifier's mpmath version with the record's; only a recomputing `--verify` replay does. It does require the record to name the pinned producer mpmath (`MPMATH_PIN`, the same version as `requirements.txt`), and it reads the mode from the transcript's first line rather than from the record. + +## Negative controls + +`--mutation NAME`. `test_h3_band_v2.py` runs every one and requires each to fail at the named check. + +| mutation | caught by | +|---|---| +| `structure_pin_pole` | CB2s (a): pin V5 has a z⁻² term | +| `structure_soft_pole2` | CB2s (b): soft Y2 has a z⁻² pole | +| `structure_soft_span` | CB2s (c): the pole part of Y0 is outside the span of V⁽⁰⁾ | +| `structure_window` | CB2s (d): E zV ≠ 0 | +| `planar_moments` (planar values in the consumed table A2) | CB1d: spectral cross-check | +| `wrong_window` | CB3: limit identity μ_aM(0) = −1 | +| `signed_maxima` (B1 regression) | CB4: majorant dominance guard | +| `wide_first_cell` (frozen first cell) | CB4: R-majorant convergence on (0, 1/400] | +| `one_tail` (B3 regression in the shared event bound) | CB5: event-bound self-test, P(|X| > t) | +| `tail_helper_one_side` | CB5: tail-helper self-test | +| `no_preconditioner` | CB6: Neumann contraction | +| `cap_raised` | CC6: consumption cap U ≤ 4 | +| `drop_fail_piece` | replay against the reviewed record; the certified decimals change | +| `old_remainder` (B2 regression) | replay against the reviewed record; the certified decimals change | + +The record checks are tested too. Each of the following is rejected: +- a lowered bound in a transcript; +- a doctored record; +- a record bound to other program bytes; +- a mutated transcript, or one without `ALL_CHECKS_PASS`, even with a record made consistent with it; +- a transcript and record doctored consistently with each other. `verify-record` alone accepts that pair; only a recomputing replay exposes it; +- a transcript whose only change is a check line that carries no certified decimal. Only the transcript digest binds those lines; +- a record that names no producer mpmath or one other than the pin, or whose environment block is missing or not an object, by `verify-record` and by a recomputing replay. `verify-record`'s outcome does not depend on the verifier's mpmath version; only a recomputing replay requires the verifier to run the pinned version; +- a record relabelled to the other envelope, because the mode is read from the transcript's first line; +- a transcript whose first line names no envelope, even with a record made consistent with it; +- a record that is valid JSON but not an object, of any JSON type, by `verify-record` and by a replay before it runs. + +## Provenance + +- **`h2/cov_exact.py`:** 21,961 bytes, sha256 `f08c1c5f653f2ffd2e85d39e1112f80d15db04d15f932e5ee42db2ca3553e783`. It is byte-identical to the H2 library in the frozen archive `a2136bc0…` (DG-MIGRATION custody). Its digest is the one at line 404 of main's `incoming/claude-bundle-v3-native-copies-20260928/K3_SIDE24_LB/RETURN_06/TREE_MANIFEST.sha256.txt`. The program refuses any other bytes. + - **Known library defect, not consumed here.** The library's interval evaluators (`e_sum`, `s_sum`, `k1_image`, and `k1`, `cov` and `moment_a` in `mode='iv'`) divide by the truncated spectral sum Z_T and pad only the numerator. Their error term omits |r_T|·(Z − Z_T)/Z, so at non-default cutoffs or higher precision their intervals can exclude the true value (review 5607640770). Even-order spectral moments at 0, as in CB1d, are still enclosed, by the argument above. + - The program never calls those evaluators for a certified value. A2 and DEV2/DEV4 come from Poisson sums with their own denominator tails. + - The program's only load-bearing use of Z_T is the CB4 majorant A_n(ρ), which divides a nonnegative sum by an interval whose lower end is at most Z_T, which is at most Z. That can only make the bound larger. + - CB1d (above) and CB7 are checks only. The library itself is frozen and stays byte-identical. A repaired library would be a separate successor. +- **Records:** `results/*.json` bind the program digest, the cov_exact digest, the producer's mpmath version, which must equal the pin `MPMATH_PIN` (a recomputing replay also requires the verifier to run it), the transcript digest and every certified decimal. They also require an unmutated run that passed every check. `results/*.transcript.txt` are the matching transcripts. +- **Re-recording:** any edit to `h3_band_v2.py` requires regenerating the records with `--record`, and that regeneration is itself a change for review. diff --git a/experiments/h3_band_v2/REVIEW_GUIDE.md b/experiments/h3_band_v2/REVIEW_GUIDE.md new file mode 100644 index 000000000..a913e450f --- /dev/null +++ b/experiments/h3_band_v2/REVIEW_GUIDE.md @@ -0,0 +1,769 @@ +# Review guide: H3 band certificate successor v2 + +**Status: a draft for review. Scientific effect NONE. Review status REVIEW_REQUIRED.** This guide is for the package `experiments/h3_band_v2/` of PR #447. It was written by the Anthropic author's session (`session_01Fij9F1cYrAUDW7PhWQze1G`), and same-provider agents checked each section. It therefore carries zero independence credit. It is a map for a reviewer from another model family, not a review. For each claim that rests on a written derivation rather than on a run, it gives the statement, the derivation with line references, what the run checks, and what it does not check. Where the guide and the code disagree, the code governs, and the disagreement is a finding. + +**References.** +- `h3_band_v2.py:N`, or a bare line number inside a scope section, means line N of `experiments/h3_band_v2/h3_band_v2.py`. That file has sha256 `7776d78a0df5ff04d6aabd9f0e853c1d70c69ed039007919e77b64c182860431`, the digest bound in `results/ceil.json` and `results/floor.json`. +- `README:N` is the package README.md, `test:N` is `test_h3_band_v2.py`, and `cov_exact.py:N` is `h2/cov_exact.py`. +- `ceil.txt:N` and `floor.txt:N` are the frozen `h3_band_ceil.py.txt` and `h3_band_floor.py.txt`, in the frozen directory named at README:11. The frozen `H3_BAND_CEIL.md` and `H3_BAND_FLOOR.md` sit beside them. + +**Provenance of the sections.** The four scope sections were written and adversarially checked against an earlier program, sha256 `99eceb98…`. The final program differs from it by fail-closed guards, comments, and record checks that never feed the computation. Its certified decimals and both transcripts are byte-identical: transcript sha256 `d0def90d…` (ceil) and `c80b7a0f…` (floor). The differences: +- `ck(up(var_iv) > 0)` in `t0_lower` (414); +- `ck(x.b < 1)` in `gram_tail` (884); +- `ck(var.a > 0)` in `trunc_mom_iv` (890); +- `ck(oneside is not None or a == -b)` in `piece_outside`, before the two-sided tail (1234); +- `structure` mode under `--mutation`, which now prints `H3_BAND_V2 MUTATED STRUCTURE RUN (NAME): …; NOT A CERTIFICATE` and exits 3 (1451-1454); +- the comments at 1149-1151 and 1203; +- the record checks: + - the runtime gate in `compare_record` (1406-1427); + - the producer pin `MPMATH_PIN` (1366); + - the transcript-header mode `transcript_mode` (1369-1375); + - `verify-record`'s checks for an unknown header and a non-object record (1458-1477); + - the replay's record check before the run (1484-1488). + + These read only records and transcripts, so they change no certified decimal. + +Every line reference below was remapped to the final file and re-read there. Where a section had described one of these guards as missing, it now describes the final code. + +## How to use this guide + +### What the run checks by itself + +- **Run commands:** README:97-113. They cover `structure`, `ceil --verify`, `floor --verify` and the unit tests. A run prints one summary line, and the exit codes are listed at README:109-112. `structure` needs only the standard library (the import guard at 38-42). +- **Fail-closed checks.** Every `ck` (101-104) that fails stops the run with `H3_BAND_V2 FAIL: …` (in `structure` mode, `STRUCTURE FAIL`) and exit 1. No PASS line is printed, and no transcript or record is written: those are written only after `run_certificate` returns (1489-1501). +- **Replay.** `--verify` requires the fresh record to equal the reviewed one: the program digest, the cov_exact digest, the mpmath version, the transcript digest and every certified decimal (`compare_record`, 1406-1427). `verify-record` compares the same fields, without recomputing, so its outcome does not depend on the verifier's mpmath version, and it needs only the standard library (1458-1477). + - **Producer pin.** Every check requires the reviewed record to name the pinned producer mpmath, `MPMATH_PIN` (1366, 1420-1424). A test binds that constant to `requirements.txt`. Only a recomputing replay also compares the verifier's own mpmath with it. + - **Mode.** `verify-record` reads the mode from the transcript's first line (`transcript_mode`, 1369-1375), which the transcript digest binds, so a record relabelled to the other envelope fails. + - **Malformed inputs.** A record that is valid JSON but not an object, and a transcript whose first line names no envelope, give `VERIFY_FAIL` instead of passing or a traceback (1462-1471). The replay checks the record's type before it runs (1484-1488). + - **Known limits.** None of these lets the values the program checks disagree with the transcript or the source: + - A record that Python's `json` module rejects stops with a traceback and exit 1. That covers most invalid JSON, and valid JSON beyond the module's limits, such as an integer of more than 4,300 digits or very deep nesting. The module accepts the non-standard literals `NaN` and `Infinity`, so a record containing them is read and compared, not stopped. A non-UTF-8 transcript stops with a traceback too. + - The transcript is read with universal newlines, and the digest is taken over the text as read. So a copy with CRLF or lone-CR line endings still matches, although a tool that splits only on `\n` may see different lines. + - A record is read last-key-wins. If a key appears twice, only its last value is checked, while a person or another parser may read the first. + - Fields are compared by value, so a record with `"schema": true` or `1.0` matches `1`. + - `make_record` reads the verifier's `mpmath.__version__` whenever the import guard succeeds (38-42, 1397). An mpmath module that provides `mp`, `iv` and `mpf` but no `__version__` stops `verify-record` with a traceback. + + The tests run `verify-record` with no site-packages and with a faked mpmath version. They check that: + - a changed check line is caught by the transcript digest; + - a record naming no or another mpmath, or with a missing or non-object environment, is rejected, by `compare_record` in both paths too; + - a record relabelled to the other mode is rejected, and a header behind other text, behind a blank line or indented names no envelope; + - a non-object record of every JSON type is rejected, and the replay rejects it before running, even with a fake mpmath that cannot compute; + - a replay still rejects a record made with another version (test:131-261, 360-385). +- **Negative controls** (README:115-146). There are fourteen mutations. Twelve are caught inside the run at a named check, and `drop_fail_piece` and `old_remainder` only by the replay. The record controls are listed there too. The tests are at test:74-98 (structure mode) and test:407-463 (full runs). +- **What the controls do and do not show.** They show that each named check can fire, and that the checks are not vacuous against those regressions. They do not show that the inequalities used between the checks are true. Those inequalities are the subject of this guide. + +### What rests on a written derivation + +| scope | claim | entry points (final lines) | imports | +|---|---|---|---| +| 1. CB2s pole lemma and the first cell | the Gram Laurent series have only the poles (a)-(b) allow; μ and Σ have no pole, exactly | `structure_lemma` 212-237 (tables 154-170, `laurent_expand` 173-183, `solve_in_span` 186-209), run at 268 and by `structure` (1444-1457); the Gram series 460-568; used by `first_cell_law` 1118-1122 | CB1 moments | +| 2. CB4 R-majorant induction | on (0, 2⁻¹⁴], μ and Σ lie within MT and ST of their order-60 Taylor polynomials, and G(t) is invertible | 705-885: A_n(ρ) 706-726, block sups 728-743, θ 778-783, R 787-800, `Maj` 806-851, `mu_tail`/`Sig_tail` 861-873; used by `first_cell_law` 1118-1122 | scope 1 | +| 3. B2 Neumann remainder | on the seven Neumann cells, G(t) is invertible, and μ, Σ lie within T_mu, T_Sig of the same polynomials | `neumann_law` 1048-1116; `gram_tail` 882-885; the CB3 polynomials 617-676 and 1042-1046 | scope 1; the block sups of scope 2 (706-742) | +| 4. Envelope algebra | given the cell law, E[G_r] ≤ ⌈hi⌉ and E[G_r] ≥ lo.a | `ceil_hi` 1148-1198, `band_lo` 1200-1247; helpers 397-433, 888-998, 1132-1146 | scope 2 on the first cell, scope 3 elsewhere | + +**Lanes.** The scopes can be reviewed separately, one lane per scope. Each section states the hypotheses it imports from the others and does not re-prove them. The sections below follow dependency order, 1 to 4. + +**Outside all four scopes:** +- the definition of G_r and the reduction to it (README:16, README:91); +- CB0 and CB1, including CB1d; +- `iv_Phi` itself; +- the CB7 cross-check; +- trust in mpmath `iv` (README:88). + +## Scope 1. CB2s: the pole-structure lemma and its use on the first cell + +### 1.1 Notation +- **The field.** f is the centred Gaussian field with covariance K(x−y) = K1(x1−y1)·K1(x2−y2), where K1(s) = Z⁻¹ Σ_{k∈(π/12)ℤ} e^{−k²/2} e^{iks}. K and K1 are defined in the module docstring, cov_exact.py:9-13. + - The covariance by multi-index is `cov`, cov_exact.py:417-431: Cov(∂^α f(x), ∂^β f(y)) = (−1)^{|β|} K1^{(α1+β1)}(dx) K1^{(α2+β2)}(dy), with (dx, dy) = x − y. + - K1 is entire, because Σ_k e^{−k²/2} e^{|k|R} < ∞ for every R. +- **Kernel derivatives at 0.** κ(n) := K1^{(n)}(0). It is 0 for odd n and (−1)^{n/2} a_n for even n, where a_n = Z⁻¹ Σ e^{−k²/2} kⁿ are the torus moments. `A2[m]` encloses a_{2m} (comment 293); that is CB1, outside this scope. +- **z** is the real separation r. It also serves as the Laurent variable. +- **Stations.** M = 0 and S = z·e_x (comment 150-151). `DISP[(a,b)]` (541) is x_a − x_b in units of z·e_x. +- **Table terms.** A term (st, α, w, k) (154-169) is the random variable w·z^k·∂^α f(x_st). A functional is a list of terms. +- **Formal space 𝓛.** 𝓛 is the ℚ-span of the formal symbols ∂^β (β ∈ ℕ²). Each symbol stands for ∂^β f(0). +- **Bilinear form B.** B(∂^β, ∂^γ) = (−1)^{|γ|} κ(β1+γ1) κ(β2+γ2), which is Cov(∂^β f(0), ∂^γ f(0)). B is symmetric: a nonzero value needs both index sums even, so |β| ≡ |γ| mod 2. +- **Formal Laurent coefficients F^{(p)} ∈ 𝓛.** An M-term adds w·∂^α at p = k. An S-term adds (w/q!)·∂^{α+(q,0)} at p = k+q, for each q ≥ 0 (Taylor expansion in x1). This is `laurent_expand` (173-183), which is complete for every p ≤ `top`. +- **The tables.** + - V = (V0..V5) are the pins, `PINS_EXACT` (154-161). + - Y = (Y0..Y5) = (q_M, q_S, a_M, a_S, b_M, b_S) are the soft variables, `SOFT_EXACT` (162-169). + - zV = (6/5, 0, 0, 0, 0, −1/6), `ZV_EXACT` (170), is the conditioning value, a constant. **zV is a name, not z·V.** +- **Gram series.** G = Cov(V,V), A = Cov(Y,V) and GY = Cov(Y,Y) are `GP_ser` (551), `GYV_ser` (553) and `GY_ser` (552). Their Laurent coefficients are G_p, A_p and GY_p. +- **Derived series.** + - H = G⁻¹ = Σ H_j z^j, w = H·zV = Σ W_j z^j, μ = A·w = Σ X_k z^k, and Σ = GY − A·H·Aᵀ = Σ D_k z^k (624-664). + - `mu_ser` = X[0..J] (648) and `Sig_ser` = D[0..J] (664), with J = 60 and NK = J+2 (455-456). +- 𝔼 denotes expectation; E denotes the matrix below. + +### 1.2 The tables, expanded (each line can be checked by hand) +- **Pins.** + - V0 = ∂00, V1 = ∂10 and V2 = ∂01. + - V3 = (S10−M10)/z = ∂20 + (z/2)∂30 + O(z²), and V4 = (S01−M01)/z = ∂11 + (z/2)∂21 + O(z²). + - V5 = (S00−M00)/z³ − (S10+M10)/(2z²). Its two parts are (S00−M00)/z³ = z⁻²∂10 + z⁻¹∂20/2 + ∂30/6 + z∂40/24 + … and −(S10+M10)/(2z²) = −z⁻²∂10 − z⁻¹∂20/2 − ∂30/4 − z∂40/12 − …. So V5 = −∂30/12 − (z/24)∂40 + …, and the z⁻² and z⁻¹ terms cancel. +- **Soft variables.** Y0 = ∂02 and Y1 = ∂02 + z∂12 + …; Y2 = z⁻¹∂20 and Y3 = z⁻¹∂20 + ∂30 + …; Y4 = z⁻¹∂11 and Y5 = z⁻¹∂11 + ∂21 + …. +- **Hence:** + - V^{(0)} = (∂00, ∂10, ∂01, ∂20, ∂11, −∂30/12). These are linearly independent in 𝓛, so E is unique. The proof does not need uniqueness: any E satisfying (c) and (d) works. + - Y^{(−1)} = (0, 0, ∂20, ∂20, ∂11, ∂11). + - E has rows (0, 0, e3, e3, e4, e4), with columns 0-indexed; rows 0 and 1 are zero. + - E·zV = (0, 0, zV3, zV3, zV4, zV4) = 0. + +`h3_band_v2.py structure` prints exactly this E (1455-1456). A re-run on the final program printed `E = [0 0 0 0 0 0] [0 0 0 0 0 0] [0 0 0 1 0 0] [0 0 0 1 0 0] [0 0 0 0 1 0] [0 0 0 0 1 0]`. + +### 1.3 Exact statement +**Hypotheses.** `structure_lemma` (212-237) decides each of them in `Fraction` arithmetic, on the tables that `mutated_tables()` returns (66-82). It calls `laurent_expand(·, 1)`, which is complete for p ≤ 1. +- **(a)** V_k^{(p)} = 0 for every p < 0 (223-225). +- **(b)** Y_a^{(p)} = 0 for every p < −1 (226-228). +- **(c)** Y_a^{(−1)} = Σ_k E_{ak} V_k^{(0)} in 𝓛, for an exact E ∈ ℚ^{6×6} (229-234). + - `solve_in_span` (186-209) is exact Gauss–Jordan elimination. It returns None if and only if a row of the form 0 = c with c ≠ 0 is left (204-205). Otherwise it returns the solution with the free variables set to 0. + - The returned E is not substituted back at run time, so its correctness rests on this short routine. + - Two things cover it. An independent substitution (pre-checks, scope 1) confirms E. Exactness rests on that rational identity. The numeric check at 691-695 (observed |A_{−1} − E·G0| = 8.9e−100, threshold 1e−80) is corroboration only: it catches a grossly wrong E, but a perturbation of E below about 1e−80 that still satisfies E·zV = 0 would pass it (scope-1 review 5607220361). +- **(d)** E·zV = 0 exactly, row by row (235-236). + +**Extra hypothesis (e): G_0 is invertible.** This is not part of CB2s; CB3 proves it at 612-614. +- The check bounds every entry of the iv product G0·Ginv0 − I by 1e−70; observed 8.26e−98. +- So ‖GX−I‖_∞ ≤ 6e−70 < 1 for every real G in the box G0 and every real X in the box Ginv0, and every such G is invertible. +- The true G_0 is among them provided it lies in the box G0. That rests on Step 1 below and on CB1's moment enclosures. +- If a pivot interval in `iv_solve6` contains 0, mpmath iv division returns an unbounded interval (checked under mpmath 1.4.1): [−inf, inf] when 0 is interior or the divisor is [0, 0], and a half-line such as [1/b, inf] when 0 is an endpoint. Either way `ub` of an affected entry is inf, and the check fails closed. + +**Conclusions.** +- **(i)** G_p = 0 for p < 0; A_p = 0 for p ≤ −2; GY_p = 0 for p ≤ −3. Also A_{−1} = E·G_0 and GY_{−2} = E·G_0·Eᵀ. +- **(ii)** Assuming (e), μ and Σ are analytic at 0: X_{−1} = 0 and D_{−2} = D_{−1} = 0. The result for μ needs (a)-(e); the result for Σ needs only (a)-(c) and (e). + +### 1.4 Proof +**Step 1: the formal expansion and the Gram series agree.** For one pair of terms, Cov(∂^α f(x_a), ∂^β f(x_b)) = (−1)^{|β|} K1^{(α1+β1)}(d·z)·κ(α2+β2), where d = DISP ∈ {−1, 0, 1}. In the code this is `cov_ser` (520-534). Since K1 is entire, K1^{(n)}(dz) = Σ_q κ(n+q)(dz)^q/q!. `kernel_ser` (500-511) stores coefficient j as κ(n+j)/j!, skips odd n+j, and flips the sign of odd j when d = −1. The three cases: +- **(S,M):** the term is Σ_q (z^q/q!)·B(∂^{α+(q,0)}, ∂^β). +- **(M,S):** the term is Σ_q (z^q/q!)·B(∂^α, ∂^{β+(q,0)}), because the (−1)^q comes from |β+(q,0)|. +- **(S,S) or (M,M):** d = 0, and the term is the constant B(∂^α, ∂^β). The formal double sum gives the same value, since Σ_{q+q'=N} C(N,q')(−1)^{q'} = 0 for N ≥ 1. + +Multiply by w_a z^{k_a} · w_b z^{k_b} and sum over term pairs (`gram_entry_ser`, 543-548). For every pair of functionals this gives Cov(F_a, F_b) = Σ_p z^p Σ_{i+j=p} B(F_a^{(i)}, F_b^{(j)}), which converges absolutely for z ≠ 0. By bilinearity, the same holds for any rational combination of functionals. + +**Step 2: conclusion (i).** +- G_p = Σ_{i+j=p} B(V^{(i)}, V^{(j)}) is 0 for p < 0, by (a). In particular G_0 = B(V^{(0)}, V^{(0)}). +- In A_p, (b) forces i ≥ −1 and (a) forces j ≥ 0, so A_p = 0 for p ≤ −2. Then A_{−1} = B(Y^{(−1)}, V^{(0)}) = E·B(V^{(0)}, V^{(0)}) = E·G_0, by (c). +- In the same way, GY_p = 0 for p < −2, and GY_{−2} = B(Y^{(−1)}, Y^{(−1)}) = E·G_0·Eᵀ. + +**Step 3: conclusion (ii).** Let Ỹ := Y − E·V/z. Then Ỹ^{(p)} = Y^{(p)} − E·V^{(p+1)}. For p ≤ −2 this is 0 by (b) and (a); for p = −1 it is 0 by (c). Applying Step 1 to Ỹ, both à := Cov(Ỹ, V) = A − E·G/z and GỸ := Cov(Ỹ, Ỹ) have no negative powers. G is entire with no negative orders, so by (e) H = G⁻¹ is analytic near 0. +- **For μ:** μ = (à + E·G/z)·H·zV = ÷H·zV + E·zV/z = ÷H·zV, by (d). +- **For Σ:** using GH = HG = I and G = Gᵀ, GY = GỸ + (÷Eᵀ + E·Ãᵀ)/z + E·G·Eᵀ/z² and A·H·Aᵀ = ÷H·Ãᵀ + (÷Eᵀ + E·Ãᵀ)/z + E·G·Eᵀ/z². So Σ = GỸ − ÷H·Ãᵀ. + +Both right-hand sides are analytic at 0. Laurent coefficients on a punctured disc are unique, so X_{−1} = 0 and D_{−2} = D_{−1} = 0. In coefficients, these are the quantities the code evaluates: +- X_{−1} = A_{−1}·W_0 = E·G_0·G_0⁻¹·zV = E·zV = 0. This is X[−1] at 642-647, where k = −1 admits only i = −1. +- D_{−2} = GY_{−2} − A_{−1}·H_0·A_{−1}ᵀ = E·G_0·Eᵀ − E·G_0·Eᵀ = 0. D[−2] and D[−1] are computed at 658-663. +- D_{−1} involves G_1 through H_1 = −G_0⁻¹G_1G_0⁻¹. G_1 ≠ 0 for these tables, so this cancellation is not trivial. + +This completes the proof. + +**Probabilistic reading.** For a centred Gaussian, 𝔼[Y | V = v] = A·G⁻¹·v and Cov(Y | V) = GY − A·G⁻¹·Aᵀ; this is standard conditioning. The docstring (218-220) states the resulting conclusion: 𝔼[Y | V = zV] = 𝔼[Ỹ | V = zV] + E·zV/z has no pole, and Cov(Y | V) = Cov(Ỹ | V). + +### 1.5 The lemma is about the objects the series use +- **Same objects.** At 265, `run_certificate` takes `pins_exact`, `soft_exact` and `zv_exact` from `mutated_tables()` and passes them to `structure_lemma` (268). From those same lists it builds `PinsV` and `SoftY` (539-540, through `to_iv_func` 536-537: the power is kept, and the weight becomes an iv enclosure of the rational) and `zV` (609). `E_sel` is used only at 694. +- **Same geometry.** `laurent_expand` places S at z·e_x (175, 181-182). DISP[('S','M')] = +1 (541) means x_S − x_M = +z·e_x, which matches `cov`'s convention (dx, dy) = x − y. +- **Complete kernel window.** + - For a cross-station pair of terms, an order-p coefficient needs kernel order j = p − k_a − k_b, and `kernel_ser` keeps j ≤ J+8−n1 (502). + - Over all cross-station term pairs in the tables, the maximum of n1 − k_a − k_b is 6. It is attained by V5×V5, V5×Y2, V5×Y3 and Y2×Y3. So every coefficient with p ≤ NK = J+2 is complete (comment 497-499). + - The lowest orders present are −6 (GP), −4 (GYV) and −2 (GY). All are ≥ LMIN = −8 (457), so `ser_trim` (493-495) removes nothing at the bottom. Above NK it removes incomplete orders, which the Cauchy tails replace. +- **Same as the frozen tables.** The tables equal the frozen PinsV, SoftY and zV of both ceil.txt and floor.txt (lines 172-188 and 237), compared programmatically. +- **CB7 is independent.** CB7 re-encodes the functionals by hand (1281-1305), separately from the tables, with M at −r/2 and S at +r/2; the covariance is translation-invariant. It is a cross-check at r = 1/20 against the Neumann law, not load-bearing. + +### 1.6 How the conclusions are used +**C1: every cell, using (i), that is (a) and (b).** +- For the law, the code reads only GP orders 0..NK (608, 617), GYV orders −1..NK (619-620) and GY orders −2..NK (618, 621-622). +- The GP orders {−6, −4, −3, −2, −1} and the GYV orders {−4, −3, −2} are computed (|c| ≤ 6.3e−100), checked at 555-567, and discarded, because their true values are 0. +- The recursions for H and W (624-634) and θ(R), which sums from i = 1 (778-783), solve G·H = I order by order. They are valid only because G has no negative orders. +- The Neumann law uses G_P from order 0 (1067) and A_P = A_{−1}/t + Σ_{j≥0} A_j t^j (1094-1095). + +**C2: the first cell only, using (ii).** +- `first_cell_law` (1118-1122) returns, entrywise, μ = Σ_{k=0}^{J} X_k t^k ± MT and Σ = Σ_{k=0}^{J} D_k t^k ± ST on t ∈ [0, up(1/16384)] = [0, 2⁻¹⁴]. It omits X_{−1}/t, D_{−2}/t² and D_{−1}/t. +- t can be arbitrarily small on (0, 1/16384], so any nonzero pole coefficient, however small, would make these terms unbounded. The omission is valid only because the pole coefficients are exactly 0. +- The heavy coefficients X_k and D_k for k ≤ J are complete: they need A up to J+1, W up to J+1 and H up to J+2 = NK. +- MT and ST (861-878) bound only the orders k > J. They include the pole-index terms A_{−1}·W_{k+1} and so on (864, 871-872). +- The floor's rotated law (1205-1206) is LR applied to the same series (681-683), so its pole coefficients LR·X_{−1} and so on also vanish. + +**C3: the Neumann cells do not use (c) or (d).** They add |X[−1]|/r_lo and |D[−2]|/r_lo² + |D[−1]|/r_lo (1099, 1108-1109), with r_lo > 0 checked at 1063. They do rely on C1. + +### 1.7 What the program checks at run time +- **(a)-(d), exactly**, at 268, failing closed through `ck` (101-104). + - There are four controls (72-79), tested at test:80-89 (structure mode) and test:411-414 (full runs). Each exits 1 at its own branch; `structure_window` fails at row 2. + - In `structure` mode, a mutation the lemma cannot see now ends `H3_BAND_V2 MUTATED STRUCTURE RUN (…): the lemma holds for these tables; NOT A CERTIFICATE` and exits 3 (1451-1454; test:91-98). Such mutations are `wrong_window` and those that do not touch the tables. For `wrong_window` the lemma really does hold on the mutated tables, because column 5 of E is zero. +- **(e):** the G0 inverse residual is below 1e−70 (612-614); observed 8.26e−98. +- **Gram pole pattern (555-567).** Every coefficient that `ser_trim` kept (orders LMIN..NK) and that the law then discards has |enclosure| ≤ 1e−60; observed ≤ 6.3e−100. `ser_trim` would remove any order below LMIN = −8 before this check, so the check never sees such orders. The tables have none (lowest −6), and the exact lemma runs first. +- **Law pole enclosures:** |X[−1]| < 1e−55 (687), and |D[−2]|, |D[−1]| < 1e−50 (690). Observed: 3.4e−99, 8.3e−99 and 1.0e−98. Each enclosure contains 0. +- **A_{−1} = E·G0 entrywise within 1e−80 (691-695).** This ties the exact E to the interval series; observed 8.9e−100. Conclusion (i)'s GY_{−2} = E·G_0·Eᵀ is not checked at run time, because the law does not need it; a scratch run observed 8.9e−100. +- The program has no run-time cross-check of its own for the first-cell law. CB7 runs at r = 1/20 only. + +### 1.8 What rests only on the derivation +1. **Step 1 (500-548).** The interval Gram Laurent coefficients equal B applied to `laurent_expand`'s formal coefficients. This rests on four things: K1 being entire, the DISP sign convention (541), the extra (−1)^q in the (M,S) case, and the cancellation of the S-S double sum. At run time the evidence is only numeric: the pole pattern (555-567), A_{−1} = E·G0 to 1e−80 (691-695), the CB3 limit identities (697-702), and CB7 at r = 1/20. +2. **Steps 2-3.** The algebra by which (a)-(d) plus (e) imply (i) and (ii), through Ỹ = Y − E·V/z, is not checked anywhere. The numeric checks at 555-567, 687 and 690 agree with it, but they cannot tell 0 from a coefficient below 1e−60, 1e−55 or 1e−50 respectively. +3. **That every discarded coefficient is exactly zero.** This covers GP orders < 0 and GYV orders ≤ −2 in every cell, since the law reads only orders 0..NK, −1..NK and −2..NK (608, 617-622). It also covers X_{−1}, D_{−2} and D_{−1} in the first cell (`first_cell_law`, 1118-1122). There, any nonzero pole coefficient would make the enclosure invalid for small enough t. +4. **That (c) holds for the E that `solve_in_span` returns.** E is not substituted back at run time; 691-695 corroborates it numerically. +5. **That the H and W recursions (624-634) and θ(R) (778-783, summing from i = 1) are the true order-by-order inverse.** This needs G to have no negative orders, which is lemma (a). +6. **That the kernel window is complete through order NK** (comment 497-499, range at 502). The maximum of n1 − k_a − k_b over cross-station pairs is 6. No run-time check would detect a missing kernel term at orders near NK. +7. **That the true G_0 lies in the box G0**, so that the residual check at 612-614 proves (e). This depends on Step 1 and on CB1's moment enclosures. +8. **That μ = A·G⁻¹·zV and Σ = GY − A·G⁻¹·Aᵀ are the conditional law of Y given V = zV.** This is centred-Gaussian conditioning with zV a constant vector, as in CB7's zVd (1301). +9. **That the tables (154-170) are the upstream corrected pins, soft variables and conditioning value that define G_r.** README:16 and README:91 say this is not re-derived. CB3's limit identities (697-702) and CB7's hand-coded functionals (1284-1289) are only partial evidence. +10. **That the first-cell Taylor series of μ and Σ converge to the true law on (0, 1/16384], with tails MT and ST (857-878).** That is scope 2; this lemma supplies only the absence of negative orders. + +## Scope 2. CB4: the R-majorant induction and the first-cell tails + +This scope covers 705-885 and the places where the first cell uses those results: 1118-1122, 1256-1258, 1264 and 1269, plus the 2MT/4ST widening at 1205-1206. + +### 2.1 Notation +- **Series variable.** The program identifies the series variable z with the cell variable r. `first_cell_law` builds t_iv = [0, r*] (1119), and the cell loop passes t_iv as r to `ceil_hi` and `band_lo` (1257, 1264, 1269). Here r* = R_STAR = 1/16384 = 2⁻¹⁴ (857), and t ranges over (0, r*]. +- **Orders.** J = 60 and NK = J+2 = 62 (455-456). +- **The circle.** ρ = RHO = mpf('0.11') at dps 100 (458). + - RHO is a dyadic rational, and RHO − 0.11 ≈ 2.9e−103. Every CB4 use reads this one mpf, so the Cauchy circle is exactly |z| = RHO. + - QRHO = up(1/RHO) ≥ 1/RHO (769-770). +- **Gram blocks**, all 6×6: G(z) is the pin Gram, A(z) the soft–pin cross Gram, and GY(z) the soft Gram. + - Their Laurent ranges, G = Σ_{j≥0} G_j z^j, A = Σ_{j≥−1} A_j z^j and GY = Σ_{j≥−2} GY_j z^j, are imported from scope 1, (a) and (b) (see H1). + - The program's enclosures are `Gj[j]`, `Aj[j]` and `GYj[j]` for 0 ≤ j ≤ NK, plus `A_m1`, `GY_m1` and `GY_m2` (617-622). + - `mat_coeff` (571-572) reads only these orders, so the program drops the numeric coefficients at every other order. + - `Ac(i)` returns A_{−1} for i = −1 and A_i otherwise (636-637). +- **Inverse and w.** H(z) = G(z)⁻¹ = Σ_{j≥0} H_j z^j and w(z) = H(z) zV = Σ W_j z^j. zV is ZV_EXACT = (6/5, 0, 0, 0, 0, −1/6) (170); the program's `zV` is its interval enclosure (609). +- **The law.** μ(z) = A(z) w(z) and Σ(z) = GY(z) − A(z) H(z) A(z)ᵀ. + - Laurent coefficients: μ_k = Σ_{i=−1}^{k} A_i W_{k−i}, and Σ_k = GY_k − Σ_{a+b+c=k} A_a H_b A_cᵀ with a, c ≥ −1 and b ≥ 0. + - Enclosures: X[k] and D[k] (641-664). `mu_ser` = X[0..J] (648) and `Sig_ser` = D[0..J] (664). +- **Norms and helpers.** ‖P‖_max = max_{a,b} |P_ab| and ‖v‖_max = max_a |v_a|. `ub(x)` = max(|lo|, |hi|) (241-245) and `up(x)` is the upper endpoint (248-250). `matmax` and `vecmax` return the maximum of `ub` over the entries (745-753). + +### 2.2 Exact claim +For every t ∈ (0, r*] and all indices a, b: + + μ_a(t) ∈ Horner(X[0..J]_a, [0, r*]) + [−MT, MT], Σ_ab(t) ∈ Horner(D[0..J]_ab, [0, r*]) + [−ST, ST]. + +- MT = `mu_tail(tS)` and ST = `Sig_tail(tS)` (875-876). Horner is interval Horner (1001-1008), applied at 1118-1122. +- The transcript prints MT = 1.04e−44 and ST = 1.84e−38, as 3-digit round-to-nearest `nstr` values. The certificate uses the full mpf upper endpoints, 1.03890e−44 and 1.83842e−38. +- G(t) is invertible for 0 ≤ t < 1/R, which contains [0, r*]. + +### 2.3 Hypotheses imported from other scopes +- **H1 (scope 1, exact in `Fraction`; 212-237, run at 268 before any series is built).** It proves four facts: + - (a) the pins have no negative powers, so G(z) is a power series; + - (b) the soft functionals have no power below z⁻¹, so A has at most a z⁻¹ pole and GY at most a z⁻² pole; + - (c) Y^{(−1)} = E V^{(0)} for an exact matrix E; + - (d) E zV = 0. + + How CB4 uses them. (a) and (b) fix the index ranges: the H and W recursions with G0 as leading coefficient, the θ sum from i = 1, i ≥ −1 in μ_k, and the a, c ≥ −1 pole groups in Σ_k. (c) and (d) give μ_{−1} = Σ_{−2} = Σ_{−1} = 0 exactly, which is why the first-cell law drops X[−1], D[−2] and D[−1]. + + A direct algebraic check of the three cancellations, with C_pq = cov(V^{(p)}, V^{(q)}): + - A_{−1} = cov(Y^{(−1)}, V^{(0)}) = E G0, hence μ_{−1} = A_{−1} G0⁻¹ zV = E zV = 0. + - Σ_{−2} = GY_{−2} − A_{−1} G0⁻¹ A_{−1}ᵀ = E G0 Eᵀ − E G0 Eᵀ = 0. + - A_0 = cov(Y^{(−1)}, V^{(1)}) + cov(Y^{(0)}, V^{(0)}) and H_1 = −G0⁻¹(C01 + C10)G0⁻¹. Then A_{−1}H_0A_0ᵀ + A_0H_0A_{−1}ᵀ + A_{−1}H_1A_{−1}ᵀ = GY_{−1}, so Σ_{−1} = 0. + + The docstring's route (218-220) gives the same result at all orders. Ỹ = Y − EV/z is regular by (a)-(c), Cov(Y|V) = Cov(Ỹ|V), and E[Y|V=zV] = E[Ỹ|V=zV] + E zV/z, where the last term vanishes because E zV = 0. The numeric checks at 555-567 (Gram pole pattern), 687-690 (poles of μ and Σ tiny) and 692-695 (A_{−1} = E G0) are consistency checks only. +- **H2 (CB2 window).** The Gram enclosures through order NK enclose the true coefficients. The kernel series uses `range(0, J+9-n)` (502), argued only in the comment at 497-499. + - Only cross-station pairs (d = ±1) use the truncated kernel series; same-station pairs are constants (525-528). + - A cross pair reaches order k through kernel index j = k − k_a − k_b, and `kernel_ser(n1)` keeps j ≤ J+8−n1. So order NK = J+2 is complete exactly when (−k_a + α_{a,x}) + (−k_b + α_{b,x}) ≤ 6. + - The per-term weight −k + α_x is at most 3 in the tables (154-169); pin V5's four terms and soft Y2 and Y3 reach 3. So the condition holds. + - No run-time check ties the window to the tables. A mutation that would exceed it (`structure_soft_pole2`, weight 4) is stopped by CB2s (b) before CB2 runs. +- **H3 (CB3).** `Ginv0` encloses G0⁻¹ (610-614). H[j] and W[j] for j ≤ NK enclose the true H_j and W_j (624-634). X[k] and D[k] for k ≤ J are complete: they use W only up to k+1 ≤ NK and H only up to k+2 ≤ NK (641-663). +- **H4 (covariance formula).** Each Gram entry is f(z) = Σ_pairs w_a w_b z^{k_a+k_b} (−1)^{|β|} K1^{(n1)}(dz) K1^{(n2)}(0), with d ∈ {−1, 0, 1} (520-548). Here K1(x) = Z⁻¹ Σ_{k∈(π/12)ℤ} e^{−k²/2} cos(kx), Z is the full lattice sum, and K1 is entire. The program evaluates the Gram through the Poisson moment table A2, which CB1d cross-checks against the spectral form. +- **H5.** mpmath `iv` rounds outward correctly. + +### 2.4 Derivation +1. **Kernel sup (706-726).** For complex w with |w| ≤ ρ, K1^{(n)}(w) = Z⁻¹ Σ e^{−k²/2} (ik)^n e^{ikw}, and |e^{ikw}| ≤ e^{|k||Im w|} ≤ e^{|k|ρ}. Hence |K1^{(n)}(w)| ≤ A_n(ρ) := Z⁻¹ Σ e^{−k²/2} |k|^n e^{|k|ρ}. + - **Lattice part:** |j| ≤ Jl = 116, where Jl = ceil(30/(π/12)) + 1 (cov_exact.py:152). |k| is replaced by its upper endpoint, and k = 0 contributes 1 because iv 0**0 = 1 (716-717). + - **Tail:** |j| ≥ j0 = Jl+1, where k0 = j0·h ≈ 30.63. For |k| ≥ 4ρ, −k²/2 + |k|ρ ≤ −k²/4. Put u_j = (jh)^n e^{−(jh)²/4}; then u_{j+1}/u_j = e^{−(2j+1)h²/4}(1+1/j)^n. This ratio decreases in j, so for j ≥ j0 it is at most q (721). The tail is therefore at most 2u_{j0}/(1−q) (723); the factor 2 covers both signs of k. + - **Normalization:** the program divides by Z1_iv, an interval enclosing the truncated lattice sum. Its lower endpoint is at most the exact truncated sum, which is at most Z. So the upper endpoint of the interval quotient bounds the quotient by Z (724). Only up(A_rho[n]) is consumed downstream (scope-2 review 5607205948). + - A_rho[0..8] is computed (726), but only n1, n2 ≤ 4 occur: n1 = 4 from soft (2,0)+(2,0), and n2 = 4 from soft (0,2)+(0,2). +2. **Block sup and Cauchy bound (728-743).** + - On |z| = ρ the triangle inequality gives |f(z)| ≤ Σ_pairs |w_a||w_b| ρ^{k_a+k_b} A_{n1}(ρ) A_{n2}(ρ). This uses |K1^{(n)}(0)| ≤ A_n(ρ), which also covers the d = 0 pairs. + - The block constants are the maxima over the 36 entries: M_GP = 2.9579e6, M_GY = 343.659 and M_GYV = 19473.5. + - Each f is a finite sum of Laurent monomials times entire functions, so it is analytic on ℂ\{0}. Cauchy's estimate on |z| = ρ then gives |f_j| ≤ M ρ^{−j} ≤ M·QRHO^j for every j ∈ ℤ, negative j included. +3. **Majorant sequences (755-761).** + - g_i = G_m[i] for 1 ≤ i ≤ J, and M_GP ρ^{−i} for i > J. θ uses the Cauchy bound even at J+1 and J+2, where computed values exist; that is valid and only conservative. + - a_i = A_m[i] for 0 ≤ i ≤ NK, and M_GYV·QRHO^i beyond. y_i = GY_m[i] for 0 ≤ i ≤ NK, and M_GY·QRHO^i beyond. + - α = A1_m ≥ ‖A_{−1}‖_max (α = 3), and γ = Ginv0_m ≥ ‖G0⁻¹‖_max (γ = 24). + - By H2, H3 and step 2, each sequence dominates the true ‖·‖_max. The B1 guard (762-767) checks matmax ≥ ub(entry); in unmutated code this holds by definition. +4. **The constant 36 and the induction.** For 6×6 matrices, ‖PQ‖_max ≤ 6‖P‖_max‖Q‖_max and ‖Pv‖_max ≤ 6‖P‖_max‖v‖_max. + - **Recursions.** From GH = I, H_j = −G0⁻¹ Σ_{i=1}^{j} G_i H_{j−i} for j ≥ 1 (628-632). From Gw = zV, W_j = −G0⁻¹ Σ_{i=1}^{j} G_i W_{j−i} (633-634). + - **Where 36 comes from.** Put h_j := ‖H_j‖_max. The first 6 is for G0⁻¹ and the second for G_i H_{j−i}: h_j ≤ 6γ Σ_i 6 g_i h_{j−i} = 36γ Σ_{i=1}^{j} g_i h_{j−i}. The same holds for w_j := ‖W_j‖_max. + - **θ.** θ(R) = 36γ [Σ_{i=1}^{J} G_m[i] R^{−i} + M_GP x^{J+1}/(1−x)], with x = 1/(ρR) < 1. This is exactly `theta` (778-783). + - **Constants.** C_H = max(1, max_{j≤NK} up(H_m[j]/R^j)), and C_W is defined the same way from W_m (802-804). + - **Induction hypothesis:** h_j ≤ C_H R^j and w_j ≤ C_W R^j for every j ≥ 0. + - **Base, j ≤ NK:** h_j ≤ H_m[j] ≤ C_H R^j. The second inequality holds because the iv quotient contains H_m[j]/R^j. + - **Step, j > NK (strong induction):** h_j ≤ 36γ Σ_{i=1}^{j} g_i C_H R^{j−i} ≤ C_H R^j θ(R). The same bound holds for W. + - So **θ(R) ≤ 1 suffices.** The program checks the stronger θ(R) ≤ 1/2, together with ρR > 1, in interval arithmetic at 800. The 1/2 is inherited from ceil.txt:364-381. + - **Search.** The search (787-799) updates R in round-to-nearest mp. R is just an exact mpf, and only the iv check at 800 on the final R is load-bearing. + - **Consequence.** Σ H_j t^j converges absolutely for |t| < 1/R. Its Cauchy product with G(t) is I, and G is entire. So G(t) is invertible for 0 ≤ t < 1/R, with G(t)⁻¹ = H(t) and w(t) = G(t)⁻¹ zV. The first cell has no other invertibility argument. +5. **The numbers.** + - γ = 24, G_m[1] = 1.5 and α = 3. The leading term of θ is 36·24·1.5/R = 1296/R. + - The program finds R = 2622.3492, ρR = 288.46 and θ = 0.49437, of which the leading term is 0.49421. + - The smallest R with θ ≤ 1 is 1296.83, and with θ ≤ 1/2 it is 2592.83. The 2% grid then returns 2622.35. + - C_W = 4 = W_m[0] and C_H = 24 = H_m[0]; both maxima are attained at j = 0. The next-largest ratios are 3.4e−4 for W and 1.1e−3 for H. +6. **`Maj`, P2 and P3 (806-851).** + - A `Maj` sequence takes the computed values c_j for j ≤ N = NK, and C q^j beyond. + - The suffix sums T[m] = Σ_{j≥m} (term)·t^j are sums of nonnegative terms, with T[N+1] = C(qt)^{N+1}/(1−qt). The condition qt < 1 is checked at 824 and 837. + - P2(A, B, m, t) = Σ_{i+l≥m} a_i b_l t^{i+l}; P3 is the triple analogue. + - The instances are majW = (W_m, C_W, R), majH = (H_m, C_H, R), majA = (A_m, M_GYV, QRHO) and majGY = (GY_m, M_GY, QRHO) (853-856). +7. **μ tail (861-865).** For k > J, split off the i = −1 term. Entrywise, |Σ_{k>J} μ_{k,a} t^k| ≤ 6 P2(majA, majW, J+1, t) + 6α Σ_{l≥J+2} w_l t^{l−1}. The second term is `majW.tail_from(J+2, t)/t`. +8. **Σ tail (867-873).** Each entry of A_a H_b A_cᵀ has at most 36 terms. Entrywise, for k > J, there are four groups: + - the GY part, `majGY.tail_from(J+1)`; + - a, c ≥ 0: 36·P3(majA, majH, majA, J+1); + - exactly one of a, c equal to −1. Then b + c = k+1 ≥ J+2 (or a + b = k+1 ≥ J+2 for the mirror term). The two groups have equal sums, giving 2·36·α·P2(majH, majA, J+2)/t; + - a = c = −1. Then b = k+2 ≥ J+3, giving 36α²·`majH.tail_from(J+3)`/t². +9. **Monotone in t.** After the divisions by t and t², every bound in steps 7 and 8 is a power series in t with nonnegative coefficients and exponents ≥ J+1. Its value at tS = 2⁻¹⁴ (859) therefore bounds it at every t ∈ (0, r*]. +10. **First-cell guard (858).** The program checks R·r* < 0.2; the actual value is 0.16005. + - Convergence itself needs only R·r* < 1 and QRHO·r* < 1. R·r* < 1 follows from 858 and is checked again in `tails` for majW and majH. QRHO·r* < 1 is checked only in `tails`, for majA and majGY (824, 837). + - The 0.2 is a margin inherited from the frozen program (ceil.txt:383). + - The frozen first cell (0, 1/400] gives R/400 ≈ 6.56 > 1, so the majorant genuinely diverges there. This is the `wide_first_cell` control (README:48 and README:128; test:418). +11. **Enclosure (1118-1122).** + - μ_a(t) = Σ_{k=0}^{J} μ_{k,a} t^k + (tail), with no k = −1 term by H1. Each μ_{k,a} lies in X[k]_a (H3), and t ∈ [0, r*]. So Horner encloses the polynomial, and ±MT covers the tail. The same argument applies to Σ with ST. + - `band_lo` widens the rotated law by 2MT and 4ST (1205-1206). This is valid because each row of LR (680) has ℓ¹ norm at most 2: a rotated-vector tail is ≤ 2MT and a rotated-matrix tail ≤ 2·2·ST. + - The Neumann law adds numeric pole terms (1099, 1108-1109). The first cell cannot, because 1/t is unbounded there, so it rests on the exact H1. + +### 2.5 Checked numerically at run time +- 720: k0 ≥ 4ρ. 722: the kernel tail ratio q < 1, for each n. +- 767: B1 dominance. This holds by definition in unmutated code; it is a regression guard for `signed_maxima`. +- 776: the Cauchy check, for the computed orders j ≤ NK, for j = −1 (A, GY) and for j = −2 (GY). It is not vacuous: GY_{−2} sits at 0.72 of its bound. The other blocks are loose (GY_0 at 0.022, GYV at 1.5e−4, GP at 1e−6), and it says nothing about j > NK. +- 793: θ(R) ≤ 1/2 after the upward search; not load-bearing. +- 800: ρR > 1 and θ(R) ≤ 1/2 for the final R, in iv; load-bearing. +- 824 and 837: qt < 1 for every `Maj` used. +- 858: R·r* < 0.2. +- 877-878: 0 ≤ ST < 1e−5 and 0 ≤ MT < 1e−5. +- Upstream, consistency checks only: 567 (Gram pole pattern), 614 (G0 inverse residual) and 687-695 (pole consistency). +- **Rounding.** Every load-bearing constant is an iv bound. A_rho, M_*, θ, C_W, C_H, QRHO, MT and ST are upper endpoints; G_m, GY_m, A_m, H_m, W_m, Ginv0_m and A1_m are ub (max |endpoint|). R is an exact mpf, and any value of it is acceptable once 800 passes. + +### 2.6 What rests only on the written derivation +The two items that carry the most weight are the Cauchy estimate beyond the computed orders (j > NK in majA and majGY, and i > J in θ) and the induction step. +- |K1^{(n)}(w)| ≤ A_n(ρ) for complex |w| ≤ ρ, from the spectral form of K1 and |e^{ikw}| ≤ e^{|k|ρ} (706-717). +- **The kernel lattice tail (718-723).** Only k0 ≥ 4ρ and q < 1 are checked. The rest is derivation: + - |k| ≥ 4ρ implies e^{−k²/2+|k|ρ} ≤ e^{−k²/4}; + - the term ratio is monotone in j, so it is bounded by q; + - the tail is ≤ 2u_{j0}/(1−q). +- Dividing by Z1_iv gives an upper bound, because the lower endpoint of Z1_iv is at most the exact truncated sum, which is at most Z (724). +- The triangle-inequality sup bound of each Gram entry on |z| = ρ, including |K1^{(n)}(0)| ≤ A_n(ρ) (728-737). +- **Analyticity and Cauchy's estimate.** Every Gram entry is analytic on ℂ\{0}, so Cauchy's estimate |f_j| ≤ M ρ^{−j} holds for every j. This is load-bearing for j > NK in majA and majGY and for i > J inside θ. The run-time check at 772-776 covers only the computed orders. +- **H1 (a), (b):** the Laurent ranges of G, A and GY. The program drops the numeric coefficients outside them, with only the consistency check at 555-567. +- **H1 (c), (d) ⇒ μ_{−1} = Σ_{−2} = Σ_{−1} = 0.** CB2s itself is exact in `Fraction`. The step from (a)-(d) to the cancellations is the algebra in H1 above, or equivalently the docstring at 218-220. `first_cell_law` relies on it with no numeric pole term (1118-1122). +- **H2:** the kernel window `range(0, J+9-n)` makes every Gram coefficient through NK complete (497-511). There is no run-time tie to the tables. +- **H3:** interval Gauss–Jordan and the series recursions enclose the true G0⁻¹, H_j, W_j, X[k] and D[k] (574-585, 610-664). The G0 residual check at 614 is a consistency check. +- The 6×6 max-norm product inequalities that give 36 = 6·6 in θ (783). +- **The induction.** θ(R) ≤ 1, the base j ≤ NK and the definitions of C_W and C_H imply ‖H_j‖_max ≤ C_H R^j and ‖W_j‖_max ≤ C_W R^j for all j. Only θ ≤ 1/2 is checked (800). +- G(t) is invertible, with G(t)⁻¹ = Σ H_j t^j, for 0 ≤ t < 1/R. This follows from absolute convergence and the Cauchy product; there is no separate first-cell check. +- The Laurent coefficient identities for μ_k and Σ_k, the split into the i = −1, a = −1 and c = −1 pole groups, and the index bounds J+1, J+2 and J+3 in `mu_tail` and `Sig_tail` (861-873). +- The closed form T[N+1] = C(qt)^{N+1}/(1−qt), and that P2 and P3 equal the stated sums over i+l ≥ m (819-851). +- Every tail bound is monotone in t, so evaluating it at t = 2⁻¹⁴ covers all of (0, r*] (859, 875-876). +- The 2MT/4ST widening of the rotated law in `band_lo` (1205-1206), which uses that each row of LR has ℓ¹ norm at most 2. + +## Scope 3. B2: the Neumann remainder in `neumann_law` (1048-1116) + +**Where it is used.** `neumann_law(c1, c2)` gives the law on cells 2 to 8 (1260). Two other calls are not load-bearing: one at the degenerate cell [1/20, 1/20] for the CB7 cross-check (1308), and, under `--mc` only, calls at 1/20 and 1/100 for the non-certifying Monte Carlo (1333). + +**What it consumes:** +- helpers: `horner` 1001-1005, `ser_eval` 1007-1008, `mat_poly_eval` 1010-1016, `vec_poly_eval` 1018-1023, `frob_up` 1025-1026, `l2_up` 1028-1029, `mp_inverse6` 1031-1040, `gram_tail` 882-885, `ck` 101-104, and `ub`/`up`/`iv_frac` 241-255; +- CB3 data: 617-676, and `Xhi`, `Dhi`, `RWl`, `RHl`, `SQRT6` (1042-1046); +- CB4 constants: `M_GP`, `M_GY`, `M_GYV` (740-742) and `RHOi` (769). `NK` (456) fixes the tail order. + +**What it does not use.** +- `mu_tail` and `Sig_tail` (861-873). They are evaluated once at 875-876, and their stored values `MT`, `ST` are used only by `first_cell_law` (1118-1122) on the first cell. +- `polys` (1132-1146). `ceil_hi` (1156) and `band_lo` (1212) call it later on the same `t_iv`. +- `H_m` (758). It appears only in the `old_remainder` control branch (1077-1087), which runs only under that mutation and is not load-bearing. + +### 3.1 Notation +- J = 60 and NK = J+2 = 62 (455-456). ρ = `RHO` = mpf('0.11') (458). The same binary number is used for the M_* constants (through `iv.mpf(RHO)` in 712-737) and for `gram_tail`. +- z is the complex Laurent variable. On a cell, z = t = r is real, with t ∈ [c1, c2] and 0 < c1. +- V = (V_0..V_5) are the pin functionals (`PINS_EXACT`, 154-161), and Y = (q_M, q_S, a_M, a_S, b_M, b_S) are the soft functionals (`SOFT_EXACT`, 162-169). +- **zV** is the constant conditioning value `ZV_EXACT` = (6/5, 0, 0, 0, 0, −1/6) (170, 609). It means "the value of V", not z·V. +- **Gram blocks.** G(t) = Cov(V, V) is 6×6. A(t) = Cov(Y, V) is soft × pin (code `GYV`, 553, 619). GY(t) = Cov(Y, Y) (552). +- **True Laurent coefficients:** G_j (j ≥ 0), A_j (j ≥ −1), GY_j (j ≥ −2). Their enclosures for j ≤ NK are `Gj[j]`, `Aj[j]`, `A_m1`, `GYj[j]`, `GY_m1` and `GY_m2` (617-622). +- **Truncations:** G_P = Σ_{j=0}^{NK} G_j t^j, A_P = Σ_{j=−1}^{NK} A_j t^j and GY_P = Σ_{j=−2}^{NK} GY_j t^j. Then δG = G − G_P, δA = A − A_P and δGY = GY − GY_P. +- **Taylor points.** `W[j]` and `H[j]` (j ≤ NK) are the interval vectors and matrices from the recursion at 624-634. + - Fix any points Ŵ_j ∈ W[j] and Ĥ_j ∈ H[j], once for the whole cell, and set W_P = Σ Ŵ_j t^j and H_P = Σ Ĥ_j t^j. + - Nothing below needs W_P or H_P to be the Taylor polynomials of G⁻¹zV and G⁻¹; they only have to be polynomials. Being close to Taylor affects only how tight the bound is. +- w = G⁻¹zV, e_w = w − W_P and E_H = G⁻¹ − H_P. μ = A w and Σ = GY − A G⁻¹ Aᵀ (Gaussian conditioning, centred field). +- **Product polynomials:** + - X(t) = A_P W_P = Σ_{k=−1}^{2NK} X_k t^k (`X`, 641-647); + - D(t) = GY_P − A_P H_P A_Pᵀ = Σ_{k=−2}^{3NK} D_k t^k (`AH`, `D`, 650-663); + - RW(t) = zV − G_P W_P and RH(t) = I − G_P H_P, each Σ_{k=0}^{2NK} (`RW`, `RH`, 666-676). +- **Index ranges, all complete:** + - X (644): i ∈ [−1, NK] with l ∈ [0, NK]; + - AH (653): a ∈ [−1, NK] with b ∈ [0, NK]; + - D (660): c ∈ [−1, NK] with m ∈ [−1, 2NK], so every triple a ∈ [−1, NK], b ∈ [0, NK], c ∈ [−1, NK] appears; + - RW and RH (671): i ∈ [0, NK] with l ∈ [0, NK]. +- **Heavy parts:** μ_h = Σ_{k=0}^{J} X_k t^k (`mu_ser`, 648) and Σ_h = Σ_{k=0}^{J} D_k t^k (`Sig_ser`, 664). +- **Preconditioner.** B is the mpf point matrix ≈ (midpoint of the cell enclosure `GP_iv`)⁻¹ (1069-1071); under the `no_preconditioner` control, B = I. E(t) = I − B G(t). +- **Norms:** ‖·‖₂ is the spectral norm for matrices and the Euclidean norm for vectors, ‖·‖_F is the Frobenius norm, and |·|_max is the largest entry. + +### 3.2 Exact claim +For every Neumann cell and every t ∈ [c1, c2] (indeed every t ∈ `t_iv`): +- G(t) is invertible; +- |μ_a(t) − μ_h,a(t)| ≤ T_mu for all a; +- |Σ_ab(t) − Σ_h,ab(t)| ≤ T_Sig for all a, b. + +Hence μ_a(t) ∈ `mu[a]` and Σ_ab(t) ∈ `Sig[a][b]`, the values returned at 1114-1116. + +### 3.3 Hypotheses consumed from other scopes (not re-proved here) +- **H1 (scope 1, CB2).** In the exact Laurent structure, G has no negative powers, A has none below t⁻¹ and GY has none below t⁻². The dropped coefficients are exactly 0, and μ, Σ are given by the conditioning formula for these functionals. CB2s is an exact `Fraction` lemma that is run and checked (268). The passage from the functional expansions to the Gram expansions is derivation-only. +- **H2 (CB2 kernel window, 497-511, with the CB1 moment table).** The enclosures in 617-622 contain the true coefficients for every j ≤ NK. +- **H3 (scope 2, 706-742).** M_GP ≥ sup_{|z|=ρ} max_ab |G_ab(z)|, and likewise M_GYV for A and M_GY for GY. K1 is entire, so each block is meromorphic on ℂ, and its only possible pole is at 0 (removable for G by H1). So the Laurent series converge on 0 < |z| < ∞, the Cauchy estimate |c_j| ≤ M ρ^{−j} holds for every j, and the series values at real t ∈ (0, r_hi] are the Gram matrices at r = t. +- **H4.** mpmath `iv` is inclusion-monotone with outward rounding. + +### 3.4 Derivation, step by step +**S0. The cell (1058-1063).** `t_iv` = [lo.a, hi.b] ⊇ [c1, c2], with `r_lo` = lo.a ≤ c1 and `r_hi` = hi.b ≥ c2. The code checks `r_lo > 0` (1063). + +**S1. Gram tails (1064-1066; `gram_tail`, 882-885).** +- For 0 < t ≤ r_hi and x = r_hi/ρ, H3 gives |δG_ab(t)| ≤ Σ_{j≥NK+1} M_GP (t/ρ)^j ≤ M_GP x^{NK+1}/(1−x) = `tailG`. The bound holds for all t in the cell, because every term is nonnegative and increasing in t. +- `tailA` (from M_GYV) and `tailY` (from M_GY) are the same bound, applied entrywise to δA and δGY. +- **x < 1 is checked.** The geometric sum needs x < 1. `gram_tail` now checks it on the interval x = r_hi/ρ, as `ck(x.b < 1, 'CB6 Gram tail needs r < rho')` (884). The actual values are x = 0.0227, 0.0455, 0.0909, 0.182, 0.273, 0.364 and 0.4545 on cells 2-8, and 0.4545 at the 1/20 calls, so x ≤ 0.4546 on every call. + +**S2. Consistent point choice (lemma P).** +- Take the true G_j, A_j, GY_j, which lie in their enclosures by H2, and the fixed Ŵ_j, Ĥ_j. +- Every interval that `neumann_law` reads is an interval extension of the corresponding point expression. Those intervals are `X[k]`, `D[k]`, `RW[k]`, `RH[k]`, `GP_iv`, `WP`, `HP` and `AP`, and each is built from the same objects `Gj`, `Aj`, `A_m1`, `GY*`, `W`, `H`. +- So by H4, each point quantity at t lies in its interval at `t_iv`. The dependency problem only widens intervals. +- This is what lets one write exact identities between point quantities and then bound each term by an interval quantity. + +**S3. Neumann bound on G⁻¹ (1067-1076).** +- `GP_iv` = G_P(t_iv) ± tailG, so it contains G(t) by S1 and S2. +- `Emat` ⊇ {I − B G(t)}, and `rhoE` = `frob_up(Emat)` ≥ ‖E(t)‖_F ≥ ‖E(t)‖₂. +- **Checked:** `rhoE < 0.5` (1074). The lemma itself needs only < 1. +- Then BG = I − E is invertible, with ‖(BG)⁻¹‖₂ ≤ 1/(1 − ‖E‖₂). Hence G(t) is invertible on the whole cell. This is load-bearing, because the conditional law exists only then. Also G⁻¹ = (BG)⁻¹B. +- So ‖G⁻¹‖₂ ≤ ‖B‖₂/(1 − ‖E‖₂) ≤ ‖B‖_F/(1 − rhoE) ≤ `Gnorm`, where `BF` = `frob_up(Biv)`. +- Any B works: the midpoint inverse affects only tightness. + +**S4. Exact identities.** These hold for any polynomials W_P, H_P; each line follows by expanding. +- (I1) G e_w = zV − G W_P = RW − δG W_P. +- (I2) G E_H = I − G H_P = RH − δG H_P. +- (I3) μ = (A_P + δA)(W_P + e_w) = A_P W_P + δA W_P + A e_w. The last term uses the full A. +- (I4) A G⁻¹ Aᵀ = (A_P + δA)(H_P + E_H)(A_P + δA)ᵀ. Therefore Σ = D + δGY − [A_P H_P δAᵀ + δA H_P A_Pᵀ + δA H_P δAᵀ] − A E_H Aᵀ. + +**S5. Bound on e_w (1089-1093).** ‖e_w‖₂ ≤ ‖G⁻¹‖₂ (‖RW‖₂ + ‖δG‖₂ ‖W_P‖₂), using ‖δG‖₂ ≤ ‖δG‖_F ≤ 6·tailG. In code, `Ew` = Gnorm·(rWn + 6·tailG·WPn), with `rWn` = `l2_up`(RW(t_iv)) and `WPn` = `l2_up`(W_P(t_iv)). + +**S6. Bound on μ (1094-1100).** `AP` = A_P(t_iv), including `A_m1/t_iv`. `AmaxP` ≥ |A_P|_max, and `Amax` = AmaxP + tailA ≥ |A|_max. For each a: +- (i) the high orders: |Σ_{k=J+1}^{2NK} X_k t^k| ≤ `ext_mu` (Horner of `Xhi` times t^{J+1}); +- (ii) the pole: |X_{−1}/t| ≤ max|X_{−1}|/r_lo = `pole_mu`; +- (iii) |(δA W_P)_a| ≤ Σ_l |δA_al||W_P,l| ≤ 6·tailA·WPmax; +- (iv) |(A e_w)_a| ≤ ‖row_a A‖₂ ‖e_w‖₂ ≤ √6·Amax·Ew. + +By (I3), μ − μ_h is exactly the sum of the four quantities bounded in (i)-(iv). The sum of their bounds is `T_mu`. + +**S7. Bound on E_H (1101-1104).** From (I2), ‖E_H‖₂ ≤ Gnorm·(‖RH‖₂ + ‖δG‖₂ ‖H_P‖₂) ≤ Gnorm·(RHF + 6·tailG·HPF) = `EH`. This uses ‖XY‖₂ ≤ ‖X‖₂‖Y‖₂ and ‖·‖₂ ≤ ‖·‖_F. + +**S8. Bound on Σ (1105-1112).** For each a, b: +- the high orders of D: ≤ `ext_S` (orders J+1..3NK; GY_k is included there for k ≤ NK); +- the poles: |D_{−2}/t² + D_{−1}/t| ≤ `pole_S`; +- |δGY_ab| ≤ `tailY`; +- each cross term: |Σ_{k,l} δA_ak H_kl A_P,bl| ≤ tailA·AmaxP·Σ_{kl}|H_kl| ≤ tailA·AmaxP·6‖H_P‖_F, by Cauchy–Schwarz over the 36 entries. The two cross terms together give `12·tA·HPF·AmaxP`; +- |(δA H_P δAᵀ)_ab| ≤ 6·tA²·HPF; +- |(A E_H Aᵀ)_ab| = |row_aᵀ E_H row_b| ≤ (√6·Amax)²·‖E_H‖₂ = 6·Amax²·EH. + +By (I4), Σ − Σ_h is exactly the sum of these terms, and the sum of their bounds is `T_Sig`. + +**S9. Assembly (1114-1115).** μ_a ∈ Horner(`mu_ser`) ± T_mu and Σ_ab ∈ Horner(`Sig_ser`) ± T_Sig. Downstream, `band_lo` uses 2·T_mu and 4·T_Sig for the δ = q_S − q_M transform (1205-1206; `LR` at 680, `mu2_ser`/`Sig2_ser` at 681-683). That is valid because every row of `LR` has absolute sum ≤ 2, so |(LR δΣ LRᵀ)_ab| ≤ 2·2·T_Sig. + +**Poles.** +- With the exact Taylor coefficients (W_0 = G_0⁻¹zV, H_0 = G_0⁻¹, H_1), CB2s makes X_{−1} = E G_0 G_0⁻¹ zV = E zV = 0, and likewise D_{−2} = D_{−1} = 0. +- For the arbitrary points of S2 they are only tiny, not zero. That is why `pole_mu` and `pole_S` are a required part of the bound. The enclosures `X[-1]`, `D[-2]` and `D[-1]` all contain 0 (checker scratch). +- The Gram-level coefficients below the kept orders (G below t⁰, A below t⁻¹) are dropped, not bounded (H1). GY has no orders below t⁻² to drop. Their 25 stored enclosures all contain 0, with max |c| = 6.3e−100 (checker scratch). +- The run-time guard at 555-567 only requires < 1e−60. That is far too weak to protect T_mu: a 1e−60 coefficient at t⁻⁶ on cell 2 would add about 1.9e−35. So validity here rests on CB2s being exact, not on that guard. + +**Changes from the frozen remainder (`old_remainder`, 1077-1087).** +- The frozen code (ceil.txt:604 and 618-627, with the same lines in floor.txt) applies G⁻¹ − G_P⁻¹ = −G⁻¹ δG G_P⁻¹. Its comment drops the sign, which is immaterial for norms. +- But it applies the identity to the truncated Taylor polynomial H_P (`GinvPnorm`, from the `H_m` majorants), not to G_P⁻¹. For H_P the exact relation is G⁻¹ − H_P = G⁻¹(RH − δG H_P). +- The frozen bound therefore omits RH and, for the mean, the residual RW. It also omits the high orders `ext_mu`/`ext_S`: orders J+1..2NK of A_P W_P and J+1..3NK of D. + +### 3.5 Checked numerically at run time +- `r_lo > 0` (1063), x < 1 in each `gram_tail` call (884) and `rhoE < 0.5` (1074). These are the only hypotheses of the derivation that are checked at run time. The `no_preconditioner` control fails the contraction check (test:421). +- `T_mu, T_Sig < 1e−2` (1113). This is a sanity bound only; no inequality depends on it. +- **Upstream consistency checks, not proofs:** + - CB2 shows that the dropped Gram pole enclosures are < 1e−60 (555-567). This is far weaker than the scale that matters here; see Poles. + - CB3 shows that the μ/Σ pole enclosures are tiny and that A_{−1} = E·G0 (686-695). + - CB4 checks |c_j| ≤ M ρ^{−j} on the computed coefficients of orders −2..NK only (772-776). The orders above NK, which are the ones the tails need, are not checked. +- **CB7 (1307-1312)** checks that the enclosure at r = 1/20 contains a dps-100 direct law. + - The printed deviation, 1.37e−91 for μ and 2.06e−91 for Σ, is exactly that reference's own rounding error, |direct dps-100 − direct dps-200|. Against dps 200, the Neumann midpoint is off by 6.7e−101 and 1.1e−100. + - The corrections there are T_mu = 1.49e−10 and T_Sig = 1.23e−7. So CB7 cannot detect an under-estimated remainder anywhere between about 1e−91 and 1e−2. + - The heavy part's half-width (≈ 9e−98) is below the reference error, so CB7 would also reject a valid remainder below about 1e−91. It tests the heavy polynomial at about 1e−91 resolution, not the remainder. +- The `old_remainder` control is caught only by the replay digests (test:424), not inside the run. + +### 3.6 What rests only on this derivation +Everything from S1 to S9 except the checks just listed: +- identities I1-I4 and lemma P; +- the Neumann lemma and every norm inequality; +- the geometric Cauchy tail (its condition x < 1 is now checked at 884; the tail bound itself is not); +- monotonicity of the tail in t; +- the pole bound with r_lo; +- H1-H4, including the exact vanishing of the dropped Gram-level orders and the identification of the Laurent series with the Gram matrices at real t. + +## Scope 4. Envelope algebra (`ceil_hi` 1148-1198, `band_lo` 1200-1247) + +### 4.0 Claim and imported hypotheses +**Claim.** Take a cell C from 1249-1250. For every r in C, E[G_r] ≤ ⌈hi⌉ (1194, used at 1264-1266) and E[G_r] ≥ lo.a (1243, used at 1269-1271). + +**Hypotheses.** None of these is checked in this scope. +- **H_law.** The law step returns the following. + - An interval t with C ⊂ t ⊂ [0, ∞). On Neumann cells, t comes from 1060, and the check at 1063 makes t.a > 0. On the first cell, t comes from 1119, with t.a = 0. + - Boxes `mu` (6 entries) and `Sig` (6×6) with μ(r) ∈ mu and Σ(r) ∈ Sig, entrywise, for every r in C. + - Numbers Tm, Ts that bound every entry of the error of the order-J Taylor polynomials, whose coefficients are enclosed by `mu_ser` and `Sig_ser`. +- **H_iv.** mpmath `iv` is outward-rounded and inclusion-isotone. Also mp.prec = iv.prec = 336: both are set to dps 100 at 271-272, cov_exact sets them again, and 286 checks the dps. So every conversion mpf(x.a), mpf(x.b) (in `up`, `ub` and `t0_lower`) and every iv.mpf(mpf) of a constant is exact. +- **H_Φ.** `iv_Phi` (367-395) encloses Φ, and `iv_phi` (364-365) encloses φ. +- **H_def.** G_r = p·1_T with T = {q_M < 0, D_M > 0, D_S < 0} (README:12-14). This is the upstream definition and is not re-derived here. + +### 4.1 Notation +- **Variables.** X = (q_M, q_S, a_M, a_S, b_M, b_S), indices 0..5 as in `p_affine` (1133), with X ~ N(μ(r), Σ(r)). +- **Polynomials (1134-1143).** D_M = a_M q_M − r b_M², D_S = a_S q_S − r b_S², and p = −D_M D_S. u_M = ½(a_M q_M)² + ½ + r b_M², u_S is the same with the S variables, and p̄ = u_M u_S. +- **Rotated coordinates.** δ = q_S − q_M, L = `LR` (680), and LX = (q_M, δ, a_M, a_S, b_M, b_S). p̄₂ is p̄ written with q_S = q_M + δ (1144-1145), so p̄₂(LX) = p̄(X). +- **Constants (1124-1129).** These are mpf values at dps 100: I_LO = −3.9, I_HI = −0.4, η = 0.3, BB = 2, AB = 0.8, EDGE_AM = −0.6. Each event and each margin check uses the same binary number. The events convert it with iv.mpf, which is exact by H_iv; the margins convert it with `mpf_to_fraction`, which is exact. +- **Functions.** Q(t) = 1 − Φ(t), ⌈v⌉ = `up(v)` (248), and `ub(v)` = max(|v.a|, |v.b|) (241). +- **Shorthand.** TE(P; c; [α, β]) = E[P(X)·1{α ≤ X_c ≤ β}]. E2p̄ = ⌈E[p̄²]⌉ (1162, 1222) and E2p = ⌈E[p²]⌉ (1163) are upper endpoints of the Wick enclosures. + +### 4.2 Gaussian primitives +Each primitive is an exact identity, evaluated by interval extension on (t, mu, Sig). +- **G1. Wick (907-929).** E[X_i m(X)] = μ_i E[m] + Σ_j Σ_ij E[∂_j m]. This is Gaussian integration by parts, applied recursively to monomials (916-918). +- **G2. Conditional law (971-998).** X | X_c = x ~ N(μ + Σ_{·c}(x − μ_c)/Σ_cc, Σ − Σ_{·c}Σ_{c·}/Σ_cc) (974-975). The mean is stored as the affine pair [μ_i − Σ_ic μ_c/Σ_cc, Σ_ic/Σ_cc]. The G1 recursion with polynomial-in-x coefficients (985-987) gives E[P | X_c = x] = Σ_k π_k x^k. +- **G3. Truncated moments (888-905).** Write X_c = μ_c + sZ, α' = (α − μ_c)/s, β' = (β − μ_c)/s and I_j = ∫_{α'}^{β'} z^j φ. + - I_0 = Φ(β') − Φ(α') (895) and I_1 = φ(α') − φ(β') (896). + - I_j = α'^{j−1}φ(α') − β'^{j−1}φ(β') + (j−1)I_{j−2} (898). This is integration by parts with φ' = −zφ. + - E[X_c^k 1{α ≤ X_c ≤ β}] = Σ_j C(k, j) s^j μ_c^{k−j} I_j (903). +- **G4. Truncated expectation.** TE = Σ_k π_k E[X_c^k 1{…}], by the tower property (1166-1170, 1214-1217, 1228-1231). +- **Enclosure.** By H_iv, H_Φ and H_law, G1-G4 contain the true value at (r, μ(r), Σ(r)) for every r in C. + - G2-G4 are used only for c ∈ {0, 1, 4}, and each G2 use is followed by a G3 call on the same Sig[c][c]. + - Σ_cc(r) > 0 for these c follows from H_law once Sig[c][c].a > 0, and that is now checked: `trunc_mom_iv` requires `var.a > 0` (890) before its square root (891), at every truncated moment. In the recorded runs the lower endpoints are 2.0, 2.0 and at least 0.4997. + - If Sig[c][c].a ≤ 0, the divisions at 974-975 run first and do not raise. When the enclosure contains 0 they give unbounded intervals (in mpmath 1.4.1, 1/[0,1] = [1, +∞] and 1/[−1,1] = [−∞, +∞]). Either way the run then stops at the check at 890 with `H3_BAND_V2 FAIL`. +- **G5. Tail bound (397-433).** Take X ~ N(m, v) with m ∈ mu_x and 0 ≤ v ≤ v̄ = ⌈var_x⌉. + - `t0_lower` (413-415) first requires v̄ > 0 (414). It then takes t0 as the lower endpoint of an interval quotient over iv.sqrt(v̄) (415). + - Only v̄ enters, so a negative lower endpoint of var_x is harmless. Such endpoints do occur: on the first cell, Sig[2][2].a = Sig[3][3].a ≈ −1.8e−38 and Sig2[1][1].a ≈ −7.4e−38. + - The check at 432 enforces t0 > margin. Production margins are 2 (1177, 1236) and 20 (1182, 1185, 1189); margin 1 is used only in the self-test at 450. + - **'gt' (424).** t0 ≤ (e − ⌈mu_x⌉)/√v̄. Since t0 > 0, e > m and (e − m)/√v ≥ t0. So P(X > e) ≤ Q(t0). + - **'lt' (427).** The same argument with mu_x.a gives P(X < e) ≤ Q(t0). + - **'abs' (430-431).** M = ub(mu_x) ≥ |m| and t0 ≤ (e − M)/√v̄, so e > M. Then P(|X| > e) = Q((e − m)/√v) + Q((e + m)/√v) ≤ 2Q((e − |m|)/√v) ≤ 2Q(t0). This uses sides = 2. + - The function returns ⌈iv_Phi(−t0)⌉·sides (399, 402). Multiplying by 2 is exact. +- **G6. Cauchy–Schwarz.** Since p̄ ≥ 0, 0 ≤ E[p̄ 1_A] ≤ √(E[p̄²]·P(A)) ≤ √(E2p̄·P̄_A), where P̄_A is the G5 bound. The code returns [0, ⌈√·⌉] (1172-1173, 1237). + - RESp (1182) uses |E[p 1_A]| ≤ √(E[p²]·P(A)). + - The checks at 1164 and 1223 make sure the square-root argument is nonnegative. + - The check labels at 1177 and 1236 call these pieces "Chernoff". A label is printed only in the message of a failing margin check, so neither reviewed transcript contains the word. No exponential-moment bound is used, and the comment at 1150-1151 now says so. + +### 4.3 The majorant (derivation only) +- Take r ≥ 0. Then |D_M| ≤ |a_M q_M| + r b_M² ≤ u_M, using |x| ≤ (x² + 1)/2. Likewise |D_S| ≤ u_S. +- Hence |p| ≤ p̄ and p̄ ≥ ¼. So for every event A, ±E[p 1_A] ≤ E[p̄ 1_A]; if A ⊂ ∪A_i, then E[p̄ 1_A] ≤ Σ_i E[p̄ 1_{A_i}]. +- P and P̄ are built with the interval r = t (1156, 1212). The condition r ≥ 0 comes from t ⊂ [0, ∞) (H_law). + +### 4.4 Ceiling (`ceil_hi`) +- **U1.** {q_M < 0} \ T = A_M ∪ A_S, with A_M = {q_M < 0, D_M ≤ 0} and A_S = {q_M < 0, D_S ≥ 0}. So E[G_r] = E[p 1{q_M < 0}] − E[p 1_{A_M ∪ A_S}] ≤ E[p 1{q_M < 0}] + E[p̄ 1_{A_M}] + E[p̄ 1_{A_S}]. +- **U2.** A_M ⊂ {a_M > EDGE_AM} ∪ {|b_M| > BB} ∪ {I_HI ≤ q_M ≤ 0}. Suppose a point of A_M lies in none of the three sets: then a_M ≤ EDGE_AM, |b_M| ≤ BB and q_M < I_HI. In that case a_M q_M > (−EDGE_AM)(−I_HI) > ⌈t⌉·BB² ≥ r b_M², so D_M > 0, a contradiction. The middle inequality is the exact Fraction check at 1155, made per cell. +- **U3.** {D_S ≥ 0} ⊂ {q_S ≥ 0} ∪ {a_S ≤ 0}: if q_S < 0 < a_S, then D_S ≤ a_S q_S < 0. No run-time check is needed. +- **U4.** {q_M < 0} splits into the windows [I_LO, I_HI], [I_HI, 0] and [−100, I_LO], plus the residue {q_M < −100}. p is signed, so this has to be a partition, not just a cover. The shared endpoints are null sets, because Σ_00 > 0. + +| term | line | upper endpoint bounds | method | +|---|---|---|---| +| Wp1, Wp2, Wp3 | 1179-1181 | the three TE(p; 0; window) | G4, two-sided | +| RESp | 1182 | E[p 1{q_M < −100}] | G6 with E2p; 'lt', e = −100, margin 20 | +| BaM | 1183 | E[p̄ 1{a_M > −0.6}] | G6; 'gt', margin 2 (1175-1177) | +| BbM | 1184-1186 | E[p̄ 1{\|b_M\| > 2}] | TE on [−100, −2] and [2, 100]; G6 'abs' at e = 100, margin 20 (two tails) | +| BqM | 1187 | E[p̄ 1{I_HI ≤ q_M ≤ 0}] | G4 | +| BqS | 1188-1189 | E[p̄ 1{q_S ≥ 0}] | TE on [0, 100]; G6 'gt' at e = 100, margin 20 | +| BaS | 1190 | E[p̄ 1{a_S ≤ 0}] | G6; 'lt', e = 0, margin 2. t0 > 2 forces m > 0, so P(a_S < 0) = P(a_S ≤ 0) | + +hi is the sum of these nine terms (1194), and U1-U4 with 4.3 give E[G_r] ≤ ⌈hi⌉. The cell decimal is rounded up exactly at 1266, after the cap check at 1265. The uniform bound is the maximum over cells (1274-1277). + +### 4.5 Floor (`band_lo`) +- **F0. Rotated law (1205-1206).** Lμ and LΣLᵀ are built from series-level coefficients (681-683), with tails 2Tm and 4Ts, rounded up. Each row of L has ℓ¹ norm ≤ 2, so H_law bounds the errors by 2Tm for the means and by 2·2·Ts = 4Ts for the covariances. This gives δ ~ N(μ_1 − μ_0, Σ_00 − 2Σ_01 + Σ_11), enclosed by mu2[1] and Sig2[1][1]. +- **F1. The box.** Q = {I_LO ≤ q_M ≤ I_HI}, and B_cert = Q ∩ {|δ| ≤ η} ∩ {a_M ≤ −AB} ∩ {a_S ≥ AB} ∩ {|b_M| ≤ BB} (1201). +- **F2. B_cert ⊂ T** (exact checks at 1209-1211). On B_cert: + - q_M ≤ I_HI < 0; + - D_M ≥ AB(−I_HI) − ⌈t⌉BB² > 1/20 (1209); + - q_S = q_M + δ ≤ I_HI + η < 0. Since a_S ≥ AB > 0 and q_S ≤ 0, D_S ≤ a_S q_S ≤ AB·q_S ≤ AB(I_HI + η) < −1/20 (1210). + + So the fact q_S < 0 is load-bearing for D_S < 0. The separate check at 1211 is redundant: 1210 with AB > 0 already forces I_HI + η < 0. +- **F3.** p > 0 on T, so p 1_T ≥ p 1_{B_cert}. Also Q \ B_cert ⊂ {|δ| > η} ∪ {|b_M| > 2} ∪ {a_M > −0.8} ∪ {a_S < 0.8}. With 4.3 this gives + + E[G_r] ≥ E[p 1_Q] − E[p̄ 1{|δ| > η}] − E[p̄ 1{|b_M| > 2}] − E[p̄ 1{a_M > −0.8}] − E[p̄ 1{a_S < 0.8}]. + +| term | line | encloses / bounds above | method | +|---|---|---|---| +| EqI | 1214-1217 | E[p 1_Q] | G4, two-sided | +| Edelta | 1239 | E[p̄ 1{\|δ\| > η}] | forced G6; 'abs' on (mu2[1], Sig2[1][1]) at e = η, margin 2, two tails (B3). E2p̄ comes from the unrotated law, which is valid because p̄₂(LX) = p̄(X) | +| Ebeta | 1240 | E[p̄ 1{\|b_M\| > 2}] | if Sig[4][4].a > 1e−4 (1227): Ebar − TE(p̄; 4; [−2, 2]) (1231); otherwise G6 'abs' at e = 2. The first branch is taken on every recorded cell (Sig[4][4].a ≥ 0.4997) | +| EaM | 1241 | E[p̄ 1{a_M > −0.8}] | forced G6; 'gt' at e = −AB | +| EaS | 1242 | E[p̄ 1{a_S < 0.8}] | forced G6; 'lt' at e = AB | + +- lo = EqI − Edelta − Ebeta − EaM − EaS (1243). Interval subtraction gives lo.a ≤ EqI.a − Σ⌈piece⌉, which is at most the right side of F3. +- The cell decimal is rounded down exactly at 1271, after the c_Z check at 1270. +- The event map at 1235 reads only b for the 'abs' kind, so it is correct only when a = −b. The check at 1234 now enforces a = −b on the tail branch whenever no side is given. Both call sites (1239, 1240) satisfy it. +- PBAR2 and Ebar2 are passed to Edelta but not used there (`force_chern`). They serve only the consistency check at 1221. + +### 4.6 Frozen programs versus v2 +**Unchanged.** This was checked by reading the code, plus a normalized diff of the helpers. +- The events, the constants, the ±100 windows, and the margins 2 and 20. +- The nine ceiling summands (ceil.txt:693-712 against 1179-1194). +- The five floor terms (floor.txt:699-702 and 704 against 1239-1243) and their branch choices, including the 1e−4 switch (floor.txt:681 against 1227). +- p, p̄ and p̄₂ (ceil.txt:653-662; floor.txt:651-660 and 668-673). +- G1-G4. The frozen ceil.txt:472-562 and floor.txt:472-562 are byte-identical to each other. Apart from formatting and the guard at 890, the only change is that `trunc_mom_iv` now takes a and b as intervals, and its callers wrap them with `iv.mpf`. + +**Changed.** +- **B3.** The frozen floor's 'abs' branch used one tail, P = iv_Phi(−t0).b (floor.txt:692-697); v2 uses two (431). The frozen ceiling's far tail at |b_M| > 100 already used two (ceil.txt:705). Every tail now goes through one function, `event_prob_upper` (417-433). +- **Rounding of t0 and of the square-root pieces.** The frozen programs used round-to-nearest mp (ceil.txt:685-691 and 696-710; floor.txt:687-697). v2 uses an interval lower bound for t0 (415) and an interval upper bound for the square root (1173, 1237). +- **Margin checks.** The frozen checks were done in mpf: the ceiling cover used the constant 0.05 (ceil.txt:652), and the floor used the mpf decimal cell end (floor.txt:645-649). v2 checks exactly in Fraction, with ⌈t⌉ taken per cell. +- **Added:** the CB5 tail-helper and event-bound self-tests (404-411, 435-452); the second-moment checks (1164, 1223); the fail-closed guards (414, 890, 1234); and the mutation hooks `tail_helper_one_side` (400-401), `one_tail` (431) and `drop_fail_piece` (1191-1192). +- **Removed:** the ceiling nesting check (ceil.txt:724), listed at README:84; the floor check "piece enclosures sane" (floor.txt:703), now listed at README:45; and the dead non-forced path in the frozen ceiling's upper_half (ceil.txt:679-683), both of whose calls forced the tail path. +- **Indirect (H_Φ scope).** These changes reach every tail and every truncated moment. O2 removes float() in the Gordon–Mills branches (378 and 382, against frozen lines 451 and 456). `iv_phi` uses iv.exp (365, against the frozen iv.e ** at line 437). The `iv_Phi` series tail is outward (394, against frozen line 469). + +**Cross-check against the transcripts.** On cells where the law did not change at the printed precision, frozen floor − v2 floor = (1 − 1/√2)·Edelta_v2: + +| cell | frozen − v2 floor | ratio to (1 − 1/√2)·Edelta_v2 | +|---|---|---| +| (1/100, 1/50] | 1.26e−12 | 1.0004 | +| (1/50, 3/100] | 9.43e−6 | 0.9995 | +| (3/100, 1/25] | 2.540e−3 | 1.0002 | + +On (1/25, 1/20] the ratio is 0.578, because B2 changes the law there. The ceilings on (1/400, 1/200] through (1/50, 3/100] equal the frozen ones up to the last digit, which v2 now rounds outward. + +### 4.7 What the program checks at run time +- The cover margin (1155) and the type margins (1209-1211; 1211 is implied by 1210). +- t0 > margin (432) at every tail evaluation: + - per ceiling cell: RESp 1182, BaM 1183, the BbM far tail 1185, the BqS far tail 1189, and BaS 1190; + - per floor cell: Edelta 1239, EaM 1241 and EaS 1242; + - Ebeta 1240 only when Sig[4][4].a ≤ 1e−4, which never happens in the recorded runs. +- v̄ > 0 before every tail's square root (414); Sig[c][c].a > 0 before every truncated moment (890); a = −b on the two-sided tail branch of `piece_outside` (1234). +- The second moments are ≥ 0 (1164, 1223), and the two E[p̄] enclosures intersect (1221). +- The tail and event self-tests against Gordon's lower bound (406-410, 442-451). These test a necessary condition on point laws only. +- Consumption: U < 4 (1265) and lo > c_Z (1270). These two are not validity checks. + +### 4.8 What rests only on the derivation +1. H_def: E[G_r] = E[p 1_T] with T = {q_M < 0, D_M > 0, D_S < 0} (README:12-14). +2. 4.3: |p| ≤ p̄ for r ≥ 0, hence ±E[p 1_A] ≤ E[p̄ 1_A] and the union bound. The code only builds the polynomials (1132-1146). +3. p > 0 on T, used in the floor for p 1_T ≥ p 1_{B_cert}. +4. F2: the step from the exact checks at 1209-1210 to B_cert ⊂ T. It includes q_S < 0, which is needed for a_S q_S ≤ AB q_S. +5. U2: the step from 1155 to A_M ⊂ {a_M > −0.6} ∪ {|b_M| > 2} ∪ {I_HI ≤ q_M ≤ 0}. +6. U3: {D_S ≥ 0} ⊂ {q_S ≥ 0} ∪ {a_S ≤ 0}. +7. U1, U4 and F3: the decomposition identities, the completeness of the nine ceiling summands (1194) and of the five floor terms (1243), and the partition of {q_M < 0} into the windows plus the residue. This item also covers the null-set conventions: shared endpoints, a_S < 0 versus a_S ≤ 0, and |b_M| = 2 in Ebar − TE. +8. G6: Cauchy–Schwarz (1172-1173, 1182, 1237). +9. G5: Gaussian tail monotonicity for 'gt' and 'lt', and P(|X| > e) ≤ 2Q((e − |m|)/√v) for 'abs' (417-433). +10. G1-G4: Isserlis, the conditional-Gaussian formulas, the truncated-moment recurrence with the binomial expansion, and the tower property (916-918, 974-975, 895-903, 1170, 1217, 1231). +11. F0: the rotated-law tails 2Tm and 4Ts, from the ℓ¹ row norms of LR (680, 1205-1206). Also that E2p̄ of the unrotated law serves for Edelta, because p̄₂(LX) = p̄(X). +12. The enclosure principle: evaluating G1-G4 and the polynomials on (t, mu, Sig) encloses the value at the true (r, μ(r), Σ(r)). This trusts H_iv, H_Φ and H_law. + +## Same-family pre-checks (credit 0) + +Same-provider agents ran these checks: the section drafters and one adversarial checker per section. They ran on the earlier program (`99eceb98…`), whose certified decimals and transcripts the final program reproduces byte for byte. The scripts are scratch files, not part of this PR. The checks earn no independence credit. They are listed so a reviewer can see what has been tried numerically, and where each check stops. + +### Scope 1 (CB2s) +- **Exact recomputation (drafter).** An independent re-expansion equals `laurent_expand` on all 12 functionals, and E·V^{(0)} = Y^{(−1)} holds by direct substitution. The program's own lines 460-548 were run on an exact Fraction shim, with planar moments and with random rational moments. They match the formal B-Gram coefficient by coefficient on orders −8..62, with 0 mismatches; the only nonzero negative orders are GYV −1 and GY −2, −1. X_{−1}, D_{−2} and D_{−1} are exactly 0, and A_{−1} = E·G_0 and GY_{−2} = E·G_0·Eᵀ hold exactly. +- **Independent exact recomputation (checker; does not import the program).** It uses hand-transcribed tables, its own expansion and its own B. It recomputes (a)-(d), E and the window margin of 6 (attained by V5×V5, V5×Y2, V5×Y3 and Y2×Y3), and gets the three pole coefficients exactly 0 for planar and for random rational moments. + - Non-vacuity: zV3 = 1 gives X_{−1} = (0, 0, 1, 1, 0, 0); Y0 = M(0,2)/z gives X_{−1,0} = −6/5 and |D_{−2}| = 2; `wrong_window` keeps the poles at 0 and moves μ_aM(0) to −1/2. +- **Interval series against formal B (checker; unmodified source with a hook after CB3).** All 7668 interval Gram coefficients (orders −8..62, 108 entries) match independent formal-B enclosures built from the program's A2 to 2.1e−99, with no disjoint pair. Σ_{p≤NK} c_p z^p matches the direct covariance `ce.cov(mode='iv')` to 5.3e−86 at z = 1/128 and to 8.4e−88 at z = 1/64. That interval path omits a denominator term (review 5607640770) of relative size about 4e−205 at the library defaults, far below these agreements. +- **First cell (drafter).** The first-cell heavy polynomials were compared with the program's direct point law (dps 100) at t = 1/16384, 1/32768, 1e−6 and 1e−9. The deviations are 5.33e−74, 2.03e−72, 5.09e−63 and 1.18e−45 for μ, and 8.0e−74, 3.04e−72, 7.64e−63 and 1.77e−45 for Σ. They grow as t → 0 because the direct law loses precision through the r⁻⁶ cancellation, and all lie below MT and ST. **Limit:** this cannot see a pole coefficient of order 1e−99. +- **Other checks.** The frozen ceiling and floor tables equal v2's. mpmath iv behaviour on zero pivots was probed (see 1.3 (e)). The pole magnitudes quoted in 1.6-1.7 come from the unmodified source run in memory with a reporting hook. +- **Final program (this guide).** `structure` prints PASS with the E of 1.2. `structure --mutation wrong_window` prints `MUTATED STRUCTURE RUN (wrong_window) … NOT A CERTIFICATE` and exits 3. `structure --mutation structure_window` fails at "CB2s (d) E zV = 0 required, row 2" and exits 1. + +### Scope 2 (CB4) +- **θ.** θ was recomputed exactly in `Fraction` from the program's mpf inputs: 0.49437044607085834, with up(θ(R)) ≥ it. By exact bisection, the smallest R with θ ≤ 1 is 1296.833 and with θ ≤ 1/2 it is 2592.833. The 2% grid returns R = 2622.3491580327, and R·r* = 0.16005. +- **Recursion past NK.** The majorant recursion stays at most 2.7e−146·C R^j for H and for W; the checker ran j ≤ 900 and the drafter j ≤ 1600. +- **MT and ST.** Recomputed by plain convolutions at dps 120, with K = 700 terms and exact 1/ρ, they agree with the program to relative 3.3e−101 and 3.5e−101. + - MT = 1.03890e−44, of which the A_{−1} pole term is 1.03884e−44. That term equals 18·[W_m[62] t^61 + C_W (Rt)^63/((1−Rt)t)], with the second part 5.77e−46. The main part is 6.3e−49. + - ST = 1.83842e−38, of which the a = c = −1 term is 1.83820e−38. The one-pole groups contribute 2.2e−42, the product group 6.8e−47 and the GY part 2.6e−203. +- **A_n(ρ).** It exceeds a direct dps-130 sum over |j| ≤ 400 with the full Z by 1.5e−99 (n = 0), 1.9e−99 (n = 2), 2.6e−97 (n = 4) and 2.2e−91 (n = 8). +- **Sampled sup against M.** GP 3.02 against 2.96e6, GY 255.5 against 343.7, and GYV 27.5 against 1.95e4. +- **Cauchy margins on the computed orders.** max_j ‖c_j‖_max ρ^j / M is 1.0e−6 for GP, 1.5e−4 for GYV and 0.022 for GY (at j = 0); 1.4e−3 for A_{−1}, 0.026 for GY_{−1} and 0.72 for GY_{−2}. So an M_GY understated by more than a factor of about 1.39 would fail the run-time check, while the GP and GYV parts would catch only errors of 3 to 6 orders of magnitude. +- **H2, two ways.** + - The drafter widened the window to J+14−n. Every coefficient through NK, and R, C, MT and ST, stayed identical. + - The checker evaluated every GP, GY and GYV entry directly from the spectral kernel (|j| ≤ 300, dps 170) at 256 points on |z| = 3, and took a DFT. This route has no kernel truncation. Every coefficient at orders 0..62, and A_{−1}, GY_{−1} and GY_{−2}, lies inside the program's enclosure: max |DFT − mid| is 4.3e−100 and max deviation/radius is 0.79. A truncated window would show errors near |c_62|, about 1e−43. + - Three exactly-zero pole entries show about 1e−130 in the DFT. That is aliasing from order 254/255: −3·3^256/(2^127·127!) = −8.13584919065e−130 matches the GY(0,1) z⁻² value digit for digit. +- **True tail.** Σ_{k=61}^{124} max|X[k]| r*^k ≈ 1.2e−287, far below MT. The computed |H_j|^{1/j} is about 0.40 and |W_j|^{1/j} about 0.34, so R = 2622 is a very loose but valid growth bound. +- **Limit.** All of these use computed orders or sampled points. None proves the Cauchy estimate beyond NK or the induction step. + +### Scope 3 (B2) +- **Drafter.** An instrumented copy reproduced `results/ceil.transcript.txt` byte for byte (sha256 `d0def90d…`). It compared each bound with a dps-200 direct law at c1, at the midpoint and at c2 of every cell. +- **Checker, at the same 21 points.** The dps-200 interval direct Grams G(t), A(t) and GY(t) intersect the program's point enclosures (Horner at [t] ± tail). They come from the library's interval path, which omits a denominator term of relative size about 4e−205 (review 5607640770), below dps-200 resolution and far below these comparisons. The dps-200 μ and Σ lie inside both the point-cell and the cell enclosures, and ‖G(t)⁻¹‖₂ ≤ 25.56 < Gnorm. + +| cell | x = r_hi/ρ | rhoE | Gnorm | tailG | Ew | EH | T_mu | T_Sig | +|---|---|---|---|---|---|---|---|---| +| (1/16384,1/400] | 0.0227 | 0.0126 | 26.0 | 8.8e−98 | 6.3e−95 | 3.5e−94 | 7.5e−90 | 5.1e−84 | +| (1/400,1/200] | 0.0455 | 0.0129 | 26.0 | 8.3e−79 | 5.9e−76 | 3.3e−75 | 1.7e−72 | 2.9e−68 | +| (1/200,1/100] | 0.0909 | 0.0258 | 26.4 | 8.0e−60 | 5.8e−57 | 3.3e−56 | 8.5e−54 | 7.0e−50 | +| (1/100,1/50] | 0.182 | 0.0516 | 27.1 | 8.2e−41 | 6.1e−38 | 3.4e−37 | 4.5e−35 | 1.9e−31 | +| (1/50,3/100] | 0.273 | 0.0516 | 27.1 | 1.2e−29 | 8.5e−27 | 4.8e−26 | 3.1e−24 | 6.5e−21 | +| (3/100,1/25] | 0.364 | 0.0517 | 27.1 | 9.8e−22 | 7.2e−19 | 4.1e−18 | 1.8e−16 | 2.4e−13 | +| (1/25,1/20] | 0.4545 | 0.0518 | 27.1 | 1.45e−15 | 1.1e−12 | 6.1e−12 | 2.0e−10 | 2.0e−7 | + +- **Actual errors.** At all 21 points, |μ − μ_h| ≤ 6.9e−101 and |Σ − Σ_h| ≤ 1.2e−100; also ‖e_w‖₂ ≤ 3.5e−100 and ‖E_H‖₂ ≤ 2.1e−99. These are at the rounding level of the dps-100 coefficients. +- **Dominant terms.** T_mu ≈ √6·Amax·Ew and T_Sig ≈ 6·Amax²·EH. Ew and EH are driven by their 6·tailG terms, while rWn ≈ 1.5e−98 to 1.9e−98 and RHF ≈ 9.1e−98 to 9.5e−98 are set by interval width. So the printed remainder is effectively the Cauchy tail of the pin Gram, times Gnorm and |A| or |A|². A carries the 1/t pole: Amax = 4.9e4 on cell 2 and 75 on cell 8. On cell 2, pole_mu = 5.6e−95 and pole_S = 2.2e−90. +- **Why the T_Sig scale matters.** On the last cell T_Sig = 2e−7, which is large enough to move the printed 15-digit decimals there. +- **Dropped coefficients.** All 25 stored dropped Gram-level enclosures contain 0 (max |c| = 6.29e−100), and so do X[−1], D[−2] and D[−1]. +- **CB7.** The printed deviation (1.37e−91 for μ, 2.06e−91 for Σ) equals |direct dps-100 − direct dps-200|. The Neumann midpoint is within 6.74e−101 (μ) and 1.14e−100 (Σ) of the dps-200 law, and the μ enclosure's half-width without T_mu is 9.03e−98. +- **x > 1, on the earlier program, before the guard at 884 existed.** On the hypothetical cells [1/20, 3/25], [1/20, 1/5] and [1/10, 1/9], plain Python stopped with mpmath's ordering AssertionError when building `GP_iv` (now line 1068). Under -O, every call failed the contraction check (now 1074), because min_{x>1} x^63/(x−1) = 169.89 makes |tailG| ≥ 5.0e8. So that path was already fail-closed; the guard now names it. +- **README:92** reports a dps-200 comparison "at 32 points". That figure is from the earlier same-family adversarial pass during authoring, on the program as it then stood (certified decimals as now). Its skeptic recomputed the direct law with `cov_exact` at dps 200 at four points per cell, c1 + (c2 − c1)·{0.001, 0.37, 0.999, 1}, over all eight cells including the first. All 32 points lay inside their cell's enclosure, with |direct − heavy| at most 6.9e−101 for μ and 1.14e−100 for Σ. The checker above used 21 points and found the same. +- **Limit.** The actual errors (≈ 1e−100) sit at rounding level, far below the corrections on every cell. No numeric test here can show that the remainder bound is valid, only that it is not violated at the sampled points. + +### Scope 4 (envelope) +- **Replays.** `ceil --verify` and `floor --verify` both PASS. `floor --mutation one_tail` fails at "CB5 event bound P(|X| > t) self-test at 2.5", with exit 1. +- **Primitives.** `trunc_mom_iv` agrees with `mp.quad` for k ≤ 8 on three windows, with relative error < 1e−25. On every cell's actual box, `event_prob_upper` dominates the exact normal-CDF supremum over the box corners at every site with t0 ≤ 60. For Edelta on (1/400, 1/200] the ratio is about 2.001, the two-tail factor. +- **Polynomials.** P, PBAR and PBAR2 (1132-1146) enclose the closed-form p and p̄ at 20 random points, for each of r = 0, 0.0371 and 0.05. +- **G1-G4 against independent quadrature.** At each cell's midpoint point law, the program's values agree with numpy Gauss–Hermite (5-D) × Gauss–Legendre quadrature to a worst relative difference of 2.3e−15. The quantities compared were E[p], E[p̄], E[p̄²], E[p²], Wp1-3, BqM, the BqS window, the BbM windows and TE(p̄; 4; [−2, 2]). +- **Point law at r = 1/20** (the CB7 direct law), with the helpers extracted verbatim. + - The envelope gives hi = 3.66697 and lo = 2.48761. Monte Carlo with 4e6 samples gives E[G] = 3.22672 ± 0.00224, inside that range. + - BqS is 0.365102 against 0.366134 ± 0.0012, and Wp1 is 2.65847 against 2.66014 ± 0.0017. + - Edelta is valid but loose: the Cauchy–Schwarz bound is 0.120, while Monte Carlo gives E[p̄ 1{|δ| > η}] ≈ 1.9e−4. + - No sample violated p̄ ≥ |p|, B_cert ⊂ T, U2 or U3. +- **Per-cell Monte Carlo.** At each cell's midpoint law, 1e6 samples give E[G] between 3.223 and 3.235 (± 0.0045). Every value lies inside that cell's [floor, ceiling]. +- **Unbounded intervals (this guide, mpmath 1.4.1).** With a variance enclosure [−1e−38, 2], [0, 2] or [0, 0], the conditional-law divisions and the products after them return unbounded intervals without raising. So the guard at 890 is reached. +- **Limit.** Quadrature and Monte Carlo use midpoint point laws. They test the primitives and the envelope numbers, not the set inclusions for every r in a cell. + +### Confirmed lower-severity findings and their disposition +Every confirmed finding has severity minor or doc. The envelope checker reported no inequality in the wrong direction, no missing term, no norm mismatch and no rounding that fails to be outward. The CB4 checker reported no blocking or major defect. + +**Fixed in the final program:** +1. `gram_tail` had no x < 1 guard (scopes 2 and 3; minor). It now has `ck(x.b < 1, 'CB6 Gram tail needs r < rho')` (884). +2. A negative variance enclosure reached `iv.sqrt` and raised ComplexResult: a negative lower endpoint in `trunc_mom_iv`, or a negative upper endpoint in `t0_lower`. `main` does not catch that exception (1491), so the run ended without the documented FAIL line (scope 4; minor). The run still failed closed: exit 1, no PASS, no transcript, no record. Both square roots are now preceded by named checks, at 414 and 890. +3. The 'abs' event map in `piece_outside` used only b, assuming a = −b (scope 4; minor). The check at 1234 now enforces it. +4. `structure --mutation` printed `STRUCTURE PASS` and exited 0 for a mutation the lemma cannot see (scope 1; doc). It now prints `H3_BAND_V2 MUTATED STRUCTURE RUN (NAME): …; NOT A CERTIFICATE` and exits 3 (1451-1454; test:91-98). +5. The `ceil_hi` comment pointed to a direction analysis in README.md, which has none (scope 4; doc). The comment at 1149 now points to section 3 of the frozen H3_BAND_CEIL.md (its lines 59-87). +6. The "Chernoff" labels hid that these pieces are Cauchy–Schwarz combined with an exact Gaussian tail (scope 4; doc). The comment at 1150-1151 now says so. +7. The `band_lo` comment used an undefined I′ (scope 4; doc). The comment at 1203 now writes [I_LO, I_HI]. + +**Documentation corrected (README):** +- README:45 now lists the dropped floor check (floor.txt:703) and the added second-moment checks. +- README:47 calls the B1 guard a regression guard. +- README:48 says that θ(R) ≤ 1/2 and the 0.2 are margins inherited from the frozen program, and states the actual convergence conditions. It also says the first cell drops the pole coefficients on CB2s alone. +- README:63 states the `structure` behaviour under a mutation, and README:64 lists the fail-closed guards. +- README:69 calls E rational, 0/1 for the tables of record. The checker had judged the earlier "0/1 matrix" wording a nit, not a defect. +- README:72 calls the numeric pole checks consistency checks only (orders from −8, threshold 1e−60). +- README:76 restricts the Cauchy sanity check to orders −2..J+2, with the higher orders resting on the derivation. +- README:92 now calls CB7's 2e−91 the rounding error of its dps-100 reference, not the true truncation error. + +**Left as is, with the reason:** +- **The 1e−60 threshold of the CB2 pole-pattern check (555-567)**, and 1e−55 and 1e−50 at 687 and 690. These are consistency checks, and validity rests on the exact CB2s. A sharper check (0 ∈ c, or |c| < 1e−95) would also pass, but it would not change what is proved. +- **The check cannot see orders below LMIN = −8.** The tables have none, and the exact lemma runs first. +- **The B1 dominance guard (767) is identically true in unmutated code.** It is a regression guard, and it catches `signed_maxima`. +- **The Cauchy sanity check (776) is loose for GP and GYV and says nothing beyond NK.** It is a sanity check; the higher orders rest on the derivation, as 2.6 says. +- **θ(R) ≤ 1/2 and R·r* < 0.2** are margins inherited from the frozen program. Validity needs only θ(R) ≤ 1, R·r* < 1 and QRHO·r* < 1. +- **No run-time tie between the kernel window and the tables.** The DFT check confirms the window for the tables of record, and `structure_soft_pole2` is stopped by CB2s (b). +- **The first cell drops the pole coefficients with no numeric fallback.** This is not a gap: the cancellations follow from (c) and (d) by the algebra in 2.3, and README:48 now says so. +- **The check at 1211 is implied by 1210.** It is harmless. The property it states, q_S < 0 on B_cert, is load-bearing and is proved by 1210. +- **The word "Chernoff" stays** in the check labels at 1177 and 1236, and in README:78. The labels are printed only in the message of a failing margin check, so neither reviewed transcript contains the word, and the comment at 1150-1151 says what the pieces are. + +**Rejected as not a defect:** the suggestion that q_S < 0 is not needed for the README definition of T. It is needed: the step a_S q_S ≤ AB q_S uses q_S ≤ 0, since a_S ≥ AB is unbounded above. + +## What a reviewer still has to check + +These are the inequalities and identities that no run checks. Each item points to where this guide derives it. + +**Scope 1 (CB2s)** +- [ ] The Gram Laurent coefficients equal B on the formal coefficients (1.4 Step 1). This covers K1 being entire, the DISP sign (541), the (−1)^q in the (M,S) case, and the S-S cancellation Σ_{q+q'=N} C(N,q')(−1)^{q'} = 0 (500-548). +- [ ] (a)-(d) plus (e) imply (i) and (ii), through Ỹ = Y − E·V/z, including the expansions of GY and A·H·Aᵀ (1.4 Steps 2-3). +- [ ] `solve_in_span` returns an E that satisfies (c) whenever it does not return None (186-209). +- [ ] The maximum cross-station weight n1 − k_a − k_b is 6, so the window at 502 is complete through NK (1.5). +- [ ] The H and W recursions (624-634) and θ's sum from i = 1 are the true inverse, given (a). +- [ ] The true G_0 lies in the box G0, so 612-614 proves (e); this also needs CB1. +- [ ] μ and Σ are the conditional law given V = zV; and, outside this scope, the tables are the upstream functionals. + +**Scope 2 (CB4)** +- [ ] |K1^{(n)}(w)| ≤ A_n(ρ) on |w| ≤ ρ, with the lattice tail ≤ 2u_{j0}/(1−q), and the lower endpoint of Z1_iv at most the exact truncated sum, which is at most Z (706-724). +- [ ] The sup bound on |z| = ρ (728-737), and Cauchy's estimate for every j, including j > NK in majA and majGY and i > J in θ. +- [ ] ‖PQ‖_max ≤ 6‖P‖_max‖Q‖_max; the induction h_j ≤ C_H R^j θ(R) and w_j ≤ C_W R^j θ(R); θ(R) ≤ 1 suffices (778-783, 800). +- [ ] G(t) is invertible on [0, 1/R), with inverse Σ H_j t^j. +- [ ] The coefficient identities for μ_k and Σ_k, the four Σ-tail groups, and the index bounds J+1, J+2 and J+3 (861-873). +- [ ] The closed form of the `Maj` suffix sums and the P2/P3 sums (819-851), and monotonicity in t, so that evaluation at 2⁻¹⁴ suffices (859, 875-876). +- [ ] The 2MT/4ST widening (1205-1206). + +**Scope 3 (B2)** +- [ ] The Gram tail M x^{NK+1}/(1−x) for x < 1, monotone in t (S1); and H3, which identifies the Laurent series with the Gram matrices at real t. +- [ ] Lemma P: the point quantities lie in the interval quantities for any fixed Ŵ_j, Ĥ_j (S2). +- [ ] The Neumann lemma and ‖G⁻¹‖₂ ≤ ‖B‖_F/(1 − rhoE) (S3). +- [ ] Identities I1-I4 (S4). +- [ ] The norm inequalities of S5-S8, including the constants 6, √6 and 12 and the Frobenius-to-spectral steps. +- [ ] The pole bounds with r_lo, and that the dropped Gram-level orders are exactly zero (via scope 1). +- [ ] The 2T_mu/4T_Sig widening in `band_lo` (S9). + +**Scope 4 (envelope)** +- [ ] H_def, upstream: E[G_r] = E[p 1_T] (README:12-14). +- [ ] |p| ≤ p̄ for r ≥ 0 and the union bound (4.3); p > 0 on T. +- [ ] U1-U4: the ceiling decomposition, the cover U2 from 1155, the inclusion U3, and the partition with its null-set conventions. +- [ ] F2: B_cert ⊂ T from 1209-1210, using q_S < 0. F3: the floor decomposition. +- [ ] G1-G4: Isserlis, conditioning, the truncated-moment recurrence and the tower property. +- [ ] G5: the 'gt', 'lt' and 'abs' tail bounds, including the factor 2 for 'abs'. G6: Cauchy–Schwarz. +- [ ] F0: the rotated tails 2Tm and 4Ts, and p̄₂(LX) = p̄(X). +- [ ] The enclosure principle under H_iv, H_Φ and H_law. diff --git a/experiments/h3_band_v2/h2/MANIFEST.sha256 b/experiments/h3_band_v2/h2/MANIFEST.sha256 new file mode 100644 index 000000000..a2b154375 --- /dev/null +++ b/experiments/h3_band_v2/h2/MANIFEST.sha256 @@ -0,0 +1 @@ +f08c1c5f653f2ffd2e85d39e1112f80d15db04d15f932e5ee42db2ca3553e783 cov_exact.py diff --git a/experiments/h3_band_v2/h2/cov_exact.py b/experiments/h3_band_v2/h2/cov_exact.py new file mode 100644 index 000000000..a160edb4c --- /dev/null +++ b/experiments/h3_band_v2/h2/cov_exact.py @@ -0,0 +1,531 @@ +#!/usr/bin/env python3 +""" +cov_exact.py -- H2 foundations LIBRARY 1: the exact-periodic covariance library +for the 2D upper-theorem campaign (work-order items 176-197). + +THE LAW. The exact normalized periodized side-24 2D Bargmann-Fock field: +spectral lattice (pi/12)Z^2, masses w(k) = exp(-|k|^2/2), Z-normalized, + + K(x, y) = (1/Z^2) sum_{k in (pi/12)Z^2} w(k) e^{i k.(x-y)}, + Z = Z_1 = sum_{j in (pi/12)Z} e^{-j^2/2} (1D normalization). + +The kernel factorizes, K(x,y) = K1(x1-y1) K1(x2-y2), with the 1D periodized +kernel K1(s) = (1/Z) sum_{k in (pi/12)Z} w(k) e^{i k s} (K1(0) = 1 exactly by +the normalization; K1 is exactly 24-periodic). + +DERIVATIVE PARITY DECOMPOSITION (derived, not copied). For real s, + + K1(s) = (1/Z) sum_k w(k) cos(k s) (imaginary parts cancel), + +so with E_n(s) = (1/Z) sum_k w(k) k^n cos(k s) and S_n(s) = (1/Z) sum_k w(k) +k^n sin(k s), + + K1^{(n)}(s) = (1/Z) sum_k w(k) k^n cos(k s + n pi/2) + n even: cos(ks + n pi/2) = (-1)^{n/2} cos(ks) -> K1^{(n)} = (-1)^{n/2} E_n(s) + n odd : cos(ks + n pi/2) = (-1)^{(n+1)/2} sin(ks) -> K1^{(n)} = (-1)^{(n+1)/2} S_n(s) + +(the odd case: n=1 gives cos(ks+pi/2) = -sin(ks); n=3 gives cos(ks+3pi/2) = ++sin(ks)). The odd-order derivatives are NONZERO SINE SUMS at s != 0. Zeroing +them is the RB_R2 rung-table defect of record (REVIEW_LPW/RB_R2_RUNGTABLE_ERRATA.md); +this library never special-cases them away. + +COVARIANCE BY MULTI-INDEX. For multi-indices alpha, beta in N^2 and +separation (dx, dy) = x - y, + + Cov(d^alpha f(x), d^beta f(y)) = (-1)^{|beta|} d^{alpha+beta} K(x-y) + = (-1)^{|beta|} K1^{(alpha1+beta1)}(dx) * K1^{(alpha2+beta2)}(dy). + +SECOND IMPLEMENTATION (no spectral sums). Poisson summation gives the +wrapped-kernel image sum + + K1(s) = (24/(sqrt(2 pi) Z)) sum_{n in Z} exp(-(s + 24 n)^2 / 2), + +and with probabilists' Hermite polynomials He_n (He_0=1, He_1=t, +He_{n+1} = t He_n - n He_{n-1}), (d/ds)^n e^{-s^2/2} = (-1)^n He_n(s) e^{-s^2/2}, + + K1^{(n)}(s) = (24/(sqrt(2 pi) Z)) sum_{n in Z} (-1)^n He_n(s+24n) e^{-(s+24n)^2/2}. + +CERTIFIED TAILS (computed in-program, interval arithmetic; see spectral_tail_1d +and image_tail_1d for the proved bounds). At the default cutoff KMAX = 30 the +spectral truncation tail is < 1e-60 for every derivative order n = 0..8; at the +default image window JMAX = 3 the image tail is < 1e-60 for n = 0..8 on the +whole post-wrap domain |s| <= 12 (and a fortiori on the adversarial set +|s| <= 1.62). + +TWO PATHS. mode='mpf' is the fast path (mpmath mpf at the configured dps). +mode='iv' is the certified path (mpmath iv.mpf interval arithmetic) and CARRIES +THE TRUNCATION TAILS: every interval returned by the iv path encloses the exact +periodized value, sum rounding + certified tail included. Theorem constants +must be consumed from the iv path (or from the mpf path plus the exposed tail +functions); the mpf path alone is diagnostics grade. + +FINITE-TORUS MOMENT DISCIPLINE. a_{2m} = (1/Z) sum_k w(k) k^{2m} are the exact +finite-torus moments. They are NOT the planar double factorials (2m-1)!!: e.g. +1 - a_2 = 576 * (sum_n n^2 e^{-288 n^2})/(sum_n e^{-288 n^2}) > 0 (~9.65e-123). +Planar values are exposed ONLY as labeled diagnostics (planar_moment). Consumers +must call moment_a (exact torus) and must not silently substitute planar values. + +FAIL-CLOSED CONTRACT. Every public entry validates its inputs and raises +CovInputError on violation; no bare asserts anywhere; nothing in this module +depends on wall clock, iteration order of unsorted containers, or environment. +Importing this module sets mp.dps and iv.dps to DEFAULT_DPS (100 >= 60); use +configure(dps=...) to change (dps < 60 is refused). + +Self-test: selftest_cov_exact.py (fail-closed; byte-identical under python and +python -O). Reference pattern consulted for cross-checking only: +SIDE24_gap_fill/scripts/verify_wp_witness_v1.py (kL/kE/kS); this implementation +is independently derived above. +""" +import mpmath +from mpmath import mp, mpf, iv + +SIDE = 24 # torus side +MIN_DPS = 60 # house floor +DEFAULT_DPS = 100 +DEFAULT_KMAX = 30 # spectral cutoff |k| <= KMAX +DEFAULT_JMAX = 3 # image window |j| <= JMAX +MAX_DERIV = 8 # certified derivative-order ceiling +TAIL_TARGET = mpf('1e-60') # certified-tail requirement at defaults + + +class CovInputError(ValueError): + """Raised on any contract violation (fail-closed input validation).""" + + +class CertificationError(RuntimeError): + """Raised when a certified bound cannot be established (fail-closed).""" + + +# -------------------------------------------------------------------------- +# configuration and lattice state +# -------------------------------------------------------------------------- +_CFG = { + 'dps': DEFAULT_DPS, + 'kmax': DEFAULT_KMAX, + 'jmax': DEFAULT_JMAX, + 'sine_drop': False, # MUTATION HOOK (self-test only; see selftest) +} +_LATT = None # (h, jmax_lattice, ks, ws, Z1) mpf +_LATT_IV = None # interval twin + + +def configure(dps=None, kmax=None, jmax=None): + """Reconfigure the engine. dps < MIN_DPS is refused (fail-closed). After + reconfiguration the certified-tail requirement (< 1e-60 for n = 0..MAX_DERIV) + is re-verified in-program; failure raises CertificationError.""" + global _LATT, _LATT_IV + if dps is not None: + if not (isinstance(dps, int) and dps >= MIN_DPS): + raise CovInputError("dps must be an int >= %d" % MIN_DPS) + _CFG['dps'] = dps + if kmax is not None: + if not (isinstance(kmax, int) and kmax >= 1): + raise CovInputError("kmax must be a positive int") + _CFG['kmax'] = kmax + if jmax is not None: + if not (isinstance(jmax, int) and jmax >= 0): + raise CovInputError("jmax must be a nonnegative int") + _CFG['jmax'] = jmax + mp.dps = _CFG['dps'] + iv.dps = _CFG['dps'] + _LATT = None + _LATT_IV = None + _build_lattice() + # fail-closed certification of the default truncation tails + for n in range(MAX_DERIV + 1): + t = spectral_tail_1d(n) + if not t < TAIL_TARGET: + raise CertificationError( + "spectral tail %.3e >= 1e-60 at n=%d under this configuration" % (float(t), n)) + ti = image_tail_1d(n, SIDE / 2) + if not ti < TAIL_TARGET: + raise CertificationError( + "image tail %.3e >= 1e-60 at n=%d, |s|<=12 under this configuration" % (float(ti), n)) + return dict(_CFG) + + +def _build_lattice(): + global _LATT, _LATT_IV + if _LATT is not None: + return _LATT + h = mp.pi / 12 + J = int(mp.ceil(mpf(_CFG['kmax']) / h)) + 1 # lattice index bound; |k| > KMAX omitted + ks = [j * h for j in range(-J, J + 1)] + ws = [mp.e ** (-k * k / 2) for k in ks] + Z1 = sum(ws) + _LATT = (h, J, ks, ws, Z1) + hiv = iv.pi / 12 + ks_iv = [j * hiv for j in range(-J, J + 1)] + ws_iv = [iv.exp(-k * k / 2) for k in ks_iv] + Z1_iv = sum(ws_iv, iv.mpf(0)) + _LATT_IV = (hiv, J, ks_iv, ws_iv, Z1_iv) + return _LATT + + +def get_config(): + _build_lattice() + return dict(_CFG) + + +def lattice(): + """Return (h, J, ks, ws, Z1): step pi/12, index bound J, lattice points, + masses, and the exact 1D normalization Z1 (mpf).""" + return _build_lattice() + + +# -------------------------------------------------------------------------- +# certified spectral truncation tail +# -------------------------------------------------------------------------- +def spectral_tail_1d(n): + """Rigorous upper bound (returned as mpf, computed with iv.mpf) on + + T_n = (1/Z1) sum_{|k| > KMAX, k in (pi/12)Z} e^{-k^2/2} |k|^n . + + Bound: let J0 = J+1 be the first omitted lattice index and h = pi/12. + For j >= J0 the term ratio obeys + + t_{j+1}/t_j = exp(-(2j+1) h^2/2) (1 + 1/j)^n + <= exp(-(2 J0+1) h^2/2) (1 + 1/J0)^n =: rho, + + because exp(-(2j+1)h^2/2) decreases in j and (1+1/j)^n decreases in j. + The routine VERIFIES rho < 1 in interval arithmetic (fail-closed), whence + the geometric series gives + + T_n <= 2 t_{J0} / (1 - rho) / Z1_lower, + t_{J0} = exp(-(J0 h)^2/2) (J0 h)^n . + + (Z1_lower >= 1 makes the division safe.)""" + if not (isinstance(n, int) and 0 <= n <= MAX_DERIV): + raise CovInputError("spectral_tail_1d: n must be an int in [0, %d]" % MAX_DERIV) + _build_lattice() + hiv, J, _, _, Z1_iv = _LATT_IV + J0 = J + 1 + k0 = J0 * hiv # first omitted |k| + rho = iv.exp(-(2 * J0 + 1) * hiv * hiv / 2) * (1 + 1 / iv.mpf(J0)) ** n + if not rho.b < 1: + raise CertificationError("spectral tail: ratio bound rho >= 1; cutoff not certified") + t0 = iv.exp(-k0 * k0 / 2) * k0 ** n + tail = 2 * t0 / (1 - rho) / Z1_iv.a # Z1_iv.a >= 1 is the safe divisor + return mpf(tail.b) + + +# -------------------------------------------------------------------------- +# certified image-sum truncation tail +# -------------------------------------------------------------------------- +def _he_absum(n): + """Sum of absolute values of the coefficients of the probabilists' Hermite + polynomial He_n, so |He_n(t)| <= _he_absum(n) * |t|^n for |t| >= 1. + Computed from the recurrence on absolute-coefficient polynomials: + A_0 = 1, A_1 = t, A_{k+1} = t A_k + k A_{k-1}; return A_n(1).""" + a, b = 1, 1 + if n == 0: + return 1 + for k in range(1, n): + a, b = b, b + k * a + return b + + +def image_tail_1d(n, smax): + """Rigorous upper bound (mpf, computed with iv.mpf) on the omitted-image + contribution to K1^{(n)}(s), uniform over |s| <= smax <= 12: + + T = (24/(sqrt(2 pi) Z1)) sum_{|j| > JMAX} |He_n(s+24j)| e^{-(s+24j)^2/2}. + + With A = 24 (JMAX+1) - smax every omitted image has |t| = |s+24j| >= A. + |He_n(t)| <= C_n |t|^n for |t| >= 1 with C_n = _he_absum(n). Consecutive + images on one side differ by 24, so for t >= A the term ratio obeys + + ((t+24)/t)^n exp(-((t+24)^2 - t^2)/2) <= (1 + 24/A)^n exp(-24 A) =: rho + + (since (t+24)^2 - t^2 = 48 t + 576 >= 48 A). The routine VERIFIES A >= 1 + and rho < 1 in interval arithmetic (fail-closed), whence + + T <= (24/(sqrt(2 pi) Z1_lower)) * 2 * C_n A^n e^{-A^2/2} / (1 - rho).""" + if not (isinstance(n, int) and 0 <= n <= MAX_DERIV): + raise CovInputError("image_tail_1d: n must be an int in [0, %d]" % MAX_DERIV) + smax = mpf(smax) + if not (0 <= smax <= SIDE / 2): + raise CovInputError("image_tail_1d: require 0 <= smax <= 12 (post-wrap domain)") + _build_lattice() + JIMG = _CFG['jmax'] + A = iv.mpf(24 * (JIMG + 1)) - iv.mpf(smax) + if not A.a >= 1: + raise CertificationError("image tail: A < 1; image window not certified") + Cn = iv.mpf(_he_absum(n)) + rho = (1 + 24 / A) ** n * iv.exp(-24 * A) + if not rho.b < 1: + raise CertificationError("image tail: ratio bound rho >= 1; window not certified") + first = Cn * A ** n * iv.exp(-A * A / 2) + Z1_iv = _LATT_IV[4] + norm = 24 / (iv.sqrt(2 * iv.pi) * Z1_iv.a) + tail = norm * 2 * first / (1 - rho) + return mpf(tail.b) + + +# -------------------------------------------------------------------------- +# 1D kernel derivatives: implementation 1 (spectral lattice sums) +# -------------------------------------------------------------------------- +def _check_order_s(n, s): + if not (isinstance(n, int) and 0 <= n <= MAX_DERIV): + raise CovInputError("derivative order n must be an int in [0, %d]" % MAX_DERIV) + return mpf(s) + + +def e_sum(n, s, mode='mpf'): + """E_n(s) = (1/Z) sum_k w(k) k^n cos(k s) (raw cosine sum, n >= 0).""" + s = _check_order_s(n, s) + _build_lattice() + if mode == 'mpf': + _, _, ks, ws, Z1 = _LATT + return sum(wj * (kj ** n) * mp.cos(kj * s) for kj, wj in zip(ks, ws)) / Z1 + if mode == 'iv': + _, _, ks, ws, Z1 = _LATT_IV + sv = iv.mpf(s) + tot = sum((wj * (kj ** n) * iv.cos(kj * sv) for kj, wj in zip(ks, ws)), iv.mpf(0)) + tau = spectral_tail_1d(n) + return tot / Z1 + iv.mpf([-tau, tau]) + raise CovInputError("mode must be 'mpf' or 'iv'") + + +def s_sum(n, s, mode='mpf'): + """S_n(s) = (1/Z) sum_k w(k) k^n sin(k s) (raw sine sum, n >= 0).""" + s = _check_order_s(n, s) + _build_lattice() + if mode == 'mpf': + _, _, ks, ws, Z1 = _LATT + return sum(wj * (kj ** n) * mp.sin(kj * s) for kj, wj in zip(ks, ws)) / Z1 + if mode == 'iv': + _, _, ks, ws, Z1 = _LATT_IV + sv = iv.mpf(s) + tot = sum((wj * (kj ** n) * iv.sin(kj * sv) for kj, wj in zip(ks, ws)), iv.mpf(0)) + tau = spectral_tail_1d(n) + return tot / Z1 + iv.mpf([-tau, tau]) + raise CovInputError("mode must be 'mpf' or 'iv'") + + +def k1_spectral(n, s, mode='mpf'): + """K1^{(n)}(s) by the spectral lattice sum, with the derived parity + decomposition: even n -> (-1)^{n/2} E_n(s); odd n -> (-1)^{(n+1)/2} S_n(s). + The odd (sine) branch is NEVER dropped (RB_R2 errata of record).""" + s = _check_order_s(n, s) + if n % 2 == 0: + if n % 4 == 0: + return e_sum(n, s, mode) + return -e_sum(n, s, mode) + # odd branch: sine sum with sign (-1)^{(n+1)/2} + if _CFG['sine_drop']: + # MUTATION HOOK (self-test only): the RB_R2 defect. Never active in + # library service; the self-test flips it to prove the guard rails. + return mpf(0) if mode == 'mpf' else iv.mpf(0) + if (n + 1) % 4 == 0: + return s_sum(n, s, mode) + return -s_sum(n, s, mode) + + +# -------------------------------------------------------------------------- +# 1D kernel derivatives: implementation 2 (wrapped-kernel image sums) +# -------------------------------------------------------------------------- +def hermite_he(n, t, mode='mpf'): + """Probabilists' Hermite polynomial He_n(t) by the recurrence + He_0 = 1, He_1 = t, He_{k+1} = t He_k - k He_{k-1}.""" + if not (isinstance(n, int) and n >= 0): + raise CovInputError("hermite_he: n must be a nonnegative int") + if mode == 'mpf': + t = mpf(t) + elif mode == 'iv': + t = iv.mpf(t) + else: + raise CovInputError("mode must be 'mpf' or 'iv'") + if n == 0: + return t ** 0 + a, b = t ** 0, t + if n == 1: + return b + for k in range(1, n): + a, b = b, t * b - k * a + return b + + +def _wrap24(s, mode): + """Reduce s modulo 24 into [-12, 12). K1 is exactly 24-periodic (every + lattice frequency k in (pi/12)Z satisfies e^{i k 24} = 1), so this is exact. + In iv mode the integer shift is verified by interval containment (fail-closed) + and the subtraction s - 24 n is exact in interval arithmetic.""" + if mode == 'mpf': + s = mpf(s) + n = int(mp.floor((s + 12) / 24)) + return s - 24 * n + sv = iv.mpf(s) + mid = mpf((sv.a + sv.b) / 2) + n = int(mp.floor((mid + 12) / 24)) + frac = (sv + 12) / 24 + if not (frac.a >= n and frac.b < n + 1): + raise CertificationError("wrap24: cannot certify the integer shift in iv mode") + return sv - 24 * n + + +def k1_image(n, s, mode='mpf'): + """K1^{(n)}(s) by the wrapped-kernel image sum (NO spectral sums): + (24/(sqrt(2 pi) Z)) sum_{|j| <= JMAX} (-1)^n He_n(t) e^{-t^2/2}, + t = s_wrapped + 24 j. The iv path adds the certified image tail.""" + s = _check_order_s(n, s) + _build_lattice() + JIMG = _CFG['jmax'] + if mode == 'mpf': + sw = _wrap24(s, 'mpf') + Z1 = _LATT[4] + tot = mpf(0) + for j in range(-JIMG, JIMG + 1): + t = sw + 24 * j + tot += hermite_he(n, t) * mp.e ** (-t * t / 2) + if n % 2 == 1: + tot = -tot # (d/ds)^n e^{-s^2/2} = (-1)^n He_n(s) e^{-s^2/2} + return tot * 24 / (mp.sqrt(2 * mp.pi) * Z1) + if mode == 'iv': + sw = _wrap24(s, 'iv') + Z1 = _LATT_IV[4] + tot = iv.mpf(0) + for j in range(-JIMG, JIMG + 1): + t = sw + 24 * j + tot = tot + hermite_he(n, t, 'iv') * iv.exp(-t * t / 2) + if n % 2 == 1: + tot = -tot + tau = image_tail_1d(n, SIDE / 2) + return tot * 24 / (iv.sqrt(2 * iv.pi) * Z1) + iv.mpf([-tau, tau]) + raise CovInputError("mode must be 'mpf' or 'iv'") + + +def k1(n, s, impl='spectral', mode='mpf'): + """K1^{(n)}(s); impl in {'spectral', 'image'}, mode in {'mpf', 'iv'}.""" + if impl == 'spectral': + return k1_spectral(n, s, mode) + if impl == 'image': + return k1_image(n, s, mode) + raise CovInputError("impl must be 'spectral' or 'image'") + + +# -------------------------------------------------------------------------- +# 2D covariance by multi-index +# -------------------------------------------------------------------------- +def _check_multi(idx, name): + if not (isinstance(idx, (tuple, list)) and len(idx) == 2 + and all(isinstance(v, int) and v >= 0 for v in idx)): + raise CovInputError("%s must be a pair of nonnegative ints" % name) + return (idx[0], idx[1]) + + +def cov(alpha, beta, dx, dy, impl='spectral', mode='mpf'): + """Cov(d^alpha f(x), d^beta f(y)) at separation (dx, dy) = x - y: + (-1)^{|beta|} K1^{(alpha1+beta1)}(dx) K1^{(alpha2+beta2)}(dy). + Cosine AND sine contributions are automatic via the parity decomposition. + Total order |alpha+beta| must be <= MAX_DERIV (certified-tail ceiling).""" + alpha = _check_multi(alpha, 'alpha') + beta = _check_multi(beta, 'beta') + n1, n2 = alpha[0] + beta[0], alpha[1] + beta[1] + if n1 + n2 > MAX_DERIV: + raise CovInputError("total order %d exceeds certified ceiling %d" + % (n1 + n2, MAX_DERIV)) + v = k1(n1, dx, impl, mode) * k1(n2, dy, impl, mode) + if (beta[0] + beta[1]) % 2 == 1: + v = -v + return v + + +# -------------------------------------------------------------------------- +# finite-torus moments (NEVER silently planar) +# -------------------------------------------------------------------------- +def moment_a(m, impl='spectral', mode='mpf'): + """Exact finite-torus moment a_{2m} = (1/Z) sum_k w(k) k^{2m} + = (-1)^m K1^{(2m)}(0). This is the ONLY moment entry point for theorem + use; planar double factorials are diagnostics (planar_moment).""" + if not (isinstance(m, int) and 0 <= 2 * m <= MAX_DERIV): + raise CovInputError("moment_a: need 0 <= 2m <= %d" % MAX_DERIV) + v = k1(2 * m, 0, impl, mode) + return v if m % 2 == 0 else -v + + +def moment_a_poisson(m, mode='mpf'): + """a_{2m} by the real-space Poisson identity + a_{2m} = (-1)^m [sum_n He_{2m}(24 n) e^{-288 n^2}] / [sum_n e^{-288 n^2}] + (|n| <= 50; the first omitted image has |t| = 1224, so omitted terms are + below e^{-749088} up to polynomial factors, certified by + moment_poisson_tail). Independent of the spectral lattice path.""" + if not (isinstance(m, int) and 0 <= 2 * m <= MAX_DERIV): + raise CovInputError("moment_a_poisson: need 0 <= 2m <= %d" % MAX_DERIV) + if mode == 'mpf': + num = mpf(0); den = mpf(0) + for j in range(-50, 51): + t = mpf(24 * j) + e = mp.e ** (-t * t / 2) + num += hermite_he(2 * m, t) * e + den += e + v = num / den + elif mode == 'iv': + num = iv.mpf(0); den = iv.mpf(0) + for j in range(-50, 51): + t = iv.mpf(24 * j) + e = iv.exp(-t * t / 2) + num = num + hermite_he(2 * m, t, 'iv') * e + den = den + e + tau = moment_poisson_tail(m) + # omitted |n|>50 terms enter BOTH num and den; |num|/den = a_{2m} and + # C_{2m} >= 1, so 2*tau rigorously covers the ratio perturbation. + v = num / den + 2 * iv.mpf([-tau, tau]) + else: + raise CovInputError("mode must be 'mpf' or 'iv'") + return v if m % 2 == 0 else -v + + +def moment_poisson_tail(m): + """Certified bound on the omitted |n| > 50 terms of the Poisson moment + sums: |He_{2m}(24n)| e^{-(24n)^2/2} <= C_{2m} |24n|^{2m} e^{-(24n)^2/2}, + first omitted |t| = 1224, ratio <= (1+24/1224)^{2m} e^{-24*1224} < 1. + The numerator tail is normalized by sum_n e^{-288 n^2} >= 1 (safe).""" + if not (isinstance(m, int) and 0 <= 2 * m <= MAX_DERIV): + raise CovInputError("moment_poisson_tail: need 0 <= 2m <= %d" % MAX_DERIV) + A = iv.mpf(1224) + Cn = iv.mpf(_he_absum(2 * m)) + rho = (1 + 24 / A) ** (2 * m) * iv.exp(-24 * A) + if not rho.b < 1: + raise CertificationError("moment Poisson tail: rho >= 1") + return mpf((2 * Cn * A ** (2 * m) * iv.exp(-A * A / 2) / (1 - rho)).b) + + +def one_minus_a2_poisson(mode='mpf'): + """The exact identity 1 - a_2 = 576 (sum_n n^2 e^{-288 n^2})/(sum_n e^{-288 n^2}) + (real-space form; > 0 strictly: the torus is never planar).""" + if mode == 'mpf': + num = mpf(0); den = mpf(0) + for j in range(-50, 51): + e = mp.e ** (mpf(-288) * j * j) + num += j * j * e + den += e + return 576 * num / den + if mode == 'iv': + num = iv.mpf(0); den = iv.mpf(0) + for j in range(-50, 51): + e = iv.exp(iv.mpf(-288) * j * j) + num = num + j * j * e + den = den + e + # omitted |n|>50 terms enter BOTH num and den; C_4 >= C_2 = 1 makes the + # m=2 tail dominate each, and 3x covers the ratio perturbation. + tau = moment_poisson_tail(2) + return 576 * num / den + 3 * iv.mpf([-tau, tau]) + raise CovInputError("mode must be 'mpf' or 'iv'") + + +def planar_moment(m): + """DIAGNOSTIC ONLY (labeled): the planar moment (2m-1)!! = E[N(0,1)^{2m}]. + Forbidden as a theorem input on the finite torus; the exact finite-torus + moment is moment_a(m).""" + if not (isinstance(m, int) and m >= 0): + raise CovInputError("planar_moment: m must be a nonnegative int") + v = 1 + for j in range(1, m + 1): + v *= (2 * j - 1) + return mpf(v) + + +# apply default configuration at import (sets mp.dps/iv.dps, builds lattice, +# and certifies the default tails fail-closed) +configure(dps=DEFAULT_DPS, kmax=DEFAULT_KMAX, jmax=DEFAULT_JMAX) diff --git a/experiments/h3_band_v2/h3_band_v2.py b/experiments/h3_band_v2/h3_band_v2.py new file mode 100644 index 000000000..aa7949396 --- /dev/null +++ b/experiments/h3_band_v2/h3_band_v2.py @@ -0,0 +1,1518 @@ +#!/usr/bin/env python3 +"""H3 band certificate, numbered successor v2: a ceiling and a floor for E[G_r], r in (0, 0.05]. + +Successor of the frozen programs on main, under +incoming/claude-bundle-v3-native-copies-20260928/K3_SIDE24_LB/UPPER2D/H3_closure/: +h3_band_ceil.py (sha256 b97c5428...) and h3_band_floor.py (sha256 a907eeed...), author Kimi. +Those two files share their law pipeline (CB0-CB6) byte for byte; this program keeps that +pipeline once, stage by stage, and runs either envelope. README.md lists every change and why. + +Arithmetic. mpmath iv (outward-rounded interval arithmetic) at dps 100, as in the frozen +programs. Every load-bearing endpoint comes from an interval. Round-to-nearest mp is used only +where the result is re-checked in iv (the midpoint preconditioner, the search for the majorant +growth R) or is not load-bearing (the CB7 direct point law, which is a cross-check). The exact +pole-structure lemma (CB2s) uses fractions.Fraction only. Certified decimals are rounded outward +exactly. The Monte Carlo diagnostic (--mc) is NON-CERTIFYING and writes to stderr only. + +Usage: + h3_band_v2.py structure the exact pole-structure lemma alone + h3_band_v2.py ceil [--transcript F] [--record R | --verify R] [--mc] + h3_band_v2.py floor [--transcript F] [--record R | --verify R] [--mc] + h3_band_v2.py verify-record --transcript F --results R +Negative controls: --mutation NAME, one of MUTATIONS below. A checker prints one summary line +and exits nonzero on failure. +""" +import argparse +import hashlib +import json +import platform +import sys +import time +import types +from collections import defaultdict +from fractions import Fraction +from functools import lru_cache +from math import comb, factorial +from pathlib import Path + +try: + import mpmath + from mpmath import mp, iv, mpf +except ImportError: # the certificate needs mpmath; the structure lemma does not + mpmath = None + +HERE = Path(__file__).resolve().parent +COV_EXACT_SHA256 = 'f08c1c5f653f2ffd2e85d39e1112f80d15db04d15f932e5ee42db2ca3553e783' +PROGRAM = 'experiments/h3_band_v2/h3_band_v2.py' + +MUTATIONS = { + 'structure_pin_pole': 'double the -1/2 weight of S(1,0) in pin V5 to -1; CB2s (a) must find a z^-2 term', + 'structure_soft_pole2': 'soft Y2 = M(2,0) z^-2; CB2s (b) must find a pole of order 2', + 'structure_soft_span': 'soft Y0 = M(0,2) z^-1, whose pole part f_yy is not in the span of V^(0); CB2s (c) must fail', + 'structure_window': 'conditioning vector with zV[3] = 1; CB2s (d) must find a mu pole', + 'planar_moments': 'planar (2m-1)!! in the consumed moment table A2; the CB1d spectral cross-check must fail', + 'wrong_window': 'zV[5] = -1/12; the CB3 limit identities must fail', + 'signed_maxima': 'B1 regression (signed upper endpoints as majorants); the CB4 dominance guard must fail', + 'wide_first_cell': 'R-majorant first cell (0, 0.0025] as in the frozen programs; the CB4 convergence guard must fail', + 'one_tail': 'B3 regression (one Gaussian tail for a two-sided event) in the event bound every call site uses; the CB5 event self-test must fail', + 'tail_helper_one_side': 'the tail helper ignores its number of sides; the CB5 Gordon self-test must fail', + 'no_preconditioner': 'identity in place of the midpoint inverse; the CB6 contraction guard must fail', + 'cap_raised': 'p and pbar doubled; the ceiling consumption cap must fail', + 'drop_fail_piece': 'ceiling without its type-fail pieces; not detectable inside the run, --verify must fail on the cells', + 'old_remainder': 'B2 regression (the frozen dGinv law correction); --verify must fail on the cells', +} + + +def mutated_tables(): + """The exact functional tables and conditioning vector, with any structure mutation applied. + Both the lemma CB2s and the interval pipeline are built from what this returns.""" + pins = [list(f) for f in PINS_EXACT] + soft = [list(f) for f in SOFT_EXACT] + zv = list(ZV_EXACT) + if MUT == 'structure_pin_pole': + pins[5] = [pins[5][0], pins[5][1], ('S', (1, 0), F(-1), -2), pins[5][3]] + if MUT == 'structure_soft_pole2': + soft[2] = [('M', (2, 0), F(1), -2)] + if MUT == 'structure_soft_span': + soft[0] = [('M', (0, 2), F(1), -1)] + if MUT == 'structure_window': + zv[3] = F(1) + if MUT == 'wrong_window': + zv[5] = F(-1, 12) + return pins, soft, zv +MUT = '' +OUT = [] + + +def emit(s): + OUT.append(s) + print(s) + + +def note(s): + # stderr only: timings and diagnostics that are not part of the transcript + print(s, file=sys.stderr) + + +class CheckFailed(Exception): + pass + + +def ck(cond, msg): + if not cond: + emit('CK_FAIL: ' + msg) + raise CheckFailed(msg) + + +# ------------------------------------------------------------------ exact helpers +def mpf_to_fraction(x): + """Exact value of a finite mpf as a Fraction (mpf is a binary float: (-1)^sign man 2^exp).""" + x = mpf(x) + if not mpmath.isfinite(x): + raise ValueError('nonfinite mpf') + sign, man, exp, _ = x._mpf_ + man = -int(man) if sign else int(man) + return Fraction(man * 2 ** exp) if exp >= 0 else Fraction(man, 2 ** (-exp)) + + +def decimal_outward(x, direction, sig=15): + """Round x (mpf or Fraction) to sig significant digits, up ('up') or down ('down'), exactly.""" + q = x if isinstance(x, Fraction) else mpf_to_fraction(x) + if q == 0: + return '0' + neg = q < 0 + a = -q if neg else q + e10 = len(str(a.numerator // a.denominator)) - 1 if a >= 1 else -len(str(a.denominator // a.numerator)) + while Fraction(10) ** e10 > a: + e10 -= 1 + while Fraction(10) ** (e10 + 1) <= a: + e10 += 1 + scale = Fraction(10) ** (sig - 1 - e10) + v = a * scale + up = (direction == 'up') != neg + n = -((-v.numerator) // v.denominator) if up else v.numerator // v.denominator + digits = str(n) + point = len(digits) - (sig - 1 - e10) + if point <= 0: + s = '0.' + '0' * (-point) + digits + elif point >= len(digits): + s = digits + '0' * (point - len(digits)) + else: + s = digits[:point] + '.' + digits[point:] + return ('-' if neg else '') + s + + +def frac_label(q): + return str(q.numerator) if q.denominator == 1 else '%d/%d' % (q.numerator, q.denominator) + + +# ------------------------------------------------- CB2s: exact Laurent pole structure +# The functionals of record, as (station, multi-index, exact weight, power of z). M sits at the +# origin and S at z e_x (DISP below: S - M = +z e_x). The interval series of CB2 are generated +# from these tables, so the lemma is about exactly the functionals the certificate uses. +F = Fraction +PINS_EXACT = [ + [('M', (0, 0), F(1), 0)], + [('M', (1, 0), F(1), 0)], + [('M', (0, 1), F(1), 0)], + [('S', (1, 0), F(1), -1), ('M', (1, 0), F(-1), -1)], + [('S', (0, 1), F(1), -1), ('M', (0, 1), F(-1), -1)], + [('S', (0, 0), F(1), -3), ('M', (0, 0), F(-1), -3), ('S', (1, 0), F(-1, 2), -2), ('M', (1, 0), F(-1, 2), -2)], +] +SOFT_EXACT = [ + [('M', (0, 2), F(1), 0)], + [('S', (0, 2), F(1), 0)], + [('M', (2, 0), F(1), -1)], + [('S', (2, 0), F(1), -1)], + [('M', (1, 1), F(1), -1)], + [('S', (1, 1), F(1), -1)], +] +ZV_EXACT = [F(6, 5), F(0), F(0), F(0), F(0), F(-1, 6)] + + +def laurent_expand(func, top): + """Exact Laurent expansion in z of a functional, through z^top: {power: {beta: weight}}. + S(alpha) = d^alpha f(z e_x) = sum_q z^q / q! d^{alpha + (q, 0)} f(0).""" + out = defaultdict(lambda: defaultdict(Fraction)) + for station, (ax, ay), w, k in func: + if station == 'M': + out[k][(ax, ay)] += w + else: + for q in range(0, top - k + 1): + out[k + q][(ax + q, ay)] += w / factorial(q) + return {p: {b: c for b, c in row.items() if c != 0} for p, row in out.items()} + + +def solve_in_span(rows, target): + """Exact coefficients x with sum_k x_k rows[k] = target (dicts beta -> Fraction), or None.""" + keys = sorted(set().union(*[set(r) for r in rows], set(target))) + n = len(rows) + A = [[rows[k].get(b, F(0)) for k in range(n)] + [target.get(b, F(0))] for b in keys] + piv_cols, r = [], 0 + for c in range(n): + p = next((i for i in range(r, len(A)) if A[i][c] != 0), None) + if p is None: + continue + A[r], A[p] = A[p], A[r] + A[r] = [v / A[r][c] for v in A[r]] + for i in range(len(A)): + if i != r and A[i][c] != 0: + f = A[i][c] + A[i] = [a - f * b for a, b in zip(A[i], A[r])] + piv_cols.append(c) + r += 1 + if any(all(v == 0 for v in A[i][:n]) and A[i][n] != 0 for i in range(len(A))): + return None + x = [F(0)] * n + for i, c in enumerate(piv_cols): + x[c] = A[i][n] + return x + + +def structure_lemma(pins, soft, zv): + """CB2s. Proves, in exact rational arithmetic, the pole structure the certificate relies on: + (a) every pin functional V_k has no negative powers of z (so the pin Gram G(z) is analytic); + (b) every soft functional Y_a has no power below z^-1; + (c) the z^-1 part of Y equals E V^(0) for an exact matrix E, where V^(0) is the z^0 part of V; + (d) E zV = 0. + Then Ytilde = Y - E V / z is regular, E[Y | V = zV] = E[Ytilde | V = zV] + E zV / z has no + pole, and Cov(Y | V) = Cov(Ytilde | V) is analytic wherever G(z) is invertible. So the + z^-1 coefficient of mu and the z^-2, z^-1 coefficients of Sigma vanish exactly.""" + pe = [laurent_expand(f, 1) for f in pins] + se = [laurent_expand(f, 1) for f in soft] + for k, e in enumerate(pe): + for p, row in e.items(): + ck(not (p < 0 and row), 'CB2s (a) pin V%d must have no z^%d term' % (k, p)) + for a, e in enumerate(se): + for p, row in e.items(): + ck(not (p < -1 and row), 'CB2s (b) soft Y%d must have no z^%d term' % (a, p)) + v0 = [e.get(0, {}) for e in pe] + E = [] + for a, e in enumerate(se): + x = solve_in_span(v0, e.get(-1, {})) + ck(x is not None, 'CB2s (c) soft Y%d pole part must lie in the span of V^(0)' % a) + E.append(x) + for a in range(len(soft)): + ck(sum(E[a][k] * zv[k] for k in range(len(zv))) == 0, 'CB2s (d) E zV = 0 required, row %d' % a) + return E + + +# ------------------------------------------------------------------ interval helpers +def ub(x): + """Exact upper bound on |x| for an iv interval or an mpf.""" + if hasattr(x, 'b'): + return max(abs(mpf(x.a)), abs(mpf(x.b))) + return abs(mpf(x)) + + +def up(x): + """Upper endpoint of an iv value, as an exact mpf.""" + return mpf(x.b) + + +def iv_frac(q): + """Enclosure of an exact rational.""" + return iv.mpf(q.numerator) / q.denominator + + +def run_certificate(mode, h2_dir, mc=False): + """Build the law pipeline (CB0-CB6) and evaluate the ceiling or the floor on every cell.""" + t_start = time.time() + emit('H3 BAND V2 %s (successor of the frozen h3_band_%s.py)' % (mode.upper(), 'ceil' if mode == 'ceil' else 'floor')) + emit('mode dps=100 MUTATION=' + (MUT if MUT else 'none')) + + # ----------------------------------------------------------- CB2s: exact structure + pins_exact, soft_exact, zv_exact = mutated_tables() + if MUT == 'wrong_window': + emit('CB3 MUTATION wrong_window applied (must be caught)') + E_sel = structure_lemma(pins_exact, soft_exact, zv_exact) + emit('CB2s PASS: exact Laurent structure: pins regular, soft poles z^-1 only, Y^(-1) = E V^(0), E zV = 0') + + mp.dps = 100 + iv.dps = 100 + # ------------------------------------------------------------------ CB0: hashes + h2 = Path(h2_dir) + src = (h2 / 'cov_exact.py').read_bytes() + digest = hashlib.sha256(src).hexdigest() + manifest = {} + for ln in (h2 / 'MANIFEST.sha256').read_text().splitlines(): + hh, fn = ln.split() + manifest[fn] = hh + ck(digest == COV_EXACT_SHA256 and manifest.get('cov_exact.py') == COV_EXACT_SHA256, 'CB0 cov_exact hash') + # execute exactly the verified bytes: no import machinery, so no bytecode cache can stand in + ce = types.ModuleType('cov_exact') + ce.__file__ = str(h2 / 'cov_exact.py') + exec(compile(src, ce.__file__, 'exec'), ce.__dict__) # sets mp.dps = iv.dps = 100, certifies its tails + ck(mp.dps == 100 and iv.dps == 100, 'CB0 precision dps=100') + emit('CB0 PASS: H2 library cov_exact.py sha256 %s... matches the pinned digest' % COV_EXACT_SHA256[:16]) + + REQ = mpf('1.615489267643502474048123284509081575382') # c_Z of record + RMAX = F(1, 20) + + # ------------------------------------------------------ CB1: exact-torus moments + # a_{2m} = (-1)^m sum_n He_{2m}(24 n) e^{-288 n^2} / sum_n e^{-288 n^2} (Poisson) + def He_iv(n, x): + if n == 0: + return iv.mpf(1) + a, b = iv.mpf(1), iv.mpf(x) + for k in range(2, n + 1): + a, b = b, x * b - (k - 1) * a + return b + + def moment_a2m(m): + # |n| <= 3 exactly, |n| >= 4 by a certified tail: + # |He_k(x)| <= (|x| + sqrt k)^k, and for n >= 4 consecutive terms (24n + c)^{2m} e^{-288 n^2} + # shrink by at least (5/4)^{2m} e^{-2592} < 1/2, so each tail is at most twice its first + # term, on both sides of 0. + num = iv.mpf(0) + den = iv.mpf(0) + for n in range(-3, 4): + w = iv.exp(-iv.mpf(288) * n * n) + num += He_iv(2 * m, iv.mpf(24 * n)) * w + den += w + ck(((iv.mpf(5) / 4) ** (2 * m) * iv.exp(-iv.mpf(2592))).b < mpf(1) / 2, 'CB1 moment tail ratio m=%d' % m) + first = iv.exp(-iv.mpf(288 * 16)) + tn = 4 * (96 + iv.sqrt(iv.mpf(2 * m))) ** (2 * m) * first + td = 4 * first + num = num + iv.mpf((-up(tn), up(tn))) + den = den + iv.mpf((0, up(td))) + return ((-1) ** m) * num / den + + MMAX = 40 + A2 = [moment_a2m(m) for m in range(MMAX + 1)] + DEV2 = ce.one_minus_a2_poisson(mode='iv') + num4 = iv.mpf(0) + den4 = iv.mpf(0) + for n in range(-3, 4): + w = iv.exp(-iv.mpf(288) * n * n) + t = iv.mpf(24 * n) + num4 += (He_iv(4, t) - 3) * w # He_4(t) - 3 = t^2 (t^2 - 6) >= 0 at t in 24Z + den4 += w + first = iv.exp(-iv.mpf(288 * 16)) + t4 = 4 * ((96 + iv.mpf(2)) ** 4 + 3) * first + DEV4 = (num4 + iv.mpf((0, up(t4)))) / (den4 + iv.mpf((0, up(4 * first)))) + if MUT == 'planar_moments': + A2 = [iv.mpf(factorial(2 * m)) / iv.mpf(factorial(m) * 2 ** m) for m in range(MMAX + 1)] + emit('CB1 MUTATION planar_moments applied to the consumed table A2 (must be caught)') + ck(DEV2.a > 0 and DEV2.b < mpf('1e-100'), 'CB1a exact-torus 0 < 1-a2 < 1e-100') + ck(DEV4.a > 0 and DEV4.b < mpf('1e-100'), 'CB1b exact-torus 0 < a4-3 < 1e-100') + ck(all(ub(A2[m] - ce.moment_a_poisson(m, mode='iv')) < mpf('1e-90') for m in range(5)), + 'CB1c Poisson moments agree with cov_exact for 2m <= 8') + # CB1d: the consumed table against an independent spectral enclosure, every m <= MMAX: + # a_{2m} = (1/Z) sum_{k in (pi/12)Z} e^{-k^2/2} k^{2m}; the lattice holds |j| <= Jl, and beyond it + # consecutive terms shrink by at most q = e^{-(2 Jl + 3) h^2 / 2} (1 + 1/(Jl+1))^{2m} < 1. + hstep_s, Jl_s, ks_s, ws_s, Z1_s = ce._LATT_IV + + def moment_spectral(m): + n = 2 * m + tot = sum((ww * kk ** n for kk, ww in zip(ks_s, ws_s)), iv.mpf(0)) + j0 = Jl_s + 1 + k0 = j0 * hstep_s + q = iv.exp(-(2 * j0 + 1) * hstep_s * hstep_s / 2) * (1 + iv.mpf(1) / j0) ** n + ck(q.b < 1, 'CB1d spectral tail ratio < 1, m=%d' % m) + tail = 2 * iv.exp(-k0 * k0 / 2) * k0 ** n / (1 - q) + return (tot + iv.mpf((0, up(tail)))) / Z1_s + + for m in range(MMAX + 1): + sp = moment_spectral(m) + ck(A2[m].a <= sp.b and sp.a <= A2[m].b, 'CB1d consumed a_{2m} meets the spectral enclosure, m=%d' % m) + emit('CB1 PASS: moments a_{2m}, m<=%d, exact finite torus with certified image tails (1-a2 = %s)' % (MMAX, mp.nstr(up(DEV2), 3))) + + # ------------------------------------------------------------ CB5 (early): Phi and tails + SQRT2PI = iv.sqrt(2 * iv.pi) + + def iv_phi(x): + return iv.exp(-(x * x) / 2) / SQRT2PI + + def iv_Phi(x): + # standard normal CDF on an interval x. |x| <= 6: power series with a certified geometric + # tail. |x| > 6: Gordon-Mills bounds (both monotone). Intervals that straddle +-6 are split + # by monotonicity. Endpoints stay mpf (O2: the frozen float() conversions are removed). + if (x.a < 6 < x.b) or (x.a < -6 < x.b): + lo = iv_Phi(iv.mpf(x.a)) + hi = iv_Phi(iv.mpf(x.b)) + return iv.mpf((mpf(lo.a), mpf(hi.b))) + if x.a >= 6: + hi = iv_phi(iv.mpf(x.a)) / iv.mpf(x.a) + lo = iv_phi(iv.mpf(x.b)) * (1 / iv.mpf(x.b) - 1 / iv.mpf(x.b) ** 3) + return iv.mpf(1) - iv.mpf((mpf(lo.a), mpf(hi.b))) + if x.b <= -6: + hi = iv_phi(iv.mpf(x.b)) / iv.mpf(-x.b) + lo = iv_phi(iv.mpf(x.a)) * (1 / iv.mpf(-x.a) - 1 / iv.mpf(-x.a) ** 3) + return iv.mpf((mpf(lo.a), mpf(hi.b))) + ck(x.a >= -6 and x.b <= 6, 'CB5 Phi series domain') + # term ratio -x^2 (2k+1) / (2 (k+1) (2k+3)): modulus < 0.9 for k >= 21 when |x| <= 6 + tot = iv.mpf(0) + K = 0 + term = x + while K < 120: + tot += term + K += 1 + term = -term * x * x * (2 * K - 1) / (2 * K * (2 * K + 1)) + if K >= 22 and abs(term.b) < 1e-90 and abs(term.a) < 1e-90: + break + tail = up(iv.mpf(ub(term)) / (1 - iv.mpf('0.9'))) + return iv.mpf(1) / 2 + (tot + iv.mpf((-tail, tail))) / SQRT2PI + + def gauss_tail_upper(t0_lo, sides): + """Upper bound on P(Z > t0) (sides=1) or P(|Z| > t0) (sides=2), Z ~ N(0,1), t0 >= t0_lo.""" + p = up(iv_Phi(iv.mpf(-t0_lo))) + if MUT == 'tail_helper_one_side': + return p + return p * sides + + # self-test against an enclosure that does not use iv_Phi: Gordon's inequality + # 1 - Phi(t) >= phi(t) t / (t^2 + 1) for t > 0 + for t in ('2.5', '4', '6.5', '9'): + tt = iv.mpf(t) + gordon = mpf((iv_phi(tt) * tt / (tt * tt + 1)).a) + ck(gauss_tail_upper(mpf(t), 1) >= gordon, 'CB5 one-sided tail helper dominates the Gordon lower bound at %s' % t) + ck(gauss_tail_upper(mpf(t), 2) >= 2 * gordon, 'CB5 two-sided tail helper dominates twice the Gordon lower bound at %s' % t) + emit('CB5 PASS: certified Phi; Gaussian tail helpers dominate the Gordon lower bounds') + + def t0_lower(num_iv, var_iv): + ck(up(var_iv) > 0, 'CB5 variance ceiling positive before its square root') + return mpf((num_iv / iv.sqrt(iv.mpf(up(var_iv)))).a) + + def event_prob_upper(mu_x, var_x, kind, edge, margin, what): + """Upper bound on P(X > edge) ('gt'), P(X < edge) ('lt') or P(|X| > edge) ('abs'), for + X ~ N(m, v) and every m in mu_x, v in var_x, the edge beyond the mean by the margin in units + of the standard deviation. Every Gaussian tail of the certificate goes through here, and the + number of tails follows from the event kind (B3: two for 'abs').""" + e = iv.mpf(edge) + if kind == 'gt': + t0 = t0_lower(e - iv.mpf(up(mu_x)), var_x) + sides = 1 + elif kind == 'lt': + t0 = t0_lower(iv.mpf(mpf(mu_x.a)) - e, var_x) + sides = 1 + else: + t0 = t0_lower(e - iv.mpf(ub(mu_x)), var_x) + sides = 2 if MUT != 'one_tail' else 1 + ck(t0 > margin, what + ' margin') + return gauss_tail_upper(t0, sides) + + # self-test of the event bounds on exact point laws, against Gordon lower bounds that do not use + # iv_Phi: P(X > t) and P(X < -t) >= g(t) for N(0, 1); P(|X| > t) >= 2 g(t); and for N(1/2, 1), + # P(|X| > t) >= g(t - 1/2) + g(t + 1/2), where g(t) = phi(t) t / (t^2 + 1) + def gordon(t): + return mpf((iv_phi(t) * t / (t * t + 1)).a) + + half = iv.mpf(1) / 2 + for t in ('2.5', '4', '6.5', '9'): + tt = iv.mpf(t) + z0 = iv.mpf(0) + one = iv.mpf(1) + g = gordon(tt) + ck(event_prob_upper(z0, one, 'gt', tt, 2, 'CB5 self-test') >= g, 'CB5 event bound P(X > t) self-test at %s' % t) + ck(event_prob_upper(z0, one, 'lt', -tt, 2, 'CB5 self-test') >= g, 'CB5 event bound P(X < -t) self-test at %s' % t) + ck(event_prob_upper(z0, one, 'abs', tt, 2, 'CB5 self-test') >= 2 * g, 'CB5 event bound P(|X| > t) self-test at %s' % t) + ck(event_prob_upper(half, one, 'abs', tt, 1, 'CB5 self-test') >= gordon(tt - half) + gordon(tt + half), + 'CB5 event bound P(|X| > t), mean 1/2, self-test at %s' % t) + emit('CB5 PASS: Gaussian event bounds dominate the Gordon lower bounds on exact point laws') + + # ------------------------------------------------------------- CB2: series machinery + J = 60 # Taylor order kept for mu, Sigma + NK = J + 2 # Gram blocks kept to order J+2; Cauchy tails start at J+3 + LMIN = -8 + RHO = mpf('0.11') # majorant circle radius (a choice; every bound below is certified on it) + + def ser_zero(): + return defaultdict(lambda: iv.mpf(0)) + + def ser_monom(c, k): + s = ser_zero() + s[k] = c + return s + + def ser_mul(A, B): + C = ser_zero() + for ja, ca in A.items(): + if ca == 0: + continue + for jb, cb in B.items(): + if cb == 0: + continue + C[ja + jb] += ca * cb + return C + + def ser_add(A, B): + C = ser_zero() + for j, c in A.items(): + C[j] += c + for j, c in B.items(): + C[j] += c + return C + + def ser_scale(A, c): + C = ser_zero() + for j, v in A.items(): + C[j] = v * c + return C + + def ser_trim(A): + return defaultdict(lambda: iv.mpf(0), + {j: c for j, c in A.items() if LMIN <= j <= NK and not (c.a == 0 and c.b == 0)}) + + # kernel series K1^{(n)}(+-z): coefficient j is (-1)^{(n+j)/2} a_{n+j} / j! (n+j even). + # A Gram entry at order k uses kernel orders up to k + 6 - n (the largest pair weight, pole + # shift plus x-order, is 6), so order NK needs j <= NK + 6 - n: range(0, J + 9 - n). + def kernel_ser(n, parity): + s = ser_zero() + for j in range(0, J + 9 - n): + m = n + j + if m % 2 == 1: + continue + ck(m // 2 <= MMAX, 'CB2 kernel order within the moment table') + val = ((-1) ** (m // 2)) * A2[m // 2] / factorial(j) + if parity < 0 and j % 2 == 1: + val = -val + s[j] = val + return s + + KSER = {} + + def get_kser(n, parity): + if (n, parity) not in KSER: + KSER[(n, parity)] = kernel_ser(n, parity) + return KSER[(n, parity)] + + def cov_ser(alpha, beta, d): + # cov(alpha, beta, d z, 0) = (-1)^{|beta|} K1^{(n1)}(d z) K1^{(n2)}(0) + n1 = alpha[0] + beta[0] + n2 = alpha[1] + beta[1] + sgn = (-1) ** (beta[0] + beta[1]) + if d == 0: + v1 = ser_zero() + if n1 % 2 == 0: + v1[0] = ((-1) ** (n1 // 2)) * A2[n1 // 2] + else: + v1 = get_kser(n1, d) + k0 = iv.mpf(0) + if n2 % 2 == 0: + k0 = ((-1) ** (n2 // 2)) * A2[n2 // 2] + return ser_scale(v1, sgn * k0) + + def to_iv_func(func): + return [(st, al, ser_monom(iv_frac(w), k)) for st, al, w, k in func] + + PinsV = [to_iv_func(f) for f in pins_exact] + SoftY = [to_iv_func(f) for f in soft_exact] + DISP = {('M', 'M'): 0, ('S', 'S'): 0, ('S', 'M'): 1, ('M', 'S'): -1} + + def gram_entry_ser(FA, FB): + tot = ser_zero() + for (xa, aa, ca) in FA: + for (xb, ab, cb) in FB: + tot = ser_add(tot, ser_mul(ser_mul(ca, cb), cov_ser(aa, ab, DISP[(xa, xb)]))) + return ser_trim(tot) + + t1 = time.time() + GP_ser = [[gram_entry_ser(PinsV[i], PinsV[j]) for j in range(6)] for i in range(6)] + GY_ser = [[gram_entry_ser(SoftY[i], SoftY[j]) for j in range(6)] for i in range(6)] + GYV_ser = [[gram_entry_ser(SoftY[i], PinsV[j]) for j in range(6)] for i in range(6)] + note('CB2 series built in %.1fs' % (time.time() - t1)) + bad = 0 + for i in range(6): + for j in range(6): + for k, c in GP_ser[i][j].items(): + if k < 0 and ub(c) > 1e-60: + bad += 1 + for k, c in GYV_ser[i][j].items(): + if k <= -2 and ub(c) > 1e-60: + bad += 1 + for k, c in GY_ser[i][j].items(): + if k <= -3 and ub(c) > 1e-60: + bad += 1 + ck(bad == 0, 'CB2 numeric Laurent pole pattern agrees with CB2s') + emit('CB2 PASS: Gram series to order %d; numeric pole pattern agrees with the exact lemma CB2s' % NK) + + # ---------------------------------------------------- CB3: series division -> mu, Sigma + def mat_coeff(Gser, j): + return [[Gser[a][b].get(j, iv.mpf(0)) for b in range(6)] for a in range(6)] + + def iv_solve6(GG, rhs): + # interval Gauss-Jordan; every step is an interval extension, so the result encloses + # the exact solution (a division by an interval containing 0 would surface downstream) + A = [GG[i][:] + [rhs[i]] for i in range(6)] + for c in range(6): + p = max(range(c, 6), key=lambda rr: abs(A[rr][c].a + A[rr][c].b)) + A[c], A[p] = A[p], A[c] + for rr in range(6): + if rr != c: + f = A[rr][c] / A[c][c] + A[rr] = [A[rr][j] - f * A[c][j] for j in range(7)] + return [A[i][6] / A[i][i] for i in range(6)] + + def matvec(G, v): + return [sum((G[i][k] * v[k] for k in range(6)), iv.mpf(0)) for i in range(6)] + + def matmul(A, B): + return [[sum((A[i][k] * B[k][j] for k in range(6)), iv.mpf(0)) for j in range(6)] for i in range(6)] + + def matadd(A, B): + return [[A[i][j] + B[i][j] for j in range(6)] for i in range(6)] + + def matneg(A): + return [[-A[i][j] for j in range(6)] for i in range(6)] + + def transpose(A): + return [[A[j][i] for j in range(6)] for i in range(6)] + + def matzero(): + return [[iv.mpf(0)] * 6 for _ in range(6)] + + def eye(): + return [[iv.mpf(1) if i == j else iv.mpf(0) for j in range(6)] for i in range(6)] + + G0 = mat_coeff(GP_ser, 0) + zV = [iv_frac(q) for q in zv_exact] + cols = [iv_solve6(G0, [iv.mpf(1) if l == k else iv.mpf(0) for l in range(6)]) for k in range(6)] + Ginv0 = [[cols[k][i] for k in range(6)] for i in range(6)] + res = max(ub(sum((G0[i][k] * Ginv0[k][j] for k in range(6)), iv.mpf(0)) - (1 if i == j else 0)) + for i in range(6) for j in range(6)) + ck(res < 1e-70, 'CB3 G0 inverse residual') + emit('CB3 G0 inverse enclosure, residual %s' % mp.nstr(res, 3)) + + Gj = [mat_coeff(GP_ser, j) for j in range(NK + 1)] + GYj = [mat_coeff(GY_ser, j) for j in range(NK + 1)] + Aj = [mat_coeff(GYV_ser, j) for j in range(NK + 1)] # A = GYV, the soft-pin cross Gram + A_m1 = mat_coeff(GYV_ser, -1) + GY_m2 = mat_coeff(GY_ser, -2) + GY_m1 = mat_coeff(GY_ser, -1) + + W = [None] * (NK + 1) # Taylor coefficients of w = G^-1 zV + H = [None] * (NK + 1) # Taylor coefficients of G^-1 + W[0] = matvec(Ginv0, zV) + H[0] = Ginv0 + for j in range(1, NK + 1): + acc = matzero() + for i in range(1, j + 1): + acc = matadd(acc, matmul(Gj[i], H[j - i])) + H[j] = matneg(matmul(Ginv0, acc)) + W[j] = matvec(Ginv0, [-sum((Gj[i][a][l] * W[j - i][l] for i in range(1, j + 1) for l in range(6)), iv.mpf(0)) + for a in range(6)]) + + def Ac(i): + return A_m1 if i == -1 else Aj[i] + + t1 = time.time() + # mu = A w: X[k] = sum_{i+l=k} A_i W_l, i in [-1, NK], l in [0, NK]; mu_k = X[k] for k <= J + X = {} + for k in range(-1, 2 * NK + 1): + acc = [iv.mpf(0)] * 6 + for i in range(max(-1, k - NK), min(k, NK) + 1): + v = matvec(Ac(i), W[k - i]) + acc = [acc[a] + v[a] for a in range(6)] + X[k] = acc + mu_ser = [X[k] for k in range(J + 1)] + # Sigma = GY - A H A^T: AH[k] = sum_{a+b=k} A_a H_b; S[k] = sum_{m+c=k} AH_m A_c^T + AH = {} + for k in range(-1, 2 * NK + 1): + acc = matzero() + for a in range(max(-1, k - NK), min(k, NK) + 1): + acc = matadd(acc, matmul(Ac(a), H[k - a])) + AH[k] = acc + AcT = {i: transpose(Ac(i)) for i in range(-1, NK + 1)} + D = {} + for k in range(-2, 3 * NK + 1): + acc = matzero() + for c in range(max(-1, k - 2 * NK), min(k + 1, NK) + 1): + acc = matadd(acc, matmul(AH[k - c], AcT[c])) + gy = GY_m2 if k == -2 else GY_m1 if k == -1 else GYj[k] if k <= NK else None + D[k] = [[(gy[x][y] if gy is not None else 0) - acc[x][y] for y in range(6)] for x in range(6)] + Sig_ser = [D[k] for k in range(J + 1)] + # residual polynomials of the truncated inverse: zV - G_P W_P and I - G_P H_P (orders 0..2NK) + RW = {} + RH = {} + for k in range(0, 2 * NK + 1): + accv = [iv.mpf(0)] * 6 + accm = matzero() + for i in range(max(0, k - NK), min(k, NK) + 1): + v = matvec(Gj[i], W[k - i]) + accv = [accv[a] + v[a] for a in range(6)] + accm = matadd(accm, matmul(Gj[i], H[k - i])) + RW[k] = [(zV[a] if k == 0 else 0) - accv[a] for a in range(6)] + RH[k] = [[(1 if (k == 0 and x == y) else 0) - accm[x][y] for y in range(6)] for x in range(6)] + note('CB3 product and residual polynomials in %.1fs' % (time.time() - t1)) + emit('CB3 series division done: mu, Sigma to order %d; product polynomials to order %d' % (J, 3 * NK)) + + LR = [[1, 0, 0, 0, 0, 0], [-1, 1, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1]] + mu2_ser = [[sum((iv.mpf(LR[i][k]) * mu_ser[j][k] for k in range(6)), iv.mpf(0)) for i in range(6)] for j in range(J + 1)] + Sig2_ser = [[[sum((iv.mpf(LR[i][k]) * Sig_ser[j][k][l] * iv.mpf(LR[x][l]) for k in range(6) for l in range(6)), iv.mpf(0)) + for x in range(6)] for i in range(6)] for j in range(J + 1)] + + # the exact lemma CB2s says these vanish; the enclosures must agree + for a in range(6): + ck(ub(X[-1][a]) < 1e-55, 'CB3 mu pole z^-1 enclosure contains only tiny values, coord %d' % a) + for x in range(6): + for y in range(6): + ck(ub(D[-2][x][y]) < 1e-50 and ub(D[-1][x][y]) < 1e-50, 'CB3 Sigma pole enclosures tiny %d %d' % (x, y)) + # E V^(0) = Y^(-1) at the Gram level: A_{-1} = E G0 (the matrix E is exact, from CB2s) + for a in range(6): + for k in range(6): + ev = sum((iv_frac(E_sel[a][l]) * G0[l][k] for l in range(6)), iv.mpf(0)) + ck(ub(A_m1[a][k] - ev) < 1e-80, 'CB3 A_{-1} = E G0 entry %d %d' % (a, k)) + emit('CB3 poles of mu (z^-1) and Sigma (z^-2, z^-1): exact zeros by CB2s; enclosures agree') + ck(ub(mu_ser[0][0] + mpf(6) / 5) < 1e-60, 'CB3 mu_qM(0) = -6/5') + ck(ub(mu_ser[0][1] + mpf(6) / 5) < 1e-60, 'CB3 mu_qS(0) = -6/5') + ck(ub(mu_ser[0][2] + 1) < 1e-60, 'CB3 mu_aM(0) = -1 (window identity)') + ck(ub(mu_ser[0][3] - 1) < 1e-60, 'CB3 mu_aS(0) = +1') + ck(ub(Sig_ser[0][2][2]) < 1e-60, 'CB3 Var(aM)(0) = 0') + ck(ub(Sig_ser[0][0][0] - 2) < 1e-55, 'CB3 Var(qM)(0) = a4 - a2^2 = 2') + emit('CB3 PASS: limit identities mu(0), Sigma(0) match the independent blow-up analysis') + + # ------------------------------------------------------------------ CB4: majorants + # |K1^{(n)}(w)| <= A_n(rho) for complex |w| <= rho: + # A_n = (1/Z) sum_{k in (pi/12)Z} w(k) |k|^n e^{|k| rho}; the lattice holds |j| <= Jl, and for + # |k| >= k0 = (Jl+1) pi/12 >= 4 rho, w(k) e^{|k| rho} <= e^{-k^2/4}; consecutive terms + # |k|^n e^{-k^2/4} then shrink by at most q = e^{-(2 Jl + 3) h^2 / 4} (1 + 1/(Jl+1))^n < 1. + hstep, Jl, ks_iv, ws_iv, Z1_iv = ce._LATT_IV + + def abs_sum_A(n, rho_mpf): + rho = iv.mpf(rho_mpf) + tot = iv.mpf(0) + for kk, ww in zip(ks_iv, ws_iv): + ak = iv.mpf(ub(kk)) + tot += ww * ak ** n * iv.exp(ak * rho) + j0 = Jl + 1 + k0 = j0 * hstep + ck(k0.a >= 4 * rho_mpf, 'CB4 kernel tail cut k0 >= 4 rho') + q = iv.exp(-(2 * j0 + 1) * hstep * hstep / 4) * (1 + iv.mpf(1) / j0) ** n + ck(q.b < 1, 'CB4 kernel tail ratio < 1, n=%d' % n) + tail = 2 * iv.exp(-k0 * k0 / 4) * k0 ** n / (1 - q) + return (tot + iv.mpf((0, up(tail)))) / Z1_iv + + A_rho = [abs_sum_A(n, RHO) for n in range(9)] + + def entry_max2(FA, FB): + tot = iv.mpf(0) + for (xa, aa, ca) in FA: + for (xb, ab, cb) in FB: + n1 = aa[0] + ab[0] + n2 = aa[1] + ab[1] + ma = sum((iv.mpf(ub(c)) * iv.mpf(RHO) ** j for j, c in ca.items()), iv.mpf(0)) + mb = sum((iv.mpf(ub(c)) * iv.mpf(RHO) ** j for j, c in cb.items()), iv.mpf(0)) + tot += ma * mb * up(A_rho[n1]) * up(A_rho[n2]) + return up(tot) + + # sup of each Gram block over |z| = rho; Laurent coefficients obey |c_j| <= M rho^-j + M_GP = max(entry_max2(PinsV[i], PinsV[j]) for i in range(6) for j in range(6)) + M_GY = max(entry_max2(SoftY[i], SoftY[j]) for i in range(6) for j in range(6)) + M_GYV = max(entry_max2(SoftY[i], PinsV[j]) for i in range(6) for j in range(6)) + emit('CB4 Gram block maxima on |z|=%s: GP %s GY %s GYV %s' % (mp.nstr(RHO, 3), mp.nstr(M_GP, 6), mp.nstr(M_GY, 6), mp.nstr(M_GYV, 6))) + + def matmax(M): + if MUT == 'signed_maxima': + return max(mpf(M[i][j].b) for i in range(6) for j in range(6)) + return max(ub(M[i][j]) for i in range(6) for j in range(6)) + + def vecmax(v): + if MUT == 'signed_maxima': + return max(mpf(x.b) for x in v) + return max(ub(x) for x in v) + + G_m = [matmax(Gj[j]) for j in range(NK + 1)] + GY_m = [matmax(GYj[j]) for j in range(NK + 1)] + A_m = [matmax(Aj[j]) for j in range(NK + 1)] + H_m = [matmax(H[j]) for j in range(NK + 1)] + W_m = [vecmax(W[j]) for j in range(NK + 1)] + Ginv0_m = matmax(Ginv0) + A1_m = matmax(A_m1) + # B1 guard: every majorant dominates every entry of the matrix or vector it majorizes + dom = all(G_m[j] >= ub(Gj[j][a][b]) and GY_m[j] >= ub(GYj[j][a][b]) and A_m[j] >= ub(Aj[j][a][b]) + and H_m[j] >= ub(H[j][a][b]) for j in range(NK + 1) for a in range(6) for b in range(6)) + dom = dom and all(W_m[j] >= ub(W[j][a]) for j in range(NK + 1) for a in range(6)) + dom = dom and all(Ginv0_m >= ub(Ginv0[a][b]) and A1_m >= ub(A_m1[a][b]) for a in range(6) for b in range(6)) + ck(dom, 'CB4 majorants dominate every |entry| (B1)') + + RHOi = iv.mpf(RHO) + QRHO = up(1 / RHOi) # upper bound on 1/rho: the Cauchy ratio of every Gram-block majorant + # sanity check on the Cauchy constants: every computed Laurent coefficient obeys |c_j| <= M rho^-j + cauchy = all(G_m[j] <= up(iv.mpf(M_GP) * RHOi ** (-j)) and GY_m[j] <= up(iv.mpf(M_GY) * RHOi ** (-j)) + and A_m[j] <= up(iv.mpf(M_GYV) * RHOi ** (-j)) for j in range(NK + 1)) + cauchy = cauchy and A1_m <= up(iv.mpf(M_GYV) * RHOi) and matmax(GY_m1) <= up(iv.mpf(M_GY) * RHOi) \ + and matmax(GY_m2) <= up(iv.mpf(M_GY) * RHOi ** 2) + ck(cauchy, 'CB4 computed Laurent coefficients obey the Cauchy bounds |c_j| <= M rho^-j') + + def theta(R): + # 36 ||Ginv0|| sum_{i>=1} g_i R^-i, with g_i = G_m[i] (i <= J) and the Cauchy bound beyond + Ri = iv.mpf(R) + x = 1 / (RHOi * Ri) + s = sum((iv.mpf(G_m[i]) / Ri ** i for i in range(1, J + 1)), iv.mpf(0)) + iv.mpf(M_GP) * x ** (J + 1) / (1 - x) + return 36 * iv.mpf(Ginv0_m) * s + + # majorant growth R: |H_j|, |W_j| <= C R^j for all j once theta(R) <= 1/2 (induction on the + # recursion H_j = -Ginv0 sum_i G_i H_{j-i}, base j <= NK from the computed coefficients) + R = 1 / RHO * mpf('1.001') + for _ in range(600): + if up(theta(R)) > 0.5: + R *= mpf('1.02') + else: + break + ck(up(theta(R)) <= 0.5, 'CB4 majorant growth R found') + for _ in range(600): + Rl = R / mpf('1.02') + if (RHOi * iv.mpf(Rl)).a > 1 and up(theta(Rl)) <= 0.5: + R = Rl + else: + break + ck((RHOi * iv.mpf(R)).a > 1 and up(theta(R)) <= 0.5, 'CB4 rho R > 1 and theta(R) <= 1/2 in interval arithmetic') + emit('CB4 majorant growth R = %s (rho R = %s, theta = %s)' % (mp.nstr(R, 8), mp.nstr(R * RHO, 4), mp.nstr(up(theta(R)), 4))) + Ri = iv.mpf(R) + C_W = max([mpf(1)] + [up(iv.mpf(W_m[j]) / Ri ** j) for j in range(NK + 1)]) + C_H = max([mpf(1)] + [up(iv.mpf(H_m[j]) / Ri ** j) for j in range(NK + 1)]) + + class Maj: + """Majorant of a nonnegative sequence: c[j] for j <= N, then C q^j.""" + + def __init__(self, coeffs, C, q): + self.c = [iv.mpf(x) for x in coeffs] + self.N = len(coeffs) - 1 + self.C = iv.mpf(C) + self.q = iv.mpf(q) + self._t = None + + def coef(self, j): + return self.c[j] if j <= self.N else self.C * self.q ** j + + def tails(self, t): + # T[m] = sum_{j >= m} a_j t^j for m = 0 .. N+1 (suffix sums of nonnegative terms) + if self._t is not None and self._t[0] is t: + return self._t[1] + qt = self.q * t + ck(qt.b < 1, 'CB4 majorant series converges (q t < 1)') + T = [None] * (self.N + 2) + T[self.N + 1] = self.C * qt ** (self.N + 1) / (1 - qt) + for j in range(self.N, -1, -1): + T[j] = T[j + 1] + self.c[j] * t ** j + self._t = (t, T) + return T + + def tail_from(self, m, t): + m = max(m, 0) + if m <= self.N + 1: + return self.tails(t)[m] + qt = self.q * t + ck(qt.b < 1, 'CB4 majorant series converges (q t < 1)') + return self.C * qt ** m / (1 - qt) + + def full(self, t): + return self.tail_from(0, t) + + def P2(Aq, Bq, m, t): + # sum over i + l >= m of a_i b_l t^{i+l} + s = sum((Aq.coef(i) * t ** i * Bq.tail_from(m - i, t) for i in range(0, max(m, 0))), iv.mpf(0)) + return s + Aq.tail_from(m, t) * Bq.full(t) + + def P3(Aq, Bq, Cq, m, t): + # sum over a + b + c >= m of a_a b_b c_c t^{a+b+c} + s = sum((Aq.coef(a) * t ** a * P2(Bq, Cq, m - a, t) for a in range(0, max(m, 0))), iv.mpf(0)) + return s + Aq.tail_from(m, t) * Bq.full(t) * Cq.full(t) + + majW = Maj(W_m, C_W, R) + majH = Maj(H_m, C_H, R) + majA = Maj(A_m, M_GYV, QRHO) + majGY = Maj(GY_m, M_GY, QRHO) + R_STAR = F(1, 16384) if MUT != 'wide_first_cell' else F(1, 400) + ck(up(Ri * iv_frac(R_STAR)) < mpf('0.2'), 'CB4 R-majorant converges on the first cell (0, %s]' % frac_label(R_STAR)) + tS = iv.mpf(up(iv_frac(R_STAR))) + + def mu_tail(t): + # entries of sum_{k > J} mu_k t^k: mu_k = sum_{i=-1}^{k} A_i W_{k-i} + main = 6 * P2(majA, majW, J + 1, t) + pole = 6 * iv.mpf(A1_m) * majW.tail_from(J + 2, t) / t + return up(main + pole) + + def Sig_tail(t): + # entries of sum_{k > J} Sigma_k t^k: Sigma_k = GY_k - sum_{a+b+c=k} A_a H_b A_c^T + gy = majGY.tail_from(J + 1, t) + prod = 36 * P3(majA, majH, majA, J + 1, t) + pole1 = 2 * 36 * iv.mpf(A1_m) * P2(majH, majA, J + 2, t) / t + pole2 = 36 * iv.mpf(A1_m) ** 2 * majH.tail_from(J + 3, t) / (t * t) + return up(gy + prod + pole1 + pole2) + + ST = Sig_tail(tS) + MT = mu_tail(tS) + ck(ST < 1e-5 and ST >= 0, 'CB4 Sigma tail small on the first cell') + ck(MT < 1e-5 and MT >= 0, 'CB4 mu tail small on the first cell') + emit('CB4 PASS: first cell (0, %s]: R t = %s, Sigma tail %s, mu tail %s' % ( + frac_label(R_STAR), mp.nstr(up(Ri * tS), 4), mp.nstr(ST, 3), mp.nstr(MT, 3))) + + def gram_tail(M, r_up): + x = iv.mpf(r_up) / RHOi + ck(x.b < 1, 'CB6 Gram tail needs r < rho') + return up(iv.mpf(M) * x ** (NK + 1) / (1 - x)) + + # ------------------------------------------------------ CB5: truncated moments and Wick + def trunc_mom_iv(mu, var, a, b, kmax): + # E[X^k 1_{a <= X <= b}], X ~ N(mu, var), every argument an interval + ck(var.a > 0, 'CB5 truncated-moment variance enclosure positive') + s = iv.sqrt(var) + al = (a - mu) / s + be = (b - mu) / s + I = [iv.mpf(0)] * (kmax + 1) + I[0] = iv_Phi(be) - iv_Phi(al) + I[1] = iv_phi(al) - iv_phi(be) + for e in range(2, kmax + 1): + I[e] = al ** (e - 1) * iv_phi(al) - be ** (e - 1) * iv_phi(be) + (e - 1) * I[e - 2] + out = [] + for k in range(kmax + 1): + tot = iv.mpf(0) + for j in range(k + 1): + tot += comb(k, j) * s ** j * I[j] * mu ** (k - j) + out.append(tot) + return out + + def make_wick_iv(mu, Sig): + @lru_cache(maxsize=None) + def wick(state): + if len(state) == 0: + return iv.mpf(1) + if len(state) == 1: + return mu[state[0]] + i = state[0] + rest = state[1:] + tot = mu[i] * wick(rest) + for k, j in enumerate(rest): + tot += Sig[i][j] * wick(rest[:k] + rest[k + 1:]) + return tot + return wick + + def poly_expect_iv(P, wick): + tot = iv.mpf(0) + for e, c in P.items(): + st = [] + for i, pp in enumerate(e): + st += [i] * pp + tot += c * wick(tuple(st)) + return tot + + def p_affine(idx): + P = defaultdict(lambda: iv.mpf(0)) + P[tuple(1 if k == idx else 0 for k in range(6))] = iv.mpf(1) + return P + + def p_mul(A, B): + C = defaultdict(lambda: iv.mpf(0)) + for ea, ca in A.items(): + for eb, cb in B.items(): + C[tuple(x + y for x, y in zip(ea, eb))] += ca * cb + return C + + def p_scale(A, s): + C = defaultdict(lambda: iv.mpf(0)) + for e, c in A.items(): + C[e] = c * s + return C + + def p_add(A, B): + C = defaultdict(lambda: iv.mpf(0)) + for e, c in A.items(): + C[e] += c + for e, c in B.items(): + C[e] += c + return C + + def tp_add(A, B): + n = max(len(A), len(B)) + return [(A[i] if i < len(A) else iv.mpf(0)) + (B[i] if i < len(B) else iv.mpf(0)) for i in range(n)] + + def tp_mul(A, B): + C = [iv.mpf(0)] * (len(A) + len(B) - 1) + for i, a in enumerate(A): + for j, b in enumerate(B): + C[i + j] += a * b + return C + + def tp_scale(A, s): + return [s * a for a in A] + + def conditional_wick_iv(mu, Sig, cond): + i0 = cond + s00 = Sig[i0][i0] + means = [[mu[i] - Sig[i][i0] / s00 * mu[i0], Sig[i][i0] / s00] for i in range(6)] + Sc = [[Sig[i][j] - Sig[i][i0] * Sig[i0][j] / s00 for j in range(6)] for i in range(6)] + + @lru_cache(maxsize=None) + def wick(state): + if len(state) == 0: + return (iv.mpf(1),) + if len(state) == 1: + return tuple(means[state[0]]) + i = state[0] + rest = state[1:] + tot = tp_mul(list(means[i]), list(wick(rest))) + for k, j in enumerate(rest): + tot = tp_add(tot, tp_scale(list(wick(rest[:k] + rest[k + 1:])), Sc[i][j])) + return tuple(tot) + return wick + + def cond_expect_poly_iv(P, cw): + tot = (iv.mpf(0),) + for e, c in P.items(): + st = [] + for i, pp in enumerate(e): + st += [i] * pp + tot = tp_add(list(tot), tp_scale(list(cw(tuple(st))), c)) + return tot + + # ---------------------------------------------------------- CB6: law on each cell + def horner(coeffs, t_iv): + v = iv.mpf(0) + for c in reversed(coeffs): + v = v * t_iv + c + return v + + def ser_eval(coeffs, t_iv, tail): + return horner(coeffs, t_iv) + iv.mpf((-tail, tail)) + + def mat_poly_eval(mats, t_iv, shift=0): + # sum_j mats[j] t^(j + shift), entrywise + out = [[horner([m[a][b] for m in mats], t_iv) for b in range(6)] for a in range(6)] + if shift: + f = t_iv ** shift + out = [[out[a][b] * f for b in range(6)] for a in range(6)] + return out + + def vec_poly_eval(vecs, t_iv, shift=0): + out = [horner([v[a] for v in vecs], t_iv) for a in range(6)] + if shift: + f = t_iv ** shift + out = [x * f for x in out] + return out + + def frob_up(M): + return up(iv.sqrt(sum((iv.mpf(ub(M[i][j])) ** 2 for i in range(6) for j in range(6)), iv.mpf(0)))) + + def l2_up(v): + return up(iv.sqrt(sum((iv.mpf(ub(x)) ** 2 for x in v), iv.mpf(0)))) + + def mp_inverse6(GG): + A = [GG[i][:] + [mpf(1) if l == i else mpf(0) for l in range(6)] for i in range(6)] + for c in range(6): + p = max(range(c, 6), key=lambda rr: abs(A[rr][c])) + A[c], A[p] = A[p], A[c] + for rr in range(6): + if rr != c: + f = A[rr][c] / A[c][c] + A[rr] = [A[rr][j] - f * A[c][j] for j in range(12)] + return [[A[i][6 + j] / A[i][i] for j in range(6)] for i in range(6)] + + Xhi = [X[k] for k in range(J + 1, 2 * NK + 1)] + Dhi = [D[k] for k in range(J + 1, 3 * NK + 1)] + RWl = [RW[k] for k in range(0, 2 * NK + 1)] + RHl = [RH[k] for k in range(0, 2 * NK + 1)] + SQRT6 = iv.sqrt(iv.mpf(6)) + + def neumann_law(c1, c2): + """Certified (mu, Sigma) on the cell [c1, c2] (exact rationals, c1 > 0). + Heavy part: the order-J Taylor polynomials, interval Horner. Correction (B2): exact identities + G (w - W_P) = (zV - G_P W_P) - (G - G_P) W_P, + G (G^-1 - H_P) = (I - G_P H_P) - (G - G_P) H_P, + mu = A_P W_P + (A - A_P) W_P + A (w - W_P), + Sigma = [GY_P - A_P H_P A_P^T] + (GY - GY_P) - [cross terms in A - A_P] - A (G^-1 - H_P) A^T, + where A_P W_P and GY_P - A_P H_P A_P^T are explicit Laurent polynomials whose orders 0..J are + the heavy part; their pole and high orders, the residual polynomials and the Cauchy tails of + the Gram blocks (from order J+3) are bounded on the cell.""" + lo = iv_frac(c1) + hi = iv_frac(c2) + t_iv = iv.mpf((mpf(lo.a), mpf(hi.b))) + r_lo = mpf(lo.a) + r_hi = mpf(hi.b) + ck(r_lo > 0, 'CB6 Neumann cell away from 0') + tailG = gram_tail(M_GP, r_hi) + tailA = gram_tail(M_GYV, r_hi) + tailY = gram_tail(M_GY, r_hi) + GP_iv = mat_poly_eval(Gj, t_iv) + GP_iv = [[GP_iv[a][b] + iv.mpf((-tailG, tailG)) for b in range(6)] for a in range(6)] + mid = [[(mpf(GP_iv[i][j].a) + mpf(GP_iv[i][j].b)) / 2 for j in range(6)] for i in range(6)] + Bm = mp_inverse6(mid) if MUT != 'no_preconditioner' else [[mpf(1) if i == j else mpf(0) for j in range(6)] for i in range(6)] + Biv = [[iv.mpf(Bm[i][j]) for j in range(6)] for i in range(6)] + Emat = [[(1 if i == j else 0) - sum((Biv[i][k] * GP_iv[k][j] for k in range(6)), iv.mpf(0)) for j in range(6)] for i in range(6)] + rhoE = frob_up(Emat) + ck(rhoE < 0.5, 'CB6 Neumann contraction on [%s,%s]' % (frac_label(c1), frac_label(c2))) + BF = frob_up(Biv) + Gnorm = up(iv.mpf(BF) / (1 - iv.mpf(rhoE))) # ||G(t)^-1||_2 for every t in the cell + if MUT == 'old_remainder': + # the frozen correction, kept only as a negative control + GinvPnorm = 6 * sum(H_m[j] * r_hi ** j for j in range(len(H_m))) + Gn = BF + 6 * BF * rhoE / (1 - rhoE) + dGinv = Gn * 6 * tailG * GinvPnorm + APo = mat_poly_eval(Aj, t_iv) + APo = [[APo[a][b] + A_m1[a][b] / t_iv for b in range(6)] for a in range(6)] + GYVmax = max(ub(APo[a][b]) for a in range(6) for b in range(6)) + zVmax = max(ub(x) for x in zV) + mpf('1e-50') + T_mu = 6 * GYVmax * dGinv * zVmax + 6 * tailA * Gn * zVmax + T_Sig = tailY + 36 * GYVmax ** 2 * dGinv + 2 * 36 * tailA * Gn * GYVmax + 36 * tailA ** 2 * Gn + else: + WP = vec_poly_eval(W, t_iv) + WPn = l2_up(WP) + WPmax = max(ub(x) for x in WP) + rWn = l2_up(vec_poly_eval(RWl, t_iv)) + Ew = up(iv.mpf(Gnorm) * (iv.mpf(rWn) + 6 * iv.mpf(tailG) * iv.mpf(WPn))) + AP = mat_poly_eval(Aj, t_iv) + AP = [[AP[a][b] + A_m1[a][b] / t_iv for b in range(6)] for a in range(6)] + AmaxP = max(ub(AP[a][b]) for a in range(6) for b in range(6)) + Amax = up(iv.mpf(AmaxP) + iv.mpf(tailA)) + ext_mu = max(ub(x) for x in vec_poly_eval(Xhi, t_iv, shift=J + 1)) + pole_mu = up(iv.mpf(max(ub(x) for x in X[-1])) / iv.mpf(r_lo)) + T_mu = up(iv.mpf(ext_mu) + iv.mpf(pole_mu) + 6 * iv.mpf(tailA) * iv.mpf(WPmax) + SQRT6 * iv.mpf(Amax) * iv.mpf(Ew)) + HP = mat_poly_eval(H, t_iv) + HPF = frob_up(HP) + RHF = frob_up(mat_poly_eval(RHl, t_iv)) + EH = up(iv.mpf(Gnorm) * (iv.mpf(RHF) + 6 * iv.mpf(tailG) * iv.mpf(HPF))) + Dext = mat_poly_eval(Dhi, t_iv, shift=J + 1) + ext_S = max(ub(Dext[a][b]) for a in range(6) for b in range(6)) + rl = iv.mpf(r_lo) + pole_S = up(iv.mpf(max(ub(D[-2][a][b]) for a in range(6) for b in range(6))) / (rl * rl) + + iv.mpf(max(ub(D[-1][a][b]) for a in range(6) for b in range(6))) / rl) + tA = iv.mpf(tailA) + T_Sig = up(iv.mpf(ext_S) + iv.mpf(pole_S) + iv.mpf(tailY) + 12 * tA * iv.mpf(HPF) * iv.mpf(AmaxP) + + 6 * tA * tA * iv.mpf(HPF) + 6 * iv.mpf(Amax) ** 2 * iv.mpf(EH)) + ck(T_mu < 1e-2 and T_Sig < 1e-2, 'CB6 law correction small on [%s,%s]' % (frac_label(c1), frac_label(c2))) + mu = [ser_eval([mu_ser[j][a] for j in range(J + 1)], t_iv, T_mu) for a in range(6)] + Sig = [[ser_eval([Sig_ser[j][a][b] for j in range(J + 1)], t_iv, T_Sig) for b in range(6)] for a in range(6)] + return mu, Sig, T_mu, T_Sig, t_iv + + def first_cell_law(): + t_iv = iv.mpf((0, up(iv_frac(R_STAR)))) + mu = [ser_eval([mu_ser[j][a] for j in range(J + 1)], t_iv, MT) for a in range(6)] + Sig = [[ser_eval([Sig_ser[j][a][b] for j in range(J + 1)], t_iv, ST) for b in range(6)] for a in range(6)] + return mu, Sig, MT, ST, t_iv + + I_LO = mpf('-3.9') + I_HI = mpf('-0.4') + ETA = mpf('0.3') + BB = mpf('2') + AB = mpf('0.8') + EDGE_AM = mpf('-0.6') + CAP_CONS = mpf(4) # consumption cap of LPW_CONSTANT v4 + + def polys(r_iv): + aqM, aqS, aaM, aaS, abM, abS = [p_affine(i) for i in range(6)] + DM = p_add(p_mul(aaM, aqM), p_scale(p_mul(abM, abM), -r_iv)) + DS = p_add(p_mul(aaS, aqS), p_scale(p_mul(abS, abS), -r_iv)) + P = p_scale(p_mul(DM, DS), -1) + one = defaultdict(lambda: iv.mpf(0)) + one[(0,) * 6] = iv.mpf(1) + + def ubar(aa, qq, bb): + u = p_add(p_scale(p_mul(p_mul(aa, qq), p_mul(aa, qq)), iv.mpf(1) / 2), p_scale(one, iv.mpf(1) / 2)) + return p_add(u, p_scale(p_mul(bb, bb), r_iv)) + PBAR = p_mul(ubar(aaM, aqM, abM), ubar(aaS, aqS, abS)) + aqS2 = p_add(p_affine(0), p_affine(1)) + PBAR2 = p_mul(ubar(aaM, aqM, abM), ubar(aaS, aqS2, abS)) + return P, PBAR, PBAR2 + + def ceil_hi(mu, Sig, r_iv, label): + # UPPER envelope (direction analysis: section 3 of the frozen H3_BAND_CEIL.md; unchanged here). + # The 'Chernoff' pieces in the labels are Cauchy-Schwarz, E[pbar 1_A] <= sqrt(E[pbar^2] P(A)), + # with P(A) an exact Gaussian tail from iv_Phi; no exponential-moment bound is used. + # E[G_r] = E[p 1_type] <= E[p 1_{q_M<0}] + E[pbar 1_{q_M<0, D_M<=0}] + E[pbar 1_{D_S>=0}], + # {D_M<=0, q_M<0} in {a_M>-0.6} u {|b_M|>2} u {q_M in (-0.4,0)}, {D_S>=0} in {q_S>=0} u {a_S<=0}. + q = mpf_to_fraction + ck(q(-EDGE_AM) * q(-I_HI) > q(up(r_iv)) * q(BB) ** 2, 'CC6 D_M cover margin ' + label) + P, PBAR, _ = polys(r_iv) + if MUT == 'cap_raised': + PBAR = p_scale(PBAR, iv.mpf(2)) + P = p_scale(P, iv.mpf(2)) + emit('CC6 MUTATION cap_raised applied (must be caught at consumption)') + wick = make_wick_iv(tuple(mu), tuple(tuple(row) for row in Sig)) + Epbar2 = up(poly_expect_iv(p_mul(PBAR, PBAR), wick)) + Ep2 = up(poly_expect_iv(p_mul(P, P), wick)) + ck(Epbar2 >= 0 and Ep2 >= 0, 'CC6 second moments nonnegative ' + label) + + def trunc_expect(PP, cond, a, b): + cw = conditional_wick_iv(mu, Sig, cond) + cp = cond_expect_poly_iv(PP, cw) + tm = trunc_mom_iv(mu[cond], Sig[cond][cond], iv.mpf(a), iv.mpf(b), len(cp) - 1) + return sum((c * t for c, t in zip(cp, tm)), iv.mpf(0)) + + def cs_piece(second, prob): + return iv.mpf((0, up(iv.sqrt(iv.mpf(second) * iv.mpf(prob))))) + + def chern(cond, kind, edge): + # E[pbar 1_{event}] <= sqrt(E[pbar^2] P(event)), the edge beyond the mean + return cs_piece(Epbar2, event_prob_upper(mu[cond], Sig[cond][cond], kind, edge, 2, 'CC6 Chernoff ' + label)) + + Wp1 = trunc_expect(P, 0, I_LO, I_HI) + Wp2 = trunc_expect(P, 0, I_HI, mpf(0)) + Wp3 = trunc_expect(P, 0, mpf(-100), I_LO) + RESp = cs_piece(Ep2, event_prob_upper(mu[0], Sig[0][0], 'lt', mpf(-100), 20, 'CC6 left residue ' + label)) + BaM = chern(2, 'gt', EDGE_AM) # E[pbar 1_{a_M > -0.6}] + BbM = trunc_expect(PBAR, 4, mpf('-100'), -BB) + trunc_expect(PBAR, 4, BB, mpf('100')) + BbM = BbM + cs_piece(Epbar2, event_prob_upper(mu[4], Sig[4][4], 'abs', mpf(100), 20, + 'CC6 beta far tail ' + label)) # |b_M| > 100 + BqM = trunc_expect(PBAR, 0, I_HI, mpf(0)) # q_M in (-0.4, 0) + BqS = trunc_expect(PBAR, 1, mpf(0), mpf('100')) # q_S >= 0 + BqS = BqS + cs_piece(Epbar2, event_prob_upper(mu[1], Sig[1][1], 'gt', mpf(100), 20, 'CC6 q_S far tail ' + label)) + BaS = chern(3, 'lt', mpf('0')) # E[pbar 1_{a_S <= 0}] + if MUT == 'drop_fail_piece': + hi = Wp1 + Wp2 + Wp3 + RESp + else: + hi = Wp1 + Wp2 + Wp3 + RESp + BaM + BbM + BqM + BqS + BaS + emit('CC6 %-18s hi.b = %s (p-win %s %s %s; fails %s %s %s %s %s)' % ( + label, mp.nstr(up(hi), 20), mp.nstr(up(Wp1), 10), mp.nstr(up(Wp2), 10), mp.nstr(up(Wp3), 10), + mp.nstr(up(BaM), 3), mp.nstr(up(BbM), 6), mp.nstr(up(BqM), 6), mp.nstr(up(BqS), 6), mp.nstr(up(BaS), 3))) + return hi + + def band_lo(mu, Sig, Tm, Ts, t_iv, label): + # LOWER envelope: B_cert = {q_M in [I_LO, I_HI]} {|delta| <= eta} {a_M <= -0.8} {a_S >= 0.8} {|b_M| <= 2} + # lies inside the type event, p > 0 there, and pbar >= |p|, so + # E[G_r] >= E[p 1_{q_M in [I_LO, I_HI]}] - E[pbar 1_{|delta|>eta}] - E[pbar 1_{|b_M|>2}] + # - E[pbar 1_{a_M>-0.8}] - E[pbar 1_{a_S<0.8}]. + mu2 = [ser_eval([mu2_ser[j][a] for j in range(J + 1)], t_iv, up(2 * iv.mpf(Tm))) for a in range(6)] + Sig2 = [[ser_eval([Sig2_ser[j][a][b] for j in range(J + 1)], t_iv, up(4 * iv.mpf(Ts))) for b in range(6)] for a in range(6)] + r_hi = mpf_to_fraction(up(t_iv)) + q = mpf_to_fraction + ck(q(AB) * q(-I_HI) - r_hi * q(BB) ** 2 > F(1, 20), 'CB6 type margin D_M>0 ' + label) + ck(q(AB) * (q(I_HI) + q(ETA)) < F(-1, 20), 'CB6 type margin D_S<0 ' + label) + ck(q(I_HI) + q(ETA) < 0, 'CB6 q_S<0 ' + label) + P, PBAR, PBAR2 = polys(t_iv) + wick = make_wick_iv(tuple(mu), tuple(tuple(row) for row in Sig)) + cw0 = conditional_wick_iv(mu, Sig, 0) + cp0 = cond_expect_poly_iv(P, cw0) + tm0 = trunc_mom_iv(mu[0], Sig[0][0], iv.mpf(I_LO), iv.mpf(I_HI), len(cp0) - 1) + EqI = sum((c * t for c, t in zip(cp0, tm0)), iv.mpf(0)) + Ebar = poly_expect_iv(PBAR, wick) + wick2 = make_wick_iv(tuple(mu2), tuple(tuple(row) for row in Sig2)) + Ebar2 = poly_expect_iv(PBAR2, wick2) + ck(Ebar.a <= Ebar2.b and Ebar2.a <= Ebar.b, 'CB6 pbar enclosures intersect ' + label) + Epbar2 = up(poly_expect_iv(p_mul(PBAR, PBAR), wick)) + ck(Epbar2 >= 0, 'CB6 second moment nonnegative ' + label) + + def piece_outside(PP, muX, SigX, cond, a, b, Etot, oneside=None, force_chern=False): + s00 = SigX[cond][cond] + if (not force_chern) and s00.a > 1e-4: + cw = conditional_wick_iv(muX, SigX, cond) + cp = cond_expect_poly_iv(PP, cw) + tm = trunc_mom_iv(muX[cond], s00, iv.mpf(a), iv.mpf(b), len(cp) - 1) + return Etot - sum((c * t for c, t in zip(cp, tm)), iv.mpf(0)) + # exact Gaussian tail (certified Phi) and Cauchy-Schwarz: + # P(X > b) with b above the mean, P(X < a) with a below it, or P(|X| > b) with a = -b + ck(oneside is not None or a == -b, 'CB6 two-sided tail needs a window symmetric about 0') + kind, edge = {'gt': ('gt', b), 'lt': ('lt', a)}.get(oneside, ('abs', b)) + prob = event_prob_upper(muX[cond], s00, kind, edge, 2, 'CB6 Chernoff tail ' + label) + return iv.mpf((0, up(iv.sqrt(iv.mpf(Epbar2) * iv.mpf(prob))))) + + Edelta = piece_outside(PBAR2, mu2, Sig2, 1, -ETA, ETA, Ebar2, force_chern=True) + Ebeta = piece_outside(PBAR, mu, Sig, 4, -BB, BB, Ebar) + EaM = piece_outside(PBAR, mu, Sig, 2, mpf('-100'), -AB, Ebar, oneside='gt', force_chern=True) + EaS = piece_outside(PBAR, mu, Sig, 3, AB, mpf('100'), Ebar, oneside='lt', force_chern=True) + lo = EqI - Edelta - Ebeta - EaM - EaS + emit('CB6 %-18s lo.a = %s (EqI %s; pieces %s %s %s %s)' % ( + label, mp.nstr(mpf(lo.a), 20), mp.nstr(mpf(EqI.a), 10), + mp.nstr(up(Edelta), 3), mp.nstr(up(Ebeta), 6), mp.nstr(up(EaM), 3), mp.nstr(up(EaS), 3))) + return lo + + CELLS = [(F(0), R_STAR), (R_STAR, F(1, 400)), (F(1, 400), F(1, 200)), (F(1, 200), F(1, 100)), + (F(1, 100), F(1, 50)), (F(1, 50), F(3, 100)), (F(3, 100), F(1, 25)), (F(1, 25), RMAX)] + ck(all(CELLS[i][1] == CELLS[i + 1][0] for i in range(len(CELLS) - 1)) and CELLS[0][0] == 0 and CELLS[-1][1] == RMAX, + 'CB6 cells tile (0, 1/20]') + results = [] + for (c1, c2) in CELLS: + label = '(%s,%s]' % (frac_label(c1), frac_label(c2)) + if c1 == 0: + mu, Sig, Tm, Ts, t_iv = first_cell_law() + method = 'R-majorant' + else: + mu, Sig, Tm, Ts, t_iv = neumann_law(c1, c2) + method = 'Neumann' + emit('CB6 %-18s law %s: correction mu %s Sigma %s' % (label, method, mp.nstr(Tm, 3), mp.nstr(Ts, 3))) + if mode == 'ceil': + hi = ceil_hi(mu, Sig, t_iv, label) + ck(up(hi) < CAP_CONS, 'CC6 consumption cap (U <= 4) on ' + label) + bound = decimal_outward(up(hi), 'up') + emit('CERT cell=%s method=%s E[G_r] <= %s' % (label, method, bound)) + else: + lo = band_lo(mu, Sig, Tm, Ts, t_iv, label) + ck(lo.a > REQ, 'CB6 band lo > c_Z on ' + label) + bound = decimal_outward(mpf(lo.a), 'down') + emit('CERT cell=%s method=%s E[G_r] >= %s' % (label, method, bound)) + results.append((label, method, bound)) + pick = max if mode == 'ceil' else min + uniform = pick((Fraction(b) for _, _, b in results)) + rel = '<=' if mode == 'ceil' else '>=' + emit('CERT uniform=(0,1/20] E[G_r] %s %s' % (rel, decimal_outward(uniform, 'up' if mode == 'ceil' else 'down'))) + + # ----------------------------------------------- CB7: rung cross-check (not load-bearing) + # direct point law at r = 0.05 from ce.cov in mp arithmetic, with no series machinery + def direct_point_law(rpt): + Mx = (-rpt / 2, mpf(0)) + Sx = (rpt / 2, mpf(0)) + V = [[(Mx, (0, 0), mpf(1))], [(Mx, (1, 0), mpf(1))], [(Mx, (0, 1), mpf(1))], + [(Sx, (1, 0), 1 / rpt), (Mx, (1, 0), -1 / rpt)], [(Sx, (0, 1), 1 / rpt), (Mx, (0, 1), -1 / rpt)], + [(Sx, (0, 0), 1 / rpt ** 3), (Mx, (0, 0), -1 / rpt ** 3), (Sx, (1, 0), -1 / (2 * rpt ** 2)), (Mx, (1, 0), -1 / (2 * rpt ** 2))]] + Yv = [[(Mx, (0, 2), mpf(1))], [(Sx, (0, 2), mpf(1))], + [(Mx, (2, 0), 1 / rpt)], [(Sx, (2, 0), 1 / rpt)], + [(Mx, (1, 1), 1 / rpt)], [(Sx, (1, 1), 1 / rpt)]] + + def lc(FA, FB): + tot = mpf(0) + for (xa, aa, ca) in FA: + for (xb, ab, cb) in FB: + tot += ca * cb * ce.cov(aa, ab, xa[0] - xb[0], xa[1] - xb[1], mode='mpf') + return tot + GP = [[lc(V[i], V[j]) for j in range(6)] for i in range(6)] + GY = [[lc(Yv[i], Yv[j]) for j in range(6)] for i in range(6)] + GYV = [[lc(Yv[i], V[j]) for j in range(6)] for i in range(6)] + Gi = mp_inverse6(GP) + zVd = [mpf(6) / 5, 0, 0, 0, 0, mpf(-1) / 6] + csol = [sum(Gi[i][k] * zVd[k] for k in range(6)) for i in range(6)] + muD = [sum(GYV[i][k] * csol[k] for k in range(6)) for i in range(6)] + SigD = [[GY[i][j] - sum(GYV[i][k] * sum(Gi[k][l] * GYV[j][l] for l in range(6)) for k in range(6)) for j in range(6)] for i in range(6)] + return muD, SigD + + muD, SigD = direct_point_law(mpf(1) / 20) + muN, SigN, TmN, TsN, _ = neumann_law(RMAX, RMAX) + for a in range(6): + ck(muN[a].a <= muD[a] <= muN[a].b, 'CB7 Neumann mu contains the direct mu at the rung, coord %d' % a) + for b in range(6): + ck(SigN[a][b].a <= SigD[a][b] <= SigN[a][b].b, 'CB7 Neumann Sigma contains the direct Sigma, %d %d' % (a, b)) + dev_mu = max(abs(muD[a] - (mpf(muN[a].a) + mpf(muN[a].b)) / 2) for a in range(6)) + dev_S = max(abs(SigD[a][b] - (mpf(SigN[a][b].a) + mpf(SigN[a][b].b)) / 2) for a in range(6) for b in range(6)) + emit('CB7 PASS: Neumann (mu, Sigma) at r=1/20 contain the independent direct point law (cross-check, not load-bearing); ' + 'observed deviation mu %s Sigma %s against corrections %s %s' % (mp.nstr(dev_mu, 3), mp.nstr(dev_S, 3), mp.nstr(TmN, 3), mp.nstr(TsN, 3))) + + if mc: + monte_carlo(neumann_law, results, mode) + note('total runtime %.1fs' % (time.time() - t_start)) + return results + + +def monte_carlo(neumann_law, results, mode): + """NON-CERTIFYING diagnostic: Monte Carlo E[G_r] at r = 0.05 and 0.01, to stderr only.""" + try: + import numpy as np + except ImportError: + note('CC7 MC DIAGNOSTIC skipped: numpy absent (NON-CERTIFYING in any case)') + return + rng = np.random.default_rng(90210) + for c in (F(1, 20), F(1, 100)): + mu, Sig, _, _, _ = neumann_law(c, c) + m = np.array([float((mpf(x.a) + mpf(x.b)) / 2) for x in mu]) + S = np.array([[float((mpf(Sig[a][b].a) + mpf(Sig[a][b].b)) / 2) for b in range(6)] for a in range(6)]) + S = (S + S.T) / 2 + L = np.linalg.cholesky(S + 1e-14 * np.eye(6)) + Y = m[None, :] + rng.standard_normal((400000, 6)) @ L.T + qM, qS, aM, aS, bM, bS = Y.T + r = float(c) + DM = aM * qM - r * bM ** 2 + DS = aS * qS - r * bS ** 2 + G = np.where((qM < 0) & (DM > 0) & (DS < 0), -DM * DS, 0.0) + note('CC7 MC DIAGNOSTIC (NON-CERTIFYING) r=%g: E[G] = %.5f +- %.5f' % (r, G.mean(), G.std(ddof=1) / np.sqrt(len(G)))) + + +# ------------------------------------------------------------------ records and replay +def parse_certified(text): + cells, uniform = [], None + for ln in text.splitlines(): + if ln.startswith('CERT cell='): + head, rel, bound = ln.rsplit(' ', 2) + cell = head.split()[1][len('cell='):] + method = head.split()[2][len('method='):] + cells.append({'cell': cell, 'method': method, 'relation': rel, 'bound': bound}) + elif ln.startswith('CERT uniform='): + head, rel, bound = ln.rsplit(' ', 2) + uniform = {'interval': head.split()[1][len('uniform='):], 'relation': rel, 'bound': bound} + return cells, uniform + + +def file_sha256(path): + return hashlib.sha256(Path(path).read_bytes()).hexdigest() + + +MPMATH_PIN = '1.4.1' # the producer runtime every record must name; requirements.txt pins the same wheel + + +def transcript_mode(body): + """The envelope a transcript names in its first line, or None if that line is missing or unknown.""" + first = body.split('\n', 1)[0] + for mode in ('ceil', 'floor'): + if first.startswith('H3 BAND V2 %s ' % mode.upper()): + return mode + return None + + +def transcript_mutation(body): + """The mutation a transcript declares in its header, or None if the header is missing.""" + for ln in body.splitlines(): + if ln.startswith('mode dps=100 MUTATION='): + return ln[len('mode dps=100 MUTATION='):] + return None + + +def make_record(mode, body): + cells, uniform = parse_certified(body) + lines = body.splitlines() + return { + 'mutation': transcript_mutation(body), + 'all_checks_pass': bool(lines) and lines[-1] == 'ALL_CHECKS_PASS', + 'schema': 1, + 'program': PROGRAM, + 'program_sha256': file_sha256(__file__), + 'cov_exact_sha256': COV_EXACT_SHA256, + 'mode': mode, + 'environment': {'python': platform.python_version(), 'mpmath': mpmath.__version__ if mpmath else None}, + 'transcript_sha256': hashlib.sha256(body.encode()).hexdigest(), + 'cells': cells, + 'uniform': uniform, + 'scientific_effect': 'NONE', + 'review_status': 'REVIEW_REQUIRED', + } + + +def compare_record(record, expected, runtime=True): + """Return a list of mismatches between a fresh record and the reviewed one. The reviewed record must + name the pinned producer mpmath; only a recomputing replay (runtime=True) must also run it, since + verify-record recomputes nothing.""" + bad = [] + for key in ('schema', 'program', 'program_sha256', 'cov_exact_sha256', 'mode', 'mutation', 'all_checks_pass', + 'transcript_sha256', 'cells', 'uniform', 'scientific_effect', 'review_status'): + if record.get(key) != expected.get(key): + bad.append(key) + # a reviewed record is only ever one of an unmutated run that passed every check + if record.get('mutation') != 'none' or expected.get('mutation') != 'none': + bad.append('mutation is not none') + if record.get('all_checks_pass') is not True or expected.get('all_checks_pass') is not True: + bad.append('ALL_CHECKS_PASS missing') + env = expected.get('environment') + if not isinstance(env, dict) or env.get('mpmath') != MPMATH_PIN: + bad.append('environment.mpmath pin') + if runtime and record['environment'].get('mpmath') != (env.get('mpmath') if isinstance(env, dict) else None): + bad.append('environment.mpmath') + if not expected.get('cells') or not expected.get('uniform'): + bad.append('empty record') + return bad + + +def main(argv=None): + global MUT + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument('mode', choices=['structure', 'ceil', 'floor', 'verify-record']) + ap.add_argument('--h2-dir', default=str(HERE / 'h2')) + ap.add_argument('--mutation', default='', choices=[''] + sorted(MUTATIONS)) + ap.add_argument('--transcript', help='write the transcript here (ceil, floor) or read it (verify-record)') + rv = ap.add_mutually_exclusive_group() + rv.add_argument('--record', help='write the replay record (JSON) here; refused under --mutation') + rv.add_argument('--verify', help='compare the run with this reviewed record; exit 1 on any difference') + ap.add_argument('--results', help='the reviewed record (verify-record)') + ap.add_argument('--mc', action='store_true', help='NON-CERTIFYING Monte Carlo diagnostic to stderr (needs numpy)') + args = ap.parse_args(argv) + MUT = args.mutation + if args.mode == 'structure': + pins, soft, zv = mutated_tables() + try: + E = structure_lemma(pins, soft, zv) + except CheckFailed as e: + print('H3_BAND_V2 STRUCTURE FAIL: %s' % e) + return 1 + if MUT: + # a mutation the lemma cannot see is still never reported as a pass + print('H3_BAND_V2 MUTATED STRUCTURE RUN (%s): the lemma holds for these tables; NOT A CERTIFICATE' % MUT) + return 3 + rows = ['[' + ' '.join(frac_label(x) for x in row) + ']' for row in E] + print('H3_BAND_V2 STRUCTURE PASS: pins regular; soft poles z^-1 only; Y^(-1) = E V^(0) with E = %s; E zV = 0' % ' '.join(rows)) + return 0 + if args.mode == 'verify-record': + if not (args.transcript and args.results): + print('H3_BAND_V2 VERIFY_FAIL: verify-record needs --transcript and --results') + return 2 + body = Path(args.transcript).read_text() + expected = json.loads(Path(args.results).read_text()) + if not isinstance(expected, dict): + print('H3_BAND_V2 VERIFY_FAIL: the reviewed record is not a JSON object') + return 1 + # the mode comes from the transcript's own header, so a relabelled record is caught + mode = transcript_mode(body) + if mode is None: + print('H3_BAND_V2 VERIFY_FAIL: the transcript header names no envelope') + return 1 + bad = compare_record(make_record(mode, body), expected, runtime=False) + if bad: + print('H3_BAND_V2 VERIFY_FAIL: record mismatch in %s' % ', '.join(bad)) + return 1 + print('H3_BAND_V2 VERIFY PASS: transcript matches the reviewed %s record (%d cells)' % (mode, len(expected['cells']))) + return 0 + if mpmath is None: + print('H3_BAND_V2 ERROR: mpmath is required for the certificate') + return 2 + if MUT and args.record: + print('H3_BAND_V2 ERROR: --record is refused under --mutation; a mutated run is not a certificate') + return 2 + expected = json.loads(Path(args.verify).read_text()) if args.verify else None + if args.verify and not isinstance(expected, dict): + # checked before the run, so a malformed record costs no recomputation + print('H3_BAND_V2 VERIFY_FAIL: the reviewed record is not a JSON object') + return 1 + try: + run_certificate(args.mode, args.h2_dir, mc=args.mc) + except CheckFailed as e: + print('H3_BAND_V2 FAIL: %s' % e) + return 1 + emit('ALL_CHECKS_PASS') + body = '\n'.join(OUT) + '\n' + print('sha256(transcript-above) = ' + hashlib.sha256(body.encode()).hexdigest()) + if args.transcript: + Path(args.transcript).write_text(body) + record = make_record(args.mode, body) + if args.record: + Path(args.record).write_text(json.dumps(record, indent=2, sort_keys=True) + '\n') + if args.verify: + bad = compare_record(record, expected) + if bad: + print('H3_BAND_V2 VERIFY_FAIL: record mismatch in %s' % ', '.join(bad)) + return 1 + u = record['uniform'] + if MUT: + # a negative control that the run itself could not catch: never reported as a pass + print('H3_BAND_V2 MUTATED RUN (%s), NOT A CERTIFICATE: E[G_r] %s %s on (0, 1/20]' % (MUT, u['relation'], u['bound'])) + return 3 + print('H3_BAND_V2 PASS %s: E[G_r] %s %s on (0, 1/20]%s' % (args.mode, u['relation'], u['bound'], + '; matches the reviewed record' if args.verify else '')) + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/experiments/h3_band_v2/requirements.txt b/experiments/h3_band_v2/requirements.txt new file mode 100644 index 000000000..4825339af --- /dev/null +++ b/experiments/h3_band_v2/requirements.txt @@ -0,0 +1 @@ +mpmath==1.4.1 --hash=sha256:dc4f0ea2304480d4a9a48a94c1020571558ade522b44a6912efac63a586e140f diff --git a/experiments/h3_band_v2/results/ceil.json b/experiments/h3_band_v2/results/ceil.json new file mode 100644 index 000000000..7d0341f44 --- /dev/null +++ b/experiments/h3_band_v2/results/ceil.json @@ -0,0 +1,71 @@ +{ + "all_checks_pass": true, + "cells": [ + { + "bound": "3.65718683240099", + "cell": "(0,1/16384]", + "method": "R-majorant", + "relation": "<=" + }, + { + "bound": "3.66271691969455", + "cell": "(1/16384,1/400]", + "method": "Neumann", + "relation": "<=" + }, + { + "bound": "3.66435028046303", + "cell": "(1/400,1/200]", + "method": "Neumann", + "relation": "<=" + }, + { + "bound": "3.67360595202420", + "cell": "(1/200,1/100]", + "method": "Neumann", + "relation": "<=" + }, + { + "bound": "3.69805856644731", + "cell": "(1/100,1/50]", + "method": "Neumann", + "relation": "<=" + }, + { + "bound": "3.71138266685765", + "cell": "(1/50,3/100]", + "method": "Neumann", + "relation": "<=" + }, + { + "bound": "3.72414747013089", + "cell": "(3/100,1/25]", + "method": "Neumann", + "relation": "<=" + }, + { + "bound": "3.73661723530033", + "cell": "(1/25,1/20]", + "method": "Neumann", + "relation": "<=" + } + ], + "cov_exact_sha256": "f08c1c5f653f2ffd2e85d39e1112f80d15db04d15f932e5ee42db2ca3553e783", + "environment": { + "mpmath": "1.4.1", + "python": "3.11.17" + }, + "mode": "ceil", + "mutation": "none", + "program": "experiments/h3_band_v2/h3_band_v2.py", + "program_sha256": "7776d78a0df5ff04d6aabd9f0e853c1d70c69ed039007919e77b64c182860431", + "review_status": "REVIEW_REQUIRED", + "schema": 1, + "scientific_effect": "NONE", + "transcript_sha256": "d0def90db630d6b95a7be3f2ec61572f9e6acb547305cc6e5d14aff9520dc661", + "uniform": { + "bound": "3.73661723530033", + "interval": "(0,1/20]", + "relation": "<=" + } +} diff --git a/experiments/h3_band_v2/results/ceil.transcript.txt b/experiments/h3_band_v2/results/ceil.transcript.txt new file mode 100644 index 000000000..0b8f579f3 --- /dev/null +++ b/experiments/h3_band_v2/results/ceil.transcript.txt @@ -0,0 +1,42 @@ +H3 BAND V2 CEIL (successor of the frozen h3_band_ceil.py) +mode dps=100 MUTATION=none +CB2s PASS: exact Laurent structure: pins regular, soft poles z^-1 only, Y^(-1) = E V^(0), E zV = 0 +CB0 PASS: H2 library cov_exact.py sha256 f08c1c5f653f2ffd... matches the pinned digest +CB1 PASS: moments a_{2m}, m<=40, exact finite torus with certified image tails (1-a2 = 9.65e-123) +CB5 PASS: certified Phi; Gaussian tail helpers dominate the Gordon lower bounds +CB5 PASS: Gaussian event bounds dominate the Gordon lower bounds on exact point laws +CB2 PASS: Gram series to order 62; numeric pole pattern agrees with the exact lemma CB2s +CB3 G0 inverse enclosure, residual 8.26e-98 +CB3 series division done: mu, Sigma to order 60; product polynomials to order 186 +CB3 poles of mu (z^-1) and Sigma (z^-2, z^-1): exact zeros by CB2s; enclosures agree +CB3 PASS: limit identities mu(0), Sigma(0) match the independent blow-up analysis +CB4 Gram block maxima on |z|=0.11: GP 2.9579e+6 GY 343.659 GYV 19473.5 +CB4 majorant growth R = 2622.3492 (rho R = 288.5, theta = 0.4944) +CB4 PASS: first cell (0, 1/16384]: R t = 0.1601, Sigma tail 1.84e-38, mu tail 1.04e-44 +CB6 (0,1/16384] law R-majorant: correction mu 1.04e-44 Sigma 1.84e-38 +CC6 (0,1/16384] hi.b = 3.6571868324009834682 (p-win 2.664381206 0.004911289092 0.5617865876; fails 1.43e-27979210 0.045889 0.0245108 0.355708 1.16e-174863653) +CERT cell=(0,1/16384] method=R-majorant E[G_r] <= 3.65718683240099 +CB6 (1/16384,1/400] law Neumann: correction mu 7.53e-90 Sigma 5.13e-84 +CC6 (1/16384,1/400] hi.b = 3.6627169196945473051 (p-win 2.667733813 0.004936980369 0.5623366192; fails 3.01e-16679 0.046512 0.0246265 0.356571 3.48e-104075) +CERT cell=(1/16384,1/400] method=Neumann E[G_r] <= 3.66271691969455 +CB6 (1/400,1/200] law Neumann: correction mu 1.73e-72 Sigma 2.87e-68 +CC6 (1/400,1/200] hi.b = 3.6643502804630231072 (p-win 2.667821993 0.004938557765 0.5623482234; fails 1.67e-4185 0.0474018 0.0247446 0.357095 3.11e-25980) +CERT cell=(1/400,1/200] method=Neumann E[G_r] <= 3.66435028046303 +CB6 (1/200,1/100] law Neumann: correction mu 8.47e-54 Sigma 7.04e-50 +CC6 (1/200,1/100] hi.b = 3.6736059520241903406 (p-win 2.671373259 0.00496961038 0.5629204601; fails 1.29e-1050 0.05083 0.0249874 0.358525 6.45e-6476) +CERT cell=(1/200,1/100] method=Neumann E[G_r] <= 3.67360595202420 +CB6 (1/100,1/50] law Neumann: correction mu 4.46e-35 Sigma 1.85e-31 +CC6 (1/100,1/50] hi.b = 3.6980585664473056353 (p-win 2.678550597 0.005041167502 0.5640544962; fails 7.99e-265 0.0634664 0.0254915 0.361454 1.7e-1609) +CERT cell=(1/100,1/50] method=Neumann E[G_r] <= 3.69805856644731 +CB6 (1/50,3/100] law Neumann: correction mu 3.12e-24 Sigma 6.46e-21 +CC6 (1/50,3/100] hi.b = 3.7113826668576457924 (p-win 2.678123654 0.005045062003 0.5639102188; fails 1.54e-119 0.0747478 0.0260148 0.363541 3.63e-711) +CERT cell=(1/50,3/100] method=Neumann E[G_r] <= 3.71138266685765 +CB6 (3/100,1/25] law Neumann: correction mu 1.76e-16 Sigma 2.44e-13 +CC6 (3/100,1/25] hi.b = 3.7241474701308878906 (p-win 2.677231616 0.005037734936 0.56366194; fails 2.87e-68 0.0860769 0.0265615 0.365578 1.05e-397) +CERT cell=(3/100,1/25] method=Neumann E[G_r] <= 3.72414747013089 +CB6 (1/25,1/20] law Neumann: correction mu 1.97e-10 Sigma 2.04e-7 +CC6 (1/25,1/20] hi.b = 3.7366172353003289263 (p-win 2.675882529 0.005019452959 0.5633121066; fails 2.56e-44 0.0977019 0.0271328 0.367568 6.39e-253) +CERT cell=(1/25,1/20] method=Neumann E[G_r] <= 3.73661723530033 +CERT uniform=(0,1/20] E[G_r] <= 3.73661723530033 +CB7 PASS: Neumann (mu, Sigma) at r=1/20 contain the independent direct point law (cross-check, not load-bearing); observed deviation mu 1.37e-91 Sigma 2.06e-91 against corrections 1.49e-10 1.23e-7 +ALL_CHECKS_PASS diff --git a/experiments/h3_band_v2/results/floor.json b/experiments/h3_band_v2/results/floor.json new file mode 100644 index 000000000..52ce98c0c --- /dev/null +++ b/experiments/h3_band_v2/results/floor.json @@ -0,0 +1,71 @@ +{ + "all_checks_pass": true, + "cells": [ + { + "bound": "2.61755715030352", + "cell": "(0,1/16384]", + "method": "R-majorant", + "relation": ">=" + }, + { + "bound": "2.58373305711774", + "cell": "(1/16384,1/400]", + "method": "Neumann", + "relation": ">=" + }, + { + "bound": "2.58174158877405", + "cell": "(1/400,1/200]", + "method": "Neumann", + "relation": ">=" + }, + { + "bound": "2.54243388604045", + "cell": "(1/200,1/100]", + "method": "Neumann", + "relation": ">=" + }, + { + "bound": "2.45593199230241", + "cell": "(1/100,1/50]", + "method": "Neumann", + "relation": ">=" + }, + { + "bound": "2.44021698557766", + "cell": "(1/50,3/100]", + "method": "Neumann", + "relation": ">=" + }, + { + "bound": "2.41538548140069", + "cell": "(3/100,1/25]", + "method": "Neumann", + "relation": ">=" + }, + { + "bound": "2.28610115562878", + "cell": "(1/25,1/20]", + "method": "Neumann", + "relation": ">=" + } + ], + "cov_exact_sha256": "f08c1c5f653f2ffd2e85d39e1112f80d15db04d15f932e5ee42db2ca3553e783", + "environment": { + "mpmath": "1.4.1", + "python": "3.11.17" + }, + "mode": "floor", + "mutation": "none", + "program": "experiments/h3_band_v2/h3_band_v2.py", + "program_sha256": "7776d78a0df5ff04d6aabd9f0e853c1d70c69ed039007919e77b64c182860431", + "review_status": "REVIEW_REQUIRED", + "schema": 1, + "scientific_effect": "NONE", + "transcript_sha256": "c80b7a0fb4e1d08300d146aac09da8ab95341cdb471fe54a9316fadeb5afdb6f", + "uniform": { + "bound": "2.28610115562878", + "interval": "(0,1/20]", + "relation": ">=" + } +} diff --git a/experiments/h3_band_v2/results/floor.transcript.txt b/experiments/h3_band_v2/results/floor.transcript.txt new file mode 100644 index 000000000..3ff76c38b --- /dev/null +++ b/experiments/h3_band_v2/results/floor.transcript.txt @@ -0,0 +1,42 @@ +H3 BAND V2 FLOOR (successor of the frozen h3_band_floor.py) +mode dps=100 MUTATION=none +CB2s PASS: exact Laurent structure: pins regular, soft poles z^-1 only, Y^(-1) = E V^(0), E zV = 0 +CB0 PASS: H2 library cov_exact.py sha256 f08c1c5f653f2ffd... matches the pinned digest +CB1 PASS: moments a_{2m}, m<=40, exact finite torus with certified image tails (1-a2 = 9.65e-123) +CB5 PASS: certified Phi; Gaussian tail helpers dominate the Gordon lower bounds +CB5 PASS: Gaussian event bounds dominate the Gordon lower bounds on exact point laws +CB2 PASS: Gram series to order 62; numeric pole pattern agrees with the exact lemma CB2s +CB3 G0 inverse enclosure, residual 8.26e-98 +CB3 series division done: mu, Sigma to order 60; product polynomials to order 186 +CB3 poles of mu (z^-1) and Sigma (z^-2, z^-1): exact zeros by CB2s; enclosures agree +CB3 PASS: limit identities mu(0), Sigma(0) match the independent blow-up analysis +CB4 Gram block maxima on |z|=0.11: GP 2.9579e+6 GY 343.659 GYV 19473.5 +CB4 majorant growth R = 2622.3492 (rho R = 288.5, theta = 0.4944) +CB4 PASS: first cell (0, 1/16384]: R t = 0.1601, Sigma tail 1.84e-38, mu tail 1.04e-44 +CB6 (0,1/16384] law R-majorant: correction mu 1.04e-44 Sigma 1.84e-38 +CB6 (0,1/16384] lo.a = 2.6175571503035234583 (EqI 2.66420959; pieces 1.47e-1311526 0.0466524 1.04e-6994803 6.15e-6993523) +CERT cell=(0,1/16384] method=R-majorant E[G_r] >= 2.61755715030352 +CB6 (1/16384,1/400] law Neumann: correction mu 7.53e-90 Sigma 5.13e-84 +CB6 (1/16384,1/400] lo.a = 2.5837330571177465232 (EqI 2.660845426; pieces 4.61e-782 0.0771124 1.18e-4170 1.23e-4138) +CERT cell=(1/16384,1/400] method=Neumann E[G_r] >= 2.58373305711774 +CB6 (1/400,1/200] law Neumann: correction mu 1.73e-72 Sigma 2.87e-68 +CB6 (1/400,1/200] lo.a = 2.5817415887740580304 (EqI 2.660699227; pieces 1.29e-195 0.0789576 1.24e-1050 3.35e-1027) +CERT cell=(1/400,1/200] method=Neumann E[G_r] >= 2.58174158877405 +CB6 (1/200,1/100] law Neumann: correction mu 8.47e-54 Sigma 7.04e-50 +CB6 (1/200,1/100] lo.a = 2.5424338860404541865 (EqI 2.65693922; pieces 6.88e-49 0.114505 7.64e-265 3.9e-253) +CERT cell=(1/200,1/100] method=Neumann E[G_r] >= 2.54243388604045 +CB6 (1/100,1/50] law Neumann: correction mu 4.46e-35 Sigma 1.85e-31 +CB6 (1/100,1/50] lo.a = 2.4559319923024106092 (EqI 2.648927014; pieces 4.3e-12 0.192995 2.63e-67 1.86e-61) +CERT cell=(1/100,1/50] method=Neumann E[G_r] >= 2.45593199230241 +CB6 (1/50,3/100] law Neumann: correction mu 3.12e-24 Sigma 6.46e-21 +CB6 (1/50,3/100] lo.a = 2.4402169855776644195 (EqI 2.647491045; pieces 3.22e-5 0.207242 8.12e-31 1.77e-26) +CERT cell=(1/50,3/100] method=Neumann E[G_r] >= 2.44021698557766 +CB6 (3/100,1/25] law Neumann: correction mu 1.76e-16 Sigma 2.44e-13 +CB6 (3/100,1/25] lo.a = 2.4153854814006904944 (EqI 2.645587993; pieces 0.00867 0.22153 8.02e-18 2.13e-14) +CERT cell=(3/100,1/25] method=Neumann E[G_r] >= 2.41538548140069 +CB6 (1/25,1/20] law Neumann: correction mu 1.97e-10 Sigma 2.04e-7 +CB6 (1/25,1/20] lo.a = 2.2861011556287826452 (EqI 2.643210223; pieces 0.121 0.236142 1.07e-11 7.18e-9) +CERT cell=(1/25,1/20] method=Neumann E[G_r] >= 2.28610115562878 +CERT uniform=(0,1/20] E[G_r] >= 2.28610115562878 +CB7 PASS: Neumann (mu, Sigma) at r=1/20 contain the independent direct point law (cross-check, not load-bearing); observed deviation mu 1.37e-91 Sigma 2.06e-91 against corrections 1.49e-10 1.23e-7 +ALL_CHECKS_PASS diff --git a/experiments/h3_band_v2/test_h3_band_v2.py b/experiments/h3_band_v2/test_h3_band_v2.py new file mode 100644 index 000000000..8dd3678d1 --- /dev/null +++ b/experiments/h3_band_v2/test_h3_band_v2.py @@ -0,0 +1,467 @@ +"""Controls for the H3 band certificate successor v2. + +Every check goes through the program's command line, so a test exercises the same code path +as a replay. Structure and record checks need only the standard library; the replays and the +mutation controls need mpmath and are skipped without it. +""" +import hashlib +import importlib.util +import json +import subprocess +import sys +import tempfile +import unittest +from fractions import Fraction +from pathlib import Path + +HERE = Path(__file__).resolve().parent +PROG = HERE / 'h3_band_v2.py' +RESULTS = HERE / 'results' + +try: + import mpmath # noqa: F401 + HAVE_MPMATH = True +except ImportError: + HAVE_MPMATH = False + + +def run(*args): + # under python -O the program runs with -O too, so both modes exercise every control + opt = ['-O'] if sys.flags.optimize else [] + return subprocess.run([sys.executable, '-I', '-B', *opt, str(PROG), *map(str, args)], + capture_output=True, text=True, timeout=1800) + + +def run_with_mpmath(stub, *args): + """Run the program without site-packages (-S) and with sys.modules['mpmath'] preset: None hides mpmath, + a string fakes that version.""" + opt = ['-O'] if sys.flags.optimize else [] + code = ('import runpy, sys, types\n' + 'stub = %r\n' + 'if stub is None:\n' + ' sys.modules["mpmath"] = None\n' + 'else:\n' + ' fake = types.ModuleType("mpmath")\n' + ' fake.__version__, fake.mp, fake.iv, fake.mpf = stub, None, None, None\n' + ' sys.modules["mpmath"] = fake\n' + 'sys.argv = [%r] + %r\n' + 'sys.exit(runpy.run_path(%r, run_name="__main__"))\n') % (stub, str(PROG), [str(a) for a in args], str(PROG)) + return subprocess.run([sys.executable, '-I', '-B', '-S', *opt, '-c', code], capture_output=True, text=True, timeout=600) + + +def run_with_mpmath_version(version, *args): + """Run the program with the installed mpmath but its __version__ replaced, as if another release ran.""" + opt = ['-O'] if sys.flags.optimize else [] + code = ('import runpy, sys, mpmath\n' + 'mpmath.__version__ = %r\n' + 'sys.argv = [%r] + %r\n' + 'sys.exit(runpy.run_path(%r, run_name="__main__"))\n') % (version, str(PROG), [str(a) for a in args], str(PROG)) + return subprocess.run([sys.executable, '-I', '-B', *opt, '-c', code], capture_output=True, text=True, timeout=1800) + + +def load_program(): + spec = importlib.util.spec_from_file_location('h3_band_v2_under_test', PROG) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def summary(result): + lines = [ln for ln in result.stdout.splitlines() if ln.startswith('H3_BAND_V2')] + return lines[-1] if lines else '' + + +class StructureLemma(unittest.TestCase): + def test_lemma_holds_for_the_functionals_of_record(self): + r = run('structure') + self.assertEqual(r.returncode, 0, r.stdout + r.stderr) + self.assertIn('STRUCTURE PASS', summary(r)) + + def test_each_branch_of_the_lemma_catches_its_violation(self): + for mutation, message in (('structure_pin_pole', 'CB2s (a) pin V5 must have no z^-2 term'), + ('structure_soft_pole2', 'CB2s (b) soft Y2 must have no z^-2 term'), + ('structure_soft_span', 'CB2s (c) soft Y0 pole part must lie in the span'), + ('structure_window', 'CB2s (d) E zV = 0 required')): + with self.subTest(mutation=mutation): + r = run('structure', '--mutation', mutation) + self.assertEqual(r.returncode, 1, r.stdout) + self.assertIn('STRUCTURE FAIL', summary(r)) + self.assertIn(message, summary(r)) + + def test_a_mutation_the_lemma_cannot_see_is_not_a_pass(self): + # wrong_window moves a soft functional the lemma accepts; planar_moments leaves the tables alone + for mutation in ('wrong_window', 'planar_moments'): + with self.subTest(mutation=mutation): + r = run('structure', '--mutation', mutation) + self.assertEqual(r.returncode, 3, r.stdout) + self.assertIn('NOT A CERTIFICATE', summary(r)) + self.assertNotIn('STRUCTURE PASS', r.stdout) + + +class OutwardDecimals(unittest.TestCase): + def test_rounding_directions(self): + m = load_program() + x = Fraction('3.657186832400983468') + self.assertEqual(m.decimal_outward(x, 'up'), '3.65718683240099') + self.assertEqual(m.decimal_outward(x, 'down'), '3.65718683240098') + self.assertEqual(m.decimal_outward(Fraction('-2.5'), 'up', 3), '-2.50') + self.assertEqual(m.decimal_outward(Fraction(1, 3), 'down'), '0.333333333333333') + self.assertEqual(m.decimal_outward(Fraction(1, 3), 'up'), '0.333333333333334') + self.assertEqual(m.decimal_outward(Fraction(123456), 'up', 3), '124000') + + def test_rounded_value_brackets_the_input(self): + m = load_program() + for x in (Fraction(2, 7), Fraction('3.73661723530032999'), Fraction(-5, 3), Fraction(10 ** 20 + 1, 3)): + self.assertGreaterEqual(Fraction(m.decimal_outward(x, 'up')), x) + self.assertLessEqual(Fraction(m.decimal_outward(x, 'down')), x) + + +class ReviewedRecords(unittest.TestCase): + """verify-record compares a transcript with a reviewed record without recomputing.""" + + def check(self, transcript, record): + return run('verify-record', '--transcript', transcript, '--results', record) + + def test_committed_records_match_committed_transcripts(self): + for mode in ('ceil', 'floor'): + r = self.check(RESULTS / ('%s.transcript.txt' % mode), RESULTS / ('%s.json' % mode)) + self.assertEqual(r.returncode, 0, r.stdout) + self.assertIn('VERIFY PASS', summary(r)) + + def test_verify_record_does_not_depend_on_the_verifiers_mpmath(self): + # verify-record recomputes nothing, so neither a missing mpmath nor another version may matter + for stub in (None, '0.0'): + for mode in ('ceil', 'floor'): + with self.subTest(mpmath=stub, mode=mode): + r = run_with_mpmath(stub, 'verify-record', '--transcript', RESULTS / ('%s.transcript.txt' % mode), + '--results', RESULTS / ('%s.json' % mode)) + self.assertEqual(r.returncode, 0, r.stdout + r.stderr) + self.assertIn('VERIFY PASS', summary(r)) + + def test_compare_record_gates_mpmath_only_when_runtime(self): + # the replay's own call site is covered by Replays.test_replay_rejects_a_record_made_with_another_mpmath + m = load_program() + expected = json.loads((RESULTS / 'ceil.json').read_text()) + fresh = json.loads(json.dumps(expected)) + fresh['environment']['mpmath'] = '0.0' + self.assertIn('environment.mpmath', m.compare_record(fresh, expected, runtime=True)) + self.assertNotIn('environment.mpmath', m.compare_record(fresh, expected, runtime=False)) + # a reviewed record naming no producer, or with a malformed environment, is a mismatch in both paths, + # never a traceback + for change in ('no version', 'no mpmath key', 'no environment', 'environment not an object'): + malformed = json.loads(json.dumps(expected)) + if change == 'no version': + malformed['environment']['mpmath'] = None + elif change == 'no mpmath key': + del malformed['environment']['mpmath'] + elif change == 'no environment': + del malformed['environment'] + else: + malformed['environment'] = '1.4.1' + with self.subTest(change=change): + self.assertEqual(m.compare_record(expected, malformed, runtime=True), + ['environment.mpmath pin', 'environment.mpmath']) + self.assertEqual(m.compare_record(expected, malformed, runtime=False), ['environment.mpmath pin']) + + def test_committed_records_name_the_pinned_mpmath(self): + # MPMATH_PIN, requirements.txt and the committed records must agree; every check then requires the pin + pin = next(ln.split()[0].split('==')[1] for ln in (HERE / 'requirements.txt').read_text().splitlines() + if ln.startswith('mpmath==')) + self.assertEqual(load_program().MPMATH_PIN, pin) + for mode in ('ceil', 'floor'): + rec = json.loads((RESULTS / ('%s.json' % mode)).read_text()) + self.assertEqual(rec['environment']['mpmath'], pin, mode) + + def test_a_record_naming_another_mpmath_is_rejected(self): + # the reviewed record must name the pinned producer, whatever mpmath the verifier has + for change in ('another version', 'no version', 'no mpmath key', 'no environment', 'environment not an object'): + rec = json.loads((RESULTS / 'ceil.json').read_text()) + if change == 'another version': + rec['environment']['mpmath'] = '0.0' + elif change == 'no version': + rec['environment']['mpmath'] = None + elif change == 'no mpmath key': + del rec['environment']['mpmath'] + elif change == 'no environment': + del rec['environment'] + else: + rec['environment'] = '1.4.1' + with self.subTest(change=change), tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(json.dumps(rec)) + for r in (self.check(RESULTS / 'ceil.transcript.txt', p), + run_with_mpmath(None, 'verify-record', '--transcript', RESULTS / 'ceil.transcript.txt', '--results', p)): + self.assertEqual(r.returncode, 1, r.stdout + r.stderr) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in environment.mpmath pin') + + def test_a_relabelled_record_is_rejected(self): + # the mode is read from the transcript's header, not taken from the record + rec = json.loads((RESULTS / 'ceil.json').read_text()) + rec['mode'] = 'floor' + with tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(json.dumps(rec)) + r = self.check(RESULTS / 'ceil.transcript.txt', p) + self.assertEqual(r.returncode, 1, r.stdout + r.stderr) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in mode') + + def test_a_transcript_whose_header_names_no_envelope_is_rejected(self): + # even with a record made consistent with the transcript; only the first line may name the envelope + original = (RESULTS / 'ceil.transcript.txt').read_text() + for layout, text, mode in (('no header', original.split('\n', 1)[1], None), + ('no header', original.split('\n', 1)[1], 'ceil'), + ('header on line 2', 'copied transcript\n' + original, 'ceil'), + ('blank first line', '\n' + original, 'ceil'), + ('indented header', ' ' + original, 'ceil')): + rec = json.loads((RESULTS / 'ceil.json').read_text()) + rec['mode'] = mode + rec['transcript_sha256'] = hashlib.sha256(text.encode()).hexdigest() + with self.subTest(layout=layout, mode=mode), tempfile.TemporaryDirectory() as d: + t = Path(d) / 't.txt' + p = Path(d) / 'r.json' + t.write_text(text) + p.write_text(json.dumps(rec)) + r = self.check(t, p) + self.assertEqual(r.returncode, 1, r.stdout + r.stderr) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: the transcript header names no envelope') + + NOT_OBJECTS = ('[]', '[1]', 'null', '"x"', '1', 'true') # falsy and truthy, so the check must be on the type + + def test_a_record_that_is_not_an_object_is_rejected(self): + for payload in self.NOT_OBJECTS: + with self.subTest(record=payload), tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(payload) + r = self.check(RESULTS / 'ceil.transcript.txt', p) + self.assertEqual(r.returncode, 1, r.stdout + r.stderr) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: the reviewed record is not a JSON object') + + def test_replay_checks_the_record_type_before_running(self): + # the fake mpmath cannot compute, so only a check made before the run can produce this summary; + # Replays.test_replay_rejects_a_record_that_is_not_an_object_before_running repeats it with mpmath + pin = load_program().MPMATH_PIN + for payload in self.NOT_OBJECTS: + with self.subTest(record=payload), tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(payload) + r = run_with_mpmath(pin, 'ceil', '--verify', p) + self.assertEqual(r.returncode, 1, r.stdout + r.stderr) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: the reviewed record is not a JSON object') + + def test_a_changed_check_line_is_rejected(self): + # only the transcript digest binds the lines that carry no certified decimal + text = (RESULTS / 'ceil.transcript.txt').read_text() + line = next(ln for ln in text.splitlines() if ln.startswith('CB0 PASS:')) + with tempfile.TemporaryDirectory() as d: + t = Path(d) / 't.txt' + t.write_text(text.replace(line, 'CB0 SKIPPED: H2 library digest not checked')) + for r in (self.check(t, RESULTS / 'ceil.json'), + run_with_mpmath(None, 'verify-record', '--transcript', t, '--results', RESULTS / 'ceil.json')): + self.assertEqual(r.returncode, 1, r.stdout + r.stderr) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in transcript_sha256') + + def test_records_carry_review_status_and_no_scientific_effect(self): + for mode in ('ceil', 'floor'): + rec = json.loads((RESULTS / ('%s.json' % mode)).read_text()) + self.assertEqual(rec['scientific_effect'], 'NONE') + self.assertEqual(rec['review_status'], 'REVIEW_REQUIRED') + self.assertEqual(len(rec['cells']), 8) + + def test_floor_lies_below_ceiling_on_every_cell(self): + ceil = json.loads((RESULTS / 'ceil.json').read_text())['cells'] + floor = json.loads((RESULTS / 'floor.json').read_text())['cells'] + self.assertEqual([c['cell'] for c in ceil], [c['cell'] for c in floor]) + for c, f in zip(ceil, floor): + self.assertEqual((c['relation'], f['relation']), ('<=', '>=')) + self.assertLess(Fraction(f['bound']), Fraction(c['bound']), c['cell']) + + def test_a_lowered_ceiling_in_the_transcript_is_rejected(self): + text = (RESULTS / 'ceil.transcript.txt').read_text() + line = next(ln for ln in text.splitlines() if ln.startswith('CERT cell=(0,1/16384]')) + doctored = line[:-1] + str((int(line[-1]) + 9) % 10) + with tempfile.TemporaryDirectory() as d: + t = Path(d) / 't.txt' + t.write_text(text.replace(line, doctored)) + r = self.check(t, RESULTS / 'ceil.json') + self.assertEqual(r.returncode, 1, r.stdout) + self.assertIn('VERIFY_FAIL', summary(r)) + + def test_a_doctored_record_is_rejected(self): + rec = json.loads((RESULTS / 'floor.json').read_text()) + rec['cells'][-1]['bound'] = '2.30659559567154' + with tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(json.dumps(rec)) + r = self.check(RESULTS / 'floor.transcript.txt', p) + self.assertEqual(r.returncode, 1, r.stdout) + self.assertIn('cells', summary(r)) + + def doctored_pair(self, mode, edit): + # a transcript edited by `edit` and a record made consistent with it from the same program + text = edit((RESULTS / ('%s.transcript.txt' % mode)).read_text()) + rec = load_program().make_record(mode, text) + # the producer environment of the reviewed record, so only the intended defect differs + rec['environment'] = json.loads((RESULTS / ('%s.json' % mode)).read_text())['environment'] + return text, rec + + def test_a_mutated_transcript_with_a_matching_record_is_rejected(self): + text, rec = self.doctored_pair('ceil', lambda t: t.replace('MUTATION=none', 'MUTATION=drop_fail_piece')) + with tempfile.TemporaryDirectory() as d: + t = Path(d) / 't.txt' + p = Path(d) / 'r.json' + t.write_text(text) + p.write_text(json.dumps(rec)) + r = self.check(t, p) + self.assertEqual(r.returncode, 1, r.stdout) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in mutation is not none') + + def test_a_transcript_without_all_checks_pass_is_rejected(self): + text, rec = self.doctored_pair('floor', lambda t: t.replace('ALL_CHECKS_PASS\n', '')) + with tempfile.TemporaryDirectory() as d: + t = Path(d) / 't.txt' + p = Path(d) / 'r.json' + t.write_text(text) + p.write_text(json.dumps(rec)) + r = self.check(t, p) + self.assertEqual(r.returncode, 1, r.stdout) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in ALL_CHECKS_PASS missing') + + def test_record_and_verify_together_are_refused(self): + with tempfile.TemporaryDirectory() as d: + p = Path(d) / 'same.json' + r = run('ceil', '--record', p, '--verify', p) + self.assertFalse(p.exists()) + self.assertEqual(r.returncode, 2, r.stdout + r.stderr) + self.assertIn('not allowed with argument', r.stderr) + + def test_a_record_for_other_program_bytes_is_rejected(self): + rec = json.loads((RESULTS / 'ceil.json').read_text()) + rec['program_sha256'] = '0' * 64 + with tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(json.dumps(rec)) + r = self.check(RESULTS / 'ceil.transcript.txt', p) + self.assertEqual(r.returncode, 1, r.stdout) + self.assertIn('program_sha256', summary(r)) + + +@unittest.skipUnless(HAVE_MPMATH, 'mpmath is required for the certificate replays') +class Replays(unittest.TestCase): + def test_ceiling_replay_matches_the_reviewed_record(self): + r = run('ceil', '--verify', RESULTS / 'ceil.json') + self.assertEqual(r.returncode, 0, r.stdout[-2000:] + r.stderr[-2000:]) + self.assertIn('matches the reviewed record', summary(r)) + + def test_floor_replay_matches_the_reviewed_record(self): + r = run('floor', '--verify', RESULTS / 'floor.json') + self.assertEqual(r.returncode, 0, r.stdout[-2000:] + r.stderr[-2000:]) + self.assertIn('matches the reviewed record', summary(r)) + + def test_replay_rejects_a_record_made_with_another_mpmath(self): + # the record differs from the reviewed one only in the producer's mpmath version + rec = json.loads((RESULTS / 'ceil.json').read_text()) + rec['environment']['mpmath'] = '0.0' + with tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(json.dumps(rec)) + r = run('ceil', '--verify', p) + self.assertEqual(r.returncode, 1, r.stdout[-2000:]) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in environment.mpmath pin, environment.mpmath') + + def test_replay_rejects_another_running_mpmath_for_the_pinned_record(self): + # the record names the pin, so only the runtime gate can fire + r = run_with_mpmath_version('0.0', 'ceil', '--verify', RESULTS / 'ceil.json') + self.assertEqual(r.returncode, 1, r.stdout[-2000:] + r.stderr[-2000:]) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: record mismatch in environment.mpmath') + + def test_replay_rejects_a_record_that_is_not_an_object_before_running(self): + for payload in ReviewedRecords.NOT_OBJECTS: + with self.subTest(record=payload), tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + p.write_text(payload) + r = run('ceil', '--verify', p) + self.assertEqual(r.returncode, 1, r.stdout[-2000:]) + self.assertEqual(summary(r), 'H3_BAND_V2 VERIFY_FAIL: the reviewed record is not a JSON object') + self.assertNotIn('ALL_CHECKS_PASS', r.stdout) + + def test_replay_rejects_a_self_consistent_doctored_record(self): + # transcript and record doctored together, so verify-record alone would accept them; + # only recomputation exposes the change + text = (RESULTS / 'ceil.transcript.txt').read_text().replace( + 'E[G_r] <= 3.73661723530033', 'E[G_r] <= 3.23510806397398') + m = load_program() + rec = json.loads((RESULTS / 'ceil.json').read_text()) + fresh = m.make_record('ceil', text) + rec.update({k: fresh[k] for k in ('transcript_sha256', 'cells', 'uniform')}) + with tempfile.TemporaryDirectory() as d: + t = Path(d) / 't.txt' + p = Path(d) / 'r.json' + t.write_text(text) + p.write_text(json.dumps(rec)) + self.assertEqual(run('verify-record', '--transcript', t, '--results', p).returncode, 0) + r = run('ceil', '--verify', p) + self.assertEqual(r.returncode, 1, r.stdout[-2000:]) + self.assertIn('VERIFY_FAIL', summary(r)) + + +@unittest.skipUnless(HAVE_MPMATH, 'mpmath is required for the mutation controls') +class Mutations(unittest.TestCase): + # mutation -> (envelope, the check that must fail) + CAUGHT_IN_RUN = { + 'structure_pin_pole': ('ceil', 'CB2s (a) pin V5 must have no z^-2 term'), + 'structure_soft_pole2': ('floor', 'CB2s (b) soft Y2 must have no z^-2 term'), + 'structure_soft_span': ('ceil', 'CB2s (c) soft Y0 pole part must lie in the span'), + 'structure_window': ('floor', 'CB2s (d) E zV = 0 required'), + 'planar_moments': ('ceil', 'CB1d consumed a_{2m} meets the spectral enclosure'), + 'wrong_window': ('ceil', 'CB3 mu_aM(0) = -1'), + 'signed_maxima': ('ceil', 'CB4 majorants dominate every |entry| (B1)'), + 'wide_first_cell': ('ceil', 'CB4 R-majorant converges on the first cell (0, 1/400]'), + 'one_tail': ('floor', 'CB5 event bound P(|X| > t) self-test'), + 'tail_helper_one_side': ('ceil', 'CB5 two-sided tail helper dominates twice the Gordon lower bound'), + 'no_preconditioner': ('ceil', 'CB6 Neumann contraction'), + 'cap_raised': ('ceil', 'CC6 consumption cap'), + } + CAUGHT_BY_REPLAY = ('drop_fail_piece', 'old_remainder') + + def test_every_mutation_has_a_control(self): + self.assertEqual(set(load_program().MUTATIONS), set(self.CAUGHT_IN_RUN) | set(self.CAUGHT_BY_REPLAY)) + + def test_mutations_are_caught_inside_the_run(self): + for mutation, (mode, message) in self.CAUGHT_IN_RUN.items(): + with self.subTest(mutation=mutation): + r = run(mode, '--mutation', mutation) + self.assertEqual(r.returncode, 1, r.stdout[-2000:]) + self.assertIn('H3_BAND_V2 FAIL', summary(r)) + self.assertIn(message, summary(r)) + + def test_mutations_are_caught_by_the_replay(self): + # the mutation header alone already changes the transcript digest, so require the + # certified decimals themselves to differ from the reviewed record + reviewed = json.loads((RESULTS / 'ceil.json').read_text()) + for mutation in self.CAUGHT_BY_REPLAY: + with self.subTest(mutation=mutation): + r = run('ceil', '--mutation', mutation, '--verify', RESULTS / 'ceil.json') + self.assertEqual(r.returncode, 1, r.stdout[-2000:]) + self.assertIn('VERIFY_FAIL', summary(r)) + self.assertIn('cells', summary(r)) + cells, _ = load_program().parse_certified(r.stdout) + self.assertEqual(len(cells), len(reviewed['cells'])) + self.assertTrue(any(c['bound'] != e['bound'] for c, e in zip(cells, reviewed['cells'])), mutation) + + def test_a_mutated_run_is_never_reported_as_a_pass(self): + r = run('ceil', '--mutation', 'drop_fail_piece') + self.assertEqual(r.returncode, 3, r.stdout[-2000:]) + self.assertIn('NOT A CERTIFICATE', summary(r)) + self.assertFalse(any(ln.startswith('H3_BAND_V2 PASS') for ln in r.stdout.splitlines())) + + def test_record_is_refused_under_a_mutation(self): + with tempfile.TemporaryDirectory() as d: + p = Path(d) / 'r.json' + r = run('ceil', '--mutation', 'drop_fail_piece', '--record', p) + self.assertFalse(p.exists()) + self.assertEqual(r.returncode, 2, r.stdout) + self.assertIn('--record is refused under --mutation', summary(r)) + + +if __name__ == '__main__': + unittest.main()