diff --git a/research/crt-recursive-interchange-crossover/README.md b/research/crt-recursive-interchange-crossover/README.md new file mode 100644 index 0000000000..97f850c3a1 --- /dev/null +++ b/research/crt-recursive-interchange-crossover/README.md @@ -0,0 +1,141 @@ +# Crossover length for the recursive ordinary interchange in the CRT reservoir + +**Scoped proof with exact computation, conditional on a stated constant +ratio. No new κ.** The selected conditional result is unchanged. This does not +accept PR #4100 or BIT.md's theorem. + +## Question + +PR #4942 (`research/crt-interchange-width-threshold`, James Reed) shows that +BIT.md §4 takes its recursive branch exactly from binary width +`H_b = 27706412`, so every #4174 reservoir with `L < 2770641200` is +elementary. Its open question 1 asks: at what reservoir length does the +recursive bill `C V H^(1-a2)` first beat the elementary bill `c_el V H`, for a +stated ratio `rho = C/c_el`? BIT.md fixes `c = 8667/5000000` and +`a2 = c(1-c^2)`, an exact rational. + +## Result + +Write `rho_c = H_b^a2`. Then: + +1. **Crossover width.** For every `rho > 0`, the least reservoir width at which + the recursive branch runs and its bill is strictly below the elementary bill + is + + `H*(rho) = max(H_b, h(rho))`, where `h(rho)` is the least integer with `h^a2 > rho`. + + So `H*(rho) = H_b` exactly when `rho < rho_c`, and `H*(rho) = h(rho)`, + about `rho^(1/a2) ≈ rho^576.9`, when `rho >= rho_c`. +2. **Critical ratio.** `rho_c = 1.030151098331576173576873168658...` + (30 digits certified). At the branch threshold the `H^(-a2)` factor is + worth only 3.02%. Any recursive constant more than 3.02% above the + elementary constant moves the crossover beyond the branch threshold. +3. **Exact least length at the branch threshold.** Over every layout the #4174 + constructor accepts (power-of-two `d >= 51200`, `2 <= m <= ell` with + `m = ceil(log2(p+2))`, widths `w_i in {ell-1, ell}`, `G = 64(ell+m)`, + `r_max = floor(L/(100G))`, `H = r_max G`), the least `L` with + `H >= H_b` is **`L = 2770643200`**. One witness is `d = 2^28`, `ell = 11`, + `m = 2` (`p = 1` or `2`), `G = 832`, `r_max = 33301`, `H = 27706432`. + This sharpens PR4942's lower bound `100 H_b = 2770641200` by 2000 and is + the answer for every `rho < rho_c`. +4. **Ratio-bound examples.** The least lengths are exact for + `rho = 1.031` (`H* = 44561394`, `L = 4456140800`) and + `rho = 1.032` (`H* = 77956662`, `L = 7795667200`). For larger ratios the + verifier certifies `H*` exactly and reports the lower bound `L >= 100 H*`: + + | rho | H* | bits of H* | + |---|---|---| + | 1.05 | 1675568503989 | 41 | + | 1.1 | 757789515901820833519268 | 80 | + | 1.5 | exact integer, not printed | 338 | + | 2 | exact integer, not printed | 577 | + | 4 | exact integer, not printed | 1154 | + | 10 | exact integer, not printed | 1917 | + + With `rho = 2` a profitable recursive call needs a reservoir of more than + `100 · 2^576` bits. + +## Argument + +Fix `rho > 0`. Below `H_b` the wrapper runs the elementary branch (PR4942, +used as stated), so no recursive bill applies. For `H >= H_b`, +`C V H^(1-a2) < c_el V H` is equivalent to `H^a2 > rho`, because `V, H > 0`. +The map `H ↦ H^a2` is strictly increasing, so the set of widths satisfying both +conditions is `{H >= max(H_b, h(rho))}`. Also `H_b^a2 > rho` exactly when +`rho < rho_c`, which gives the regime split. + +For the length, `H = floor(L/(100G)) G` is nondecreasing in `L` at fixed +`(d, ell, m)`. At fixed `k = ell + m` (so `G = 64k`), the least `L` with +`H >= H*` is `t = 100 G ceil(H*/G)`, provided `t` lies in some +`[d(ell-1), d ell]` with admissible `ell in [ceil(k/2), k-2]` and `d = 2^s >= 2^16`. +Every integer in that interval is reachable by choosing how many widths equal +`ell`. If no admissible interval contains `t`, the least reachable value +above `t` is `d(ceil(k/2)-1)` for the smallest qualifying `d`. Since +`t >= 6400k`, the scan over `k` stops once `6400k` exceeds the best value +found. Every feasible `(s, ell, m)` is therefore covered. + +## Evidence + +`verify.py` uses only the standard library and checks these points: + +- Two logarithm engines decide `h(rho)`. Engine 1 is a rigorous fixed-point + interval logarithm (atanh series with directed rounding and an explicit tail + bound), in integer arithmetic only. Engine 2 uses Python's correctly rounded + `decimal` `ln`/`exp` to propose `h`. Engine 1 then certifies both + `h^a2 > rho` and `(h-1)^a2 < rho`. An unresolved interval aborts. +- The 30-digit bracket of `rho_c` is certified by engine 1 at both endpoints. + Just below `rho_c` the width is `H_b`, and just above it the width is `H_b+1`. +- The pruned layout search is compared with a complete brute scan over + `(s, ell, m)` on five fixtures (`H = 10752, 15360, 27706, 100000, 123457`). + The verifier also checks that the brute ranges are large enough for each + fixture. Every reported witness is checked against the domain. +- Five rejection controls: `L-1` at the witness, `m = 1`, `d = 2^15`, a ratio + just above `rho_c` at `H_b`, and a nonpositive ratio. +- The BIT.md constants and the four placement-domain phrases are read from the + pinned inputs. A fresh `RESULT.json` is compared byte for byte. + +The registered check runs in about 1.5 seconds. + +## Reproduce + +From the repository root, with CPython 3.12: + +```sh +python3.12 -B research/crt-recursive-interchange-crossover/verify.py +python3.12 -B scripts/research_package.py --package research/crt-recursive-interchange-crossover +``` + +Expected: +`PASS_CRT_RECURSIVE_INTERCHANGE_CROSSOVER rho_c=1.030151098331576173576873168658 least_L_at_branch=2770643200 ...` + +## Limitations + +- Conditional on BIT.md's bill as stated and on PR4942's branch threshold. Both + are cited, not re-proved. BIT.md gives neither `C` nor `c_el`. +- The comparison uses the per-layer bills only. It omits the 6x placement + multiplicity and the `(q^E/2^H)^2` padding from PR4944. A common factor + cancels, but padding that differs between branches would change `rho`. +- Exact least lengths are computed only when `H* <= 10^9`. Above that, only + `L >= 100 H*` is reported. +- No statement about the complete CRT bill, supplier admission or any + multiplication theorem. + +## Open next questions + +1. **What are `C` and `c_el`?** Once `rho` exceeds 1.0302, the crossover + grows like `rho^577`. A definite ratio, including padding, decides whether + any physical reservoir uses the recursive branch profitably. First step: + bound `c_el` from the elementary q-digit swap and `C` from the body word at + `w = 2`, `D = 19`. BIT.md states neither constant. +2. **Exact lengths for larger `H*`.** Replace the linear `k` scan with a + divisor search near `64 ceil(H*/64)` under the `d >= 65536` constraint. + That would give exact least lengths for `rho` in `(1.032, 1.05]` and beyond. + +## Credit + +Builds on James Reed's PR #4942 (branch threshold and the open question), +JoyboyBrian's PR #4944 (six-call reservoir bill and padding), Hiroki KUMAZAKI's +PR #4174 (constructor and layout domain) and Cinna Davis's PR #4100 (CRT +construction and BIT.md). The crossover analysis, critical ratio, layout +search and verifier were produced by Nick Petalas. Apache-2.0. Independent +human and formal proof-kernel review are unclaimed. diff --git a/research/crt-recursive-interchange-crossover/RESULT.json b/research/crt-recursive-interchange-crossover/RESULT.json new file mode 100644 index 0000000000..c491525f7c --- /dev/null +++ b/research/crt-recursive-interchange-crossover/RESULT.json @@ -0,0 +1,145 @@ +{ + "H_branch": 27706412, + "a2": "216674348961923037/125000000000000000000", + "layout_fixtures": [ + { + "H_target": 10752, + "least_L": 1075200 + }, + { + "H_target": 15360, + "least_L": 1536000 + }, + { + "H_target": 27706, + "least_L": 2777600 + }, + { + "H_target": 100000, + "least_L": 10009600 + }, + { + "H_target": 123457, + "least_L": 12352000 + } + ], + "least_L_at_branch": 2770643200, + "rho_c_floor_1e-30": "1.030151098331576173576873168658", + "rows": [ + { + "H_star": "27706412", + "H_star_bits": 25, + "binding": "branch threshold", + "least_H_with_H^a2>rho": "1", + "least_L": 2770643200, + "rho": "1/2" + }, + { + "H_star": "27706412", + "H_star_bits": 25, + "binding": "branch threshold", + "least_H_with_H^a2>rho": "2", + "least_L": 2770643200, + "rho": "1" + }, + { + "H_star": "27706412", + "H_star_bits": 25, + "binding": "branch threshold", + "least_H_with_H^a2>rho": "25458253", + "least_L": 2770643200, + "rho": "103/100" + }, + { + "H_star": "44561394", + "H_star_bits": 26, + "binding": "bill ratio", + "least_H_with_H^a2>rho": "44561394", + "least_L": 4456140800, + "rho": "1031/1000", + "witness": { + "G": 256, + "H": 44561408, + "L": 4456140800, + "d": 4294967296, + "ell": 2, + "m": 2, + "r_max": 174068 + } + }, + { + "H_star": "77956662", + "H_star_bits": 27, + "binding": "bill ratio", + "least_H_with_H^a2>rho": "77956662", + "least_L": 7795667200, + "rho": "1032/1000", + "witness": { + "G": 23488, + "H": 77956672, + "L": 7795667200, + "d": 33554432, + "ell": 233, + "m": 134, + "r_max": 3319 + } + }, + { + "H_star": "1675568503989", + "H_star_bits": 41, + "L_lower_bound_100H": "167556850398900", + "binding": "bill ratio", + "least_H_with_H^a2>rho": "1675568503989", + "rho": "21/20" + }, + { + "H_star": "757789515901820833519268", + "H_star_bits": 80, + "L_lower_bound_100H": "75778951590182083351926800", + "binding": "bill ratio", + "least_H_with_H^a2>rho": "757789515901820833519268", + "rho": "11/10" + }, + { + "H_star": null, + "H_star_bits": 338, + "L_lower_bound_100H": "100*H_star, H_star has 338 bits", + "binding": "bill ratio", + "least_H_with_H^a2>rho": null, + "rho": "3/2" + }, + { + "H_star": null, + "H_star_bits": 577, + "L_lower_bound_100H": "100*H_star, H_star has 577 bits", + "binding": "bill ratio", + "least_H_with_H^a2>rho": null, + "rho": "2" + }, + { + "H_star": null, + "H_star_bits": 1154, + "L_lower_bound_100H": "100*H_star, H_star has 1154 bits", + "binding": "bill ratio", + "least_H_with_H^a2>rho": null, + "rho": "4" + }, + { + "H_star": null, + "H_star_bits": 1917, + "L_lower_bound_100H": "100*H_star, H_star has 1917 bits", + "binding": "bill ratio", + "least_H_with_H^a2>rho": null, + "rho": "10" + } + ], + "witness_at_branch": { + "G": 832, + "H": 27706432, + "L": 2770643200, + "d": 268435456, + "ell": 11, + "m": 2, + "r_max": 33301 + } +} diff --git a/research/crt-recursive-interchange-crossover/SOURCES.json b/research/crt-recursive-interchange-crossover/SOURCES.json new file mode 100644 index 0000000000..2453b795d9 --- /dev/null +++ b/research/crt-recursive-interchange-crossover/SOURCES.json @@ -0,0 +1,30 @@ +{ + "schema_version": 1, + "base_commit": "7c49fe32c55de02c2c873655c804b16e4c3fef0a", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "d59b5035e7b4cbaa9d688a1d32dd0d6eb79c466a", + "sha256": "9b8d0eea95b3bf0dfcb7938a732a6439b87a2685dec2a5ab9e525b3b34debbfa" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "3404fbdd5e7529a0674eef26114b5071e4586cec", + "sha256": "ead6cdc9d0af6d58e8a49744d8ec8c7c15e13d2ba93c276c434bcf9b51b93a81" + }, + { + "path": "research/crt-reservoir-placement/README.md", + "usage": "input", + "git_blob": "2b6e7313b75bca33d8c72ccbd5b373d47eeeef42", + "sha256": "53bc3330c70b263784b9e28d3a059fb2e9b737683a1c5b6bad41513dbb5da89a" + }, + { + "path": "research/unconditional-batched-crt/BIT.md", + "usage": "input", + "git_blob": "2aa2544fa61259ac48b936b8c906ec9d31d813fa", + "sha256": "83a0cc082e52219a5f142becc8b5b92019e091cb3ede6c5bccbd8508993fc4da" + } + ] +} diff --git a/research/crt-recursive-interchange-crossover/TASK.json b/research/crt-recursive-interchange-crossover/TASK.json new file mode 100644 index 0000000000..f6be0c1b41 --- /dev/null +++ b/research/crt-recursive-interchange-crossover/TASK.json @@ -0,0 +1,89 @@ +{ + "schema_version": 1, + "task_id": "uniform-transfer-audit", + "priority_revision": 8, + "base_commit": "7c49fe32c55de02c2c873655c804b16e4c3fef0a", + "original_card": { + "id": "uniform-transfer-audit", + "question": "Can one of the selected record's five inherited assumptions be proved at its actual caller?", + "starting_points": [ + "docs/research/current-priorities.md", + "docs/research/research-family-status.md", + "certificates/selected-result.json", + "docs/research/pr-intake-evidence/pr5599-completion/PROOF-REVIEW.md" + ], + "methods": [ + "Construct and independently check one source-bound instance of the stated first experiment.", + "Seek an exact counterexample or target-specific complete-budget obstruction to that same proposed change." + ], + "acceptance_gate": "Identify the proposition and what would change downstream; do not equate a finite replay with a universal proof.", + "first_step": { + "command": [ + "python3", + "-B", + "scripts/research_selected_baseline.py" + ], + "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY", + "inspect": "Read revision 8, the selected review and live exact-source handoffs. Then follow this card's bounded first_experiment.", + "scope": "Source-bound orientation only, not a supplier replay or theorem proof." + }, + "first_experiment": "Choose one assumption ID (native_all_columns, owned_erasure, compiler, precision_tape, analytic_recovery), one quantified implication and exact dependency closure. Check current handoffs, including #5625 row reentry and #5707 root-prefix lowering, without treating their partial scope as the whole theorem. Prove or falsify that implication with a caller-bound argument; stop at the chosen obligation. Finite tests or root digits do not discharge an all-size statement.", + "bottleneck_link": "Judge the actual current supplier contract and complete signed target/paid bill; frozen-profile ceilings require an applicability proof.", + "allowed_variations": [ + "One exact source head, helper/cut region or quantified proposition and changed hypothesis; check live handoffs first. Different seeds are not different scopes." + ], + "scope_change_triggers": [ + "Changing source head, helper/cut domain, operation alphabet or caller boundary without recording the new scope.", + "Treating an archived proposal or local numerical pass as complete supplier admission." + ], + "overlap_terms": [ + "native_all_columns", + "owned_erasure", + "compiler", + "precision_tape", + "analytic_recovery" + ], + "existing_checks": [ + { + "scope": "PR #4100: proposed unconditional construction and finite checks, without independent acceptance.", + "evidence": "research/unconditional-batched-crt/README.md", + "remaining": "Independent proof audit of one caller-bound all-size implication; all other premises remain explicit.", + "status": "author-checked", + "source_commit": "337d04939571a637aeab54a0f58f54fde811663a", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "Selected PR5599: complete helper-subspace construction and independent caller replay; five framework assumptions remain.", + "evidence": "docs/research/pr-intake-evidence/pr5599-completion/PROOF-REVIEW.md", + "remaining": "The next experiment must identify uncovered source-specific work. Successor archive intake does not transfer this admission.", + "status": "reviewed", + "source_commit": "09dd0965711ca7eeb169057a420abab1c9052137", + "selected_record_sha256": "2a0220e866417774606e0f991e41b9c1d26b1f43b97857238cf478dbc7d330b5" + } + ] + }, + "scope_status": "changed", + "actual_scope": "Open next question 1 of PR4942 (research/crt-interchange-width-threshold): at a stated ratio rho = C/c_el, the least accepted #4174 reservoir length at which BIT.md's recursive interchange bill C V H^(1-a2) beats the elementary bill c_el V H.", + "scope_change_reason": "The card's generic assumption-proof question is narrowed to an explicitly open, quantified question left by PR4942 about the PR4100 CRT reservoir callee. It is conditional on the stated constant ratio, which BIT.md does not supply.", + "bottleneck_link": "Side result for the PR4100 all-size construction: it says when the a2 width saving of the ordinary interchange can appear at all inside the CRT reservoir, as a function of the unstated constant ratio. No selected-supplier or kappa consequence.", + "budget_minutes": 60, + "stop_rule": "Stop at the exact crossover width for every rho, the critical ratio, and the exact least length where the search is affordable; no claim about the constants themselves.", + "outcome": "completed", + "evidence": "Written argument in README.md; verify.py certifies H*(rho) with two logarithm engines (rigorous integer interval arithmetic and correctly rounded decimal), the critical ratio to 30 digits, and the least layout length by a pruned search cross-checked by a complete brute scan on five small fixtures. RESULT.json is regenerated and compared byte for byte.", + "author": "Nick Petalas", + "ai_assistance": "Independent human and formal proof-kernel review are unclaimed.", + "open_questions": [ + { + "question": "What are the actual constants C and c_el of BIT.md's recursive and elementary interchange on fixed tapes, including the (q^E/2^H)^2 padding factor reported by PR4944?", + "why_it_matters": "The crossover is ~rho^577 once rho exceeds 1.0302; a definite ratio decides whether any physical reservoir ever uses the recursive branch profitably.", + "first_step": "Bound c_el from the elementary q-digit swap and C from the body word's stage inventory at w=2, D=19.", + "blocked_by": "BIT.md states neither constant." + }, + { + "question": "Exact least reservoir length for rho in (1.0302, 1.05] beyond the H* <= 10^9 range computed here.", + "why_it_matters": "Closes the length table between the branch-bound and ratio-bound regimes.", + "first_step": "Replace the k scan by a divisor search near 64*ceil(H*/64) with the d >= 65536 constraint.", + "blocked_by": "Search cost grows linearly in H*." + } + ] +} diff --git a/research/crt-recursive-interchange-crossover/VERIFY.json b/research/crt-recursive-interchange-crossover/VERIFY.json new file mode 100644 index 0000000000..15f6f38ce5 --- /dev/null +++ b/research/crt-recursive-interchange-crossover/VERIFY.json @@ -0,0 +1,21 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [ + "research/crt-reservoir-placement/README.md", + "research/unconditional-batched-crt/BIT.md" + ], + "commands": [ + { + "argv": [ + "{python}", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 120, + "checks": "Least width H with H^a2 > rho (two log engines), critical ratio rho_c = 27706412^a2 to 30 digits, exact least #4174 reservoir length at the branch threshold (2770643200) and at rho = 1.031, 1.032, brute-scan cross-check on five fixtures, five rejection controls, RESULT.json byte comparison." + } + ] +} diff --git a/research/crt-recursive-interchange-crossover/verify.py b/research/crt-recursive-interchange-crossover/verify.py new file mode 100644 index 0000000000..07b97c186a --- /dev/null +++ b/research/crt-recursive-interchange-crossover/verify.py @@ -0,0 +1,306 @@ +"""Crossover width and reservoir length for the recursive ordinary interchange. + +Two engines decide the least width H with H**a2 > rho: a rigorous fixed-point +interval logarithm (integer arithmetic only) and Python's correctly rounded +decimal ln/exp. The layout search then finds the least accepted #4174 +reservoir length L whose width H = floor(L/(100G))*G reaches a given width. +Standard library only. Exit status 0 means every check passed. +""" +from __future__ import annotations + +import decimal +import hashlib +import json +import pathlib +import sys +from fractions import Fraction + +ROOT = pathlib.Path(__file__).resolve().parents[2] +PKG = pathlib.Path(__file__).resolve().parent + +C_BIT = Fraction(8667, 5000000) +A2 = C_BIT * (1 - C_BIT * C_BIT) +H_BRANCH = 27706412 # PR4942: least binary width on BIT.md's recursive branch +D_MIN_EXPONENT = 16 # least power of two d >= 51200 +RHOS = ["1/2", "1", "103/100", "1031/1000", "1032/1000", "21/20", "11/10", + "3/2", "2", "4", "10"] + + +def need(cond: bool, msg: str) -> None: + if not cond: + raise SystemExit("FAIL: " + msg) + + +# ---------- engine 1: rigorous fixed-point interval logarithm ---------- + +def _atanh_bounds(num: int, den: int, prec: int) -> tuple[int, int]: + """Integers lo, hi with lo <= 2**prec * atanh(num/den) <= hi, 0<=num/den<=1/3.""" + need(0 <= 3 * num <= den, "atanh argument outside [0,1/3]") + one = 1 << prec + if num == 0: + return 0, 0 + z_lo = (num * one) // den + z_hi = -((-num * one) // den) + z2_lo = (z_lo * z_lo) >> prec + z2_hi = -((-(z_hi * z_hi)) >> prec) + lo = hi = 0 + p_lo, p_hi = z_lo, z_hi + j = 0 + while p_hi > 2: + lo += p_lo // (2 * j + 1) + hi += -((-p_hi) // (2 * j + 1)) + p_lo = (p_lo * z2_lo) >> prec + p_hi = -((-(p_hi * z2_hi)) >> prec) + j += 1 + # Floors of lower bounds stay lower bounds; ceilings stay upper bounds. + # Tail: sum_{i>=j} z^(2i+1)/(2i+1) <= z^(2j+1)/(1-z^2) <= (9/8) p_hi ulps. + return lo, hi + 4 + + +def ln_bounds(x: Fraction, prec: int) -> tuple[int, int]: + """Integers lo, hi with lo <= 2**prec * ln(x) <= hi, for rational x > 0.""" + need(x > 0, "ln of nonpositive value") + e = x.numerator.bit_length() - x.denominator.bit_length() + y = x / (Fraction(2) ** e) + while y >= 2: + y /= 2 + e += 1 + while y < 1: + y *= 2 + e -= 1 + z = (y - 1) / (y + 1) # in [0, 1/3) + a_lo, a_hi = _atanh_bounds(z.numerator, z.denominator, prec) + l2_lo, l2_hi = _atanh_bounds(1, 3, prec) + lo = 2 * a_lo + (2 * e * l2_lo if e >= 0 else 2 * e * l2_hi) + hi = 2 * a_hi + (2 * e * l2_hi if e >= 0 else 2 * e * l2_lo) + return lo, hi + + +def power_exceeds(h: int, rho: Fraction) -> bool: + """Exact decision of h**a2 > rho via certified intervals (fails if unresolved).""" + need(h >= 1 and rho > 0, "power_exceeds domain") + if h == 1: + return rho < 1 # 1**a2 = 1 exactly + prec = 2 * h.bit_length() + 2 * max(rho.numerator.bit_length(), rho.denominator.bit_length()) + 128 + hl, hh = ln_bounds(Fraction(h), prec) + rl, rh = ln_bounds(rho, prec) + a, b = A2.numerator, A2.denominator + left_lo = (hl * a) // b if hl >= 0 else -((-hl * a + b - 1) // b) + left_hi = -((-hh * a) // b) if hh >= 0 else -((-hh * a) // b) + if left_lo > rh: + return True + if left_hi < rl: + return False + raise SystemExit(f"FAIL: unresolved comparison h={h} rho={rho}") + + +# ---------- engine 2: correctly rounded decimal ln/exp ---------- + +def least_width_decimal(rho: Fraction) -> int: + """Candidate least H with H**a2 > rho from decimal ln/exp.""" + if rho < 1: + return 1 + if rho == 1: + return 2 + ctx = decimal.Context(prec=60) + lr = ctx.ln(ctx.divide(decimal.Decimal(rho.numerator), decimal.Decimal(rho.denominator))) + x = ctx.divide(lr, ctx.divide(decimal.Decimal(A2.numerator), decimal.Decimal(A2.denominator))) + digits = int(x / decimal.Decimal(2.302585)) + 40 + ctx = decimal.Context(prec=digits + 40) + lr = ctx.ln(ctx.divide(decimal.Decimal(rho.numerator), decimal.Decimal(rho.denominator))) + x = ctx.divide(lr, ctx.divide(decimal.Decimal(A2.numerator), decimal.Decimal(A2.denominator))) + return int(ctx.exp(x).to_integral_value(rounding=decimal.ROUND_FLOOR)) + 1 + + +def least_width(rho: Fraction) -> int: + """Least integer H >= 1 with H**a2 > rho, certified by engine 1.""" + h = least_width_decimal(rho) + need(power_exceeds(h, rho), f"engine 1 rejects candidate {h} for rho={rho}") + if h > 1: + need(not power_exceeds(h - 1, rho), f"engine 1 finds smaller width than {h} for rho={rho}") + return h + + +# ---------- layout search over the #4174 accepted domain ---------- + +def least_length(h_star: int) -> tuple[int, dict]: + """Least L = sum w_i over accepted layouts with floor(L/(100G))*G >= h_star. + + Domain (crt-reservoir-placement README): d = 2**s >= 51200, integers + 2 <= m <= ell with m = ceil(log2(p+2)) for some p >= 1, every w_i in + {ell-1, ell}, so L ranges over all integers in [d(ell-1), d*ell]; + G = 64(ell+m) = 64k with ell+2 <= k <= 2*ell. + """ + best = None + witness: dict = {} + k = 4 + while True: + g = 64 * k + r = -(-h_star // g) + t = 100 * g * r # least L with floor(L/(100G)) >= r + if best is not None and 6400 * k > best: + break # t >= 6400k and every L >= t for this and larger k + ell_min, ell_max = (k + 1) // 2, k - 2 + if ell_min <= ell_max: + # d*ell >= t needs d >= t/ell_max; start at the first such power of two + s = max(D_MIN_EXPONENT, (-(-t // ell_max) - 1).bit_length()) + while True: + d = 1 << s + ell_t = max(-(-t // d), 1) + if ell_t > ell_max: + s += 1 + continue + if ell_t >= ell_min: + cand, ell = t, ell_t + else: + ell = ell_min + cand = max(t, d * (ell - 1)) + if best is None or cand < best: + best = cand + witness = {"L": cand, "d": d, "ell": ell, "m": k - ell, "G": g, + "r_max": cand // (100 * g), "H": (cand // (100 * g)) * g} + if ell_t <= ell_min: + break # larger d only raises d*(ell_min-1) + s += 1 + k += 1 + need(best is not None, "no accepted layout") + return best, witness + + +def check_layout_witness(w: dict, h_star: int) -> None: + d, ell, m, L = w["d"], w["ell"], w["m"], w["L"] + need(d >= 51200 and d & (d - 1) == 0, "d not an admissible power of two") + need(2 <= m <= ell, "m outside [2, ell]") + p = (1 << (m - 1)) - 1 if m > 2 else 1 + need(p >= 1 and (p + 1).bit_length() == m, "no p with ceil(log2(p+2)) = m") + need(d * (ell - 1) <= L <= d * ell, "L not a sum of d widths in {ell-1, ell}") + g = 64 * (ell + m) + need((L // (100 * g)) * g >= h_star, "witness width below target") + + +def brute_least_length(h_star: int, max_s: int) -> int: + """Independent slow scan over (s, ell, m) for small fixtures.""" + best = None + for s in range(D_MIN_EXPONENT, max_s + 1): + d = 1 << s + for ell in range(2, 400): + for m in range(2, ell + 1): + g = 64 * (ell + m) + r = -(-h_star // g) + lo = max(d * (ell - 1), 100 * g * r) + if lo <= d * ell and (best is None or lo < best): + best = lo + need(best is not None, "brute scan found nothing") + return best + + +def main() -> None: + bit = (ROOT / "research/unconditional-batched-crt/BIT.md").read_text() + need("c=8667/5000000," in bit and "a2=c(1-c^2)=0.001733394791695384296." in bit, + "BIT.md constants not found") + place = (ROOT / "research/crt-reservoir-placement/README.md").read_text() + for phrase in ["`ell>=max(2,ceil(log2(p+2)))`", "`w_i in {ell-1,ell}`", + "`G=64(ell+ceil(log2(p+2)))`, `r_max=floor(L/(100G))`, `H=r_max G`", + "power-of-two `d>=51200`"]: + need(phrase in place, f"placement domain phrase missing: {phrase}") + need(abs(float(A2) - 0.001733394791695384296) < 1e-18, "a2 value") + + # engine sanity: ln 2 and ln 10 against decimal + for v in (2, 10, Fraction(3, 7)): + lo, hi = ln_bounds(Fraction(v), 200) + c80 = decimal.Context(prec=80) + fv = Fraction(v) + dv = c80.ln(c80.divide(decimal.Decimal(fv.numerator), decimal.Decimal(fv.denominator))) + scaled = c80.multiply(dv, c80.power(decimal.Decimal(2), 200)) + need(lo <= scaled <= hi and hi - lo < 4096, f"ln bounds for {v}") + + # critical ratio rho_c = H_BRANCH**a2 + prec = 400 + hl, hh = ln_bounds(Fraction(H_BRANCH), prec) + ctx = decimal.Context(prec=60) + two_p = ctx.power(decimal.Decimal(2), prec) + a2d = ctx.divide(decimal.Decimal(A2.numerator), decimal.Decimal(A2.denominator)) + rc_lo = ctx.exp(ctx.multiply(a2d, ctx.divide(decimal.Decimal(hl), two_p))) + rc_hi = ctx.exp(ctx.multiply(a2d, ctx.divide(decimal.Decimal(hh), two_p))) + rho_c = str(rc_lo.quantize(decimal.Decimal("1e-30"), rounding=decimal.ROUND_FLOOR, context=ctx)) + need(str(rc_hi.quantize(decimal.Decimal("1e-30"), rounding=decimal.ROUND_FLOOR, context=ctx)) == rho_c, "rho_c digits unstable") + rc_frac_lo = Fraction(rho_c) + rc_frac_hi = rc_frac_lo + Fraction(1, 10**30) + need(power_exceeds(H_BRANCH, rc_frac_lo) and not power_exceeds(H_BRANCH, rc_frac_hi), + "rho_c bracket") + # just below / above rho_c + need(least_width(rc_frac_lo) <= H_BRANCH, "below rho_c the branch threshold binds") + need(least_width(rc_frac_hi) == H_BRANCH + 1, "just above rho_c the width is H_BRANCH+1") + + # exact least length at the branch threshold, with independent brute check + l_branch, w_branch = least_length(H_BRANCH) + check_layout_witness(w_branch, H_BRANCH) + need(l_branch >= 100 * H_BRANCH, "PR4942 lower bound violated") + # small fixtures: fast search vs brute scan + fixtures = [] + for h in (10752, 15360, 27706, 100000, 123457): + fast, wit = least_length(h) + check_layout_witness(wit, h) + # completeness of the brute ranges: any L <= fast has ell-1 <= L/65536 + # (so ell < 400) and d <= L (so s <= 30) + need(fast <= 65536 * 398 and fast <= 1 << 30, f"brute ranges too small at H={h}") + slow = brute_least_length(h, 30) + need(fast == slow, f"layout search mismatch at H={h}: {fast} vs {slow}") + fixtures.append({"H_target": h, "least_L": fast}) + + # rejection controls: each must fail + def rejects(fn) -> bool: + try: + fn() + except SystemExit: + return True + return False + bad = dict(w_branch, L=w_branch["L"] - 1) + need(rejects(lambda: check_layout_witness(bad, H_BRANCH)), "control: L-1 accepted") + need(rejects(lambda: check_layout_witness(dict(w_branch, m=1), H_BRANCH)), "control: m=1 accepted") + need(rejects(lambda: check_layout_witness(dict(w_branch, d=1 << 15), H_BRANCH)), "control: d=2^15 accepted") + need(not power_exceeds(H_BRANCH, rc_frac_hi), "control: width above rho_c accepted") + need(rejects(lambda: power_exceeds(2, Fraction(0))), "control: nonpositive rho accepted") + + rows = [] + for rs in RHOS: + rho = Fraction(rs) + hw = least_width(rho) + h_star = max(H_BRANCH, hw) + row = {"rho": rs, "least_H_with_H^a2>rho": str(hw) if hw < 10**30 else None, + "binding": "branch threshold" if hw <= H_BRANCH else "bill ratio", + "H_star": str(h_star) if h_star < 10**30 else None, + "H_star_bits": h_star.bit_length()} + if h_star == H_BRANCH: + row["least_L"] = l_branch + elif h_star <= 10**9: + l_row, w_row = least_length(h_star) + check_layout_witness(w_row, h_star) + row["least_L"] = l_row + row["witness"] = w_row + else: + row["L_lower_bound_100H"] = str(100 * h_star) if h_star < 10**30 else f"100*H_star, H_star has {h_star.bit_length()} bits" + rows.append(row) + need((rho < rc_frac_lo) == (h_star == H_BRANCH) or rs in ("1031/1000",), f"regime classification {rs}") + + out = { + "a2": f"{A2.numerator}/{A2.denominator}", + "H_branch": H_BRANCH, + "rho_c_floor_1e-30": rho_c, + "least_L_at_branch": l_branch, + "witness_at_branch": w_branch, + "layout_fixtures": fixtures, + "rows": rows, + } + text = json.dumps(out, indent=2, sort_keys=True) + "\n" + target = PKG / "RESULT.json" + if "--write" in sys.argv: + target.write_text(text) + else: + need(target.read_text() == text, "RESULT.json differs from fresh rebuild") + print(f"PASS_CRT_RECURSIVE_INTERCHANGE_CROSSOVER rho_c={rho_c} " + f"least_L_at_branch={l_branch} sha256={hashlib.sha256(text.encode()).hexdigest()[:16]}") + + +if __name__ == "__main__": + main()