diff --git a/.github/workflows/c8-elder-failure-coefficient.yml b/.github/workflows/c8-elder-failure-coefficient.yml new file mode 100644 index 000000000..c2fdf619e --- /dev/null +++ b/.github/workflows/c8-elder-failure-coefficient.yml @@ -0,0 +1,82 @@ +name: C8 elder-failure coefficient (certified replay) + +'on': + pull_request: + paths: + - 'frontiers/c8_elder_failure_coefficient_20261001/**' + - '.github/workflows/c8-elder-failure-coefficient.yml' + - 'frontiers/c6_cluster_coefficients_numerics_20260930/NOTE.md' + - 'frontiers/c6_cluster_coefficients_numerics_20260930/RESULTS.json' + - 'frontiers/local_elder_geometry_20260930/PROOF.md' + - 'imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md' + - 'frontiers/c6_cluster_law_20260929/PROOF.md' + - 'frontiers/planar_cubic_cluster_20260929/PROOF.md' + workflow_dispatch: + +permissions: + contents: read + +jobs: + certified-replay: + runs-on: ubuntu-latest + timeout-minutes: 150 + strategy: + fail-fast: false + matrix: + flags: ['-B -S', '-B -O -S'] + steps: + - uses: actions/checkout@11d5960a326750d5838078e36cf38b85af677262 + with: + persist-credentials: false + - uses: actions/setup-python@a26af69be951a213d495a4c3e4e4022e16d87065 + with: + python-version: '3.11.16' + - name: Verify manifest and main-resident pins + env: + PYTHONDONTWRITEBYTECODE: '1' + run: | + set -euo pipefail + python -B -S - <<'PY' + import hashlib, json, pathlib + root = pathlib.Path('frontiers/c8_elder_failure_coefficient_20261001').resolve() + manifest = json.loads((root / 'SOURCE_FILES.json').read_text())['files'] + present = sorted(p.name for p in root.iterdir()) + if present != sorted([e['path'] for e in manifest] + ['SOURCE_FILES.json']): + raise ValueError('packet tree differs from manifest: ' + repr(present)) + for entry in manifest: + path = root / entry['path'] + if pathlib.PurePosixPath(entry['path']).name != entry['path'] or path.is_symlink(): + raise ValueError('flat regular source required: ' + entry['path']) + data = path.read_bytes() + blob = hashlib.sha1(b'blob %d\0' % len(data) + data).hexdigest() + if len(data) != entry['bytes'] or hashlib.sha256(data).hexdigest() != entry['sha256'] or blob != entry['git_blob']: + raise ValueError('source identity mismatch: ' + entry['path']) + pins = json.loads((root / 'SOURCE_MAP.json').read_text()) + for entry in pins['sources']: + if not entry['current_required']: + continue + data = pathlib.Path(entry['path']).read_bytes() + blob = hashlib.sha1(b'blob %d\0' % len(data) + data).hexdigest() + if len(data) != entry['bytes'] or hashlib.sha256(data).hexdigest() != entry['sha256'] or blob != entry['git_blob']: + raise ValueError('consumed source identity mismatch on this tree: ' + entry['path']) + print('manifest and pins verified') + PY + - name: Self-test, full replay against RESULTS.json, mutants + env: + PYTHONDONTWRITEBYTECODE: '1' + FLAGS: ${{ matrix.flags }} + run: | + set -euo pipefail + cd frontiers/c8_elder_failure_coefficient_20261001 + python $FLAGS certificate.py --check --procs "$(nproc)" + for m in vw-sign ww-drop corner-swap mid-branch k-weight; do + set +e + python $FLAGS certificate.py --check --mutant "$m" --procs "$(nproc)" > /dev/null + rc=$? + set -e + if [ "$rc" -ne 1 ]; then echo "mutant $m not rejected (rc=$rc)"; exit 1; fi + echo "mutant $m rejected" + done + cd "$GITHUB_WORKSPACE" + git diff --exit-code + test -z "$(git status --porcelain)" diff --git a/frontiers/c8_elder_failure_coefficient_20261001/NOTE.md b/frontiers/c8_elder_failure_coefficient_20261001/NOTE.md new file mode 100644 index 000000000..9903dd6f6 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/NOTE.md @@ -0,0 +1,411 @@ +# Certified planar near cluster coefficients `alpha_1`, `alpha_2` and the elder-failure constant `C_fail^{B,K}` + +**Object:** CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1. +**Author:** Anthropic Claude (Claude Code session `session_015wNj8LPTKXsaT68G3DgPPh`), 1 October 2026. +**Disposition:** certified numerical enclosures. Binary64 interval arithmetic with outward rounding encloses explicit +finite-dimensional Gaussian integrals that the sources define. Author-side. Nonauthor reads at `b54ea71` (OpenAI / +Codex): Slices A and C ACCEPT (review 5378890193), Slice B ACCEPT (review 5379043497); v1.3 applies their three +non-blocking clarifications (the `chi^2` coordinates, outward comparison brackets, the `Phi`-series bound). +**Scientific effect:** NONE. No `STATUS`, `PROOF_INDEX`, `GRAPH`, claim, catalog, prize or source-body change. Same GitHub +account as every lane; zero organizational-independence credit. + +## 0. Statement + +Planar model, reference (continuum) kernel `K(z) = exp(-|z|^2/2)`. `(a, beta, c) ~ N(0, diag(2, 2, 6))` are the odd +contact jets and `I_1`, `I_2` the near failure integrals of [NUM] §§2–3 (recalled in §1). Write `J_j(k) = E[I_j]` and +`J_fail = J_1 + J_2`. By [NUM] §2, from [CL] Theorem N or [CUB] Theorem F, + + alpha_j(b, k) = (p_b(0)/z_0) J_j(k), p_b(0)/z_0 = phi(b/sqrt2) / (sqrt2 36 k^2 m_(2,b)), + m_(2,b) = (b^2 + 2) Phi(b/sqrt2) + sqrt2 b phi(b/sqrt2), + +with `nu(2) = alpha_2` ([NUM] §1) and `a_fail = alpha_1 + alpha_2` the selector-failure coefficient of [S] (S2). + +**Theorem (certified enclosures).** The following intervals contain the exact values for the reference kernel. By Lemma T +(§6a) they also contain them for the torus model of [LP] and [S], with every side `L >= 10` and every frame. The torus +enclosures in `RESULTS.json` agree with these at the displayed precision. The certificate's own intervals (`RESULTS.json`, +10 digits) are rounded further outward here. + +| `k` | `J_fail = J_1 + J_2` | `J_1` | `J_2` | `J_2 / J_fail` | near Palm excess `2J_2/(J_1 + 2J_2)` | +|---|---|---|---|---|---| +| `1/2` | `[36.783058, 36.786801]` | `[35.673609, 35.680652]` | `[1.10614, 1.10945]` | `[0.030069, 0.030162]` | `[0.058377, 0.058563]` | +| `1` | `[172.84636, 172.86021]` | `[169.58395, 169.60709]` | `[3.25312, 3.26241]` | `[0.018819, 0.018875]` | `[0.036941, 0.037052]` | +| `2` | `[1111.9822, 1112.0308]` | `[1103.1296, 1103.2057]` | `[8.82508, 8.85263]` | `[0.0079360, 0.0079612]` | `[0.015746, 0.015797]` | + +| `k` | `b` | `p_b(0)/z_0` | `alpha_1` | `alpha_2 = nu(2)` | `a_fail = alpha_1 + alpha_2` | +|---|---|---|---|---|---| +| `1/2` | `0` | `[0.031343862, 0.031343869]` | `[1.118148, 1.118370]` | `[0.034670, 0.034775]` | `[1.152923, 1.153041]` | +| `1/2` | `1` | `[0.0089740288, 0.0089740307]` | `[0.3201360, 0.3201993]` | `[0.0099266, 0.0099563]` | `[0.3300922, 0.3301259]` | +| `1` | `0` | `[0.0078359656, 0.0078359673]` | `[1.328854, 1.329036]` | `[0.025491, 0.025565]` | `[1.354418, 1.354527]` | +| `1` | `1` | `[0.0022435072, 0.0022435077]` | `[0.3804628, 0.3805148]` | `[0.0072984, 0.0073193]` | `[0.3877820, 0.3878132]` | +| `2` | `0` | `[0.0019589914, 0.0019589919]` | `[2.161021, 2.161171]` | `[0.017288, 0.017343]` | `[2.178363, 2.178459]` | +| `2` | `1` | `[0.00056087680, 0.00056087692]` | `[0.6187198, 0.6187626]` | `[0.0049497, 0.0049653]` | `[0.6236851, 0.6237124]` | + +**Corollary 1 (the elder-failure constant).** On `B x K = [0, 1] x [1/2, 2]`, the coefficient of [S] §9 for the reference +kernel, + + C_fail^{B,K} = integral_(B x K x S^1) A_0^ref a_fail / (3 k^(5/3)) db dk dsigma(u) = const I_b I_k, + const = 4 / (sqrt2 (2 pi)^(5/2) sqrt12), I_b = integral_0^1 exp(-b^2) db, + I_k = integral_(1/2)^2 exp(-12 k^2) k^(-5/3) J_fail(k) dk, + +satisfies + + I_k in [0.4421366, 0.4432234], C_fail^{B,K} in [0.00272445, 0.00273116]. + +By Lemma T the same interval holds for the torus model with every `L >= 10`. That is the model of [S] §9, which states +`nu_rej^{B,K}(ell) ~ C_fail^{B,K} ell^(2/3)` and `E N_rej^{B,K}(0, t] ~ (3/5) C_fail^{B,K} t^(5/3)`. + +**Corollary 2 (`d = 3`, conditional on Math-#175 and Math-#184).** +- Math-#175 (OpenAI/Codex, open, AUTHOR_SIDE / HOLD) is [FIB]: + - its Theorem F, (S3), identifies `(1 - p_r)/r^3 -> a_1 + a_2` in fixed `d >= 3`; + - its (L1) is `C_fail^{B,K} = integral_(B x K x S^(d-1)) A_0 a_fail / (3k^(5/3))`. +- Math-#184 (open) is [HD]. Its Theorem 1 gives `a_j^(d) = R_d(b) alpha_j` with `R_d = N_d m_(2,b) / m_(d,b)`. + +In `d = 3`, for the reference kernel, the derivation of §1 gives +`A_0^ref,(3) = 432 (2 pi)^-4 12^(-1/2) k^2 exp(-12k^2) exp(-3b^2/4) m_(3,b)`: the transverse jet `f_xy2` adds a factor +`(2 pi)^(-1/2)`, `p_G(0) = (2 pi)^(-3/2)` and `z_0 = 36 k^2 m_(3,b)`. Both `m_(3,b)` and `m_(2,b)` cancel. With +`N_3(b) = pi (2 pi)^(-1/2) g(b)` ([HD] (2), `m = 2`) we have `exp(-3b^2/4) phi(b/sqrt2) N_3(b) = exp(-b^2) g(b) / 2`, so + + C_fail^(3),{B,K} = 4 pi (4/sqrt2) (2 pi)^-4 12^(-1/2) I_b^(3) I_k, I_b^(3) = (1/2) integral_0^1 exp(-b^2) g(b) db, + +with `g(b) = (b^3 + 6b) Phi(b/sqrt2) + sqrt2 (b^2 + 4) phi(b/sqrt2)`. The certificate checks `N_3(0) = 2 sqrt2` ([HD] (3)) and +gives + + I_b^(3) in [1.513649, 1.513684], C_fail^(3),{B,K} in [0.00440581, 0.00441675]. + +At `b = 0`, `R_3(0) = (32 + 28 sqrt2)/17`: + +| `k` | `alpha_1^(3)` | `alpha_2^(3) = nu^(3)(2)` | `a_fail^(3)` | +|---|---|---|---| +| `1/2` | `[4.709246, 4.710177]` | `[0.14602, 0.14646]` | `[4.855704, 4.856199]` | +| `1` | `[5.596663, 5.597427]` | `[0.10736, 0.10767]` | `[5.704330, 5.704788]` | +| `2` | `[9.101458, 9.102086]` | `[0.072812, 0.073040]` | `[9.174497, 9.174898]` | + +**Consequences.** +- **`nu(2)` and `a_fail` are certified.** At `(b, k) = (0, 1)`: `alpha_2 = nu(2)` is in `[0.025491, 0.025565]` and `a_fail` in + `[1.354418, 1.354527]`. The two-point share `J_2/J_fail` decreases from + `3.01%` at `k = 1/2` to `0.79%` at `k = 2`. +- **These are the model's constants.** By Lemma T the planar enclosures hold for the torus model of [LP] and [S] with + every `L >= 10` and every frame, not only for the reference kernel. +- **[NUM]'s open digit is settled.** At `k = 2`, [NUM] observed that Gauss–Hermite and Monte Carlo disagree on `J_2` by + 2.6% and said the third digit of `alpha_2(2)` was not established. The certified `J_2(2)` is in `[8.82508, 8.85263]`. It + excludes the GH60 value `8.641` and is consistent with the Monte Carlo value `8.875 +- 0.073`. The Gauss–Hermite + sequence had not converged there. +- **The cap route overstates the lifetime difference by about `1.5 x 10^5`.** Math-#215 (open) certifies + `nu_cand - nu_eld <= C_{B,K}(r_*) ell^(2/3)` through the cap criterion, with `C_{B,K}(1/4096) = 402.5297`. As + `r_* -> 0`, that route can give no less than `Gamma_2 = 307.51` (floating). The leading constant of the difference is + `C_fail^{B,K}`, about `0.00273`: `1.47 x 10^5` times smaller than `C_{B,K}(1/4096)` and `1.13 x 10^5` times + smaller than `Gamma_2`. +- **The same holds in `d = 3`, conditionally.** `C_fail^(3),{B,K}` is about `0.00441`. That is `2.73 x 10^5` times + smaller than Math-#215's `C_{B,K}^(3)(1/4096) = 1205.7089` and `1.23 x 10^5` times smaller than + `Gamma_3 = 540.97`. + +## 1. What is evaluated + +[NUM] §3 (from [CUB] (C2)–(C3) and [CL] (3.12)) defines, with `B = beta - a^2/(12k)`, `D = (c - a beta/(4k) + +a^3/(72k^2))/2`, the typed domain `s < -|B|/2`, the weight `w_0 = 9k^2 (4s^2 - B^2)` and + + g(s) = -2s^3 + 3Bs^2 - B^3 - 12k D^2, + +the classification: `n = 2` iff `s > B` and `g(s) > 0`; `n = 1` iff `g(s) < 0`; `n = 0` iff `s <= B` and `g(s) > 0`. +Then `I_j = integral_(s < -|B|/2) w_0 1{n = j} ds`. [NUM] §2 gives the prefactor `p_b(0)/z_0`. Here `p_b(0)` is the density +at `0` of `A = f_zz | (f, f_xx, f_xz) = (b, 0, 0) ~ N(-b, 2)`, and `z_0 = 36 k^2 m_(2,b)`. + +The reference integrand of [S] §9 is `A_0^ref = 12 pi_0 z_0`. [LP] §15 factorizes `pi_0 = p_(f,V_u)(b, 0) p_G(0) +phi_tau(12k)`. With the jet covariances of [NUM] §2 these factors are: +- `(f, f_xx)` has covariance `[[1, -1], [-1, 3]]`, and `f_xz` is independent with variance `1`. So + `p_(f,V_u)(b, 0) = exp(-3b^2/4) / (2 pi sqrt2 sqrt(2 pi))`. +- `p_G(0) = 1/(2 pi)`. +- `tau^2 = Var(f_xxx | grad f = 0) = 15 - 9 = 6`. + +Hence `pi_0 = (2 pi)^-3 12^(-1/2) exp(-3b^2/4 - 12k^2)` and + + A_0^ref = 432 (2 pi)^-3 12^(-1/2) k^2 exp(-12 k^2) exp(-3b^2/4) m_(2,b). + +In `A_0^ref a_fail / (3 k^(5/3))`, `m_(2,b)` and `k^2` cancel, `exp(-3b^2/4) phi(b/sqrt2) / sqrt2 = exp(-b^2) / (2 sqrt pi)`, +and the circle gives `2 pi`. This leaves `const I_b I_k` with the `const` above. + +## 2. Reduction to one root (Lemma 1) + +Put `x = -s` and `P(x) = (x + B)^2 (2x - B)`. Expanding, `g(s) = P(x) - T` with `T = 12k D^2`. The typed domain is +`x > |B|/2`. + +- **`B >= 0`.** `P` increases on `x >= B/2` from `P(B/2) = 0`. So `n = 1` exactly on `(B/2, x_b)`, and `n = 2` never + occurs (it needs `x < -B`). +- **`B < 0`.** `P` decreases on `(|B|/2, |B|)` from `|B|^3/2` to `0`, then increases. + - If `T < |B|^3/2`: `n = 2` on `(|B|/2, x_m)` and `n = 1` on `(x_m, x_b)`. + - Otherwise: `n = 1` on `(|B|/2, x_b)`. + +Here `x_b >= max(-B, B/2)` is the **big root** and `x_m in [|B|/2, |B|]` the **middle root** of `P = T`. In every case +the failure set `{n >= 1}` is the single interval `|B|/2 < x < x_b` (up to null sets). With the antiderivative +`9k^2 (4x^3/3 - B^2 x)`: + + I_fail = I_1 + I_2 = 3k^2 H, H = |B|^3 + 4 x_b^3 - 3B^2 x_b = 2(B_-)^3 + 2T + V, V = -3B(2x_b - B)(x_b + B), + I_2 = 3k^2 H_2, H_2 = |B|^3 + 4 x_m^3 - 3B^2 x_m on {B < 0, T < |B|^3/2}, 0 elsewhere. + +The identity `2T + V = (x + B)(2x - B)^2 = 4x^3 - 3B^2 x + B^3` holds at either root. The controls `(0, -2k, 0)` (`H_2 = 16k^3`, +`I_2 = 48k^5`) and `(0, 0, c)` (`H = 2T = 24k D^2`, `I_1 = 72k^3 D^2`) agree with [NUM] §3. + +**Scaled coordinates.** `a^ = a / sqrt k ~ N(0, 2/k)`, `c^ = sqrt(k) c ~ N(0, 6k)`. Then `B = beta - a^2/12` and +`T = 3u^2` with `u = c^ + d` and `d = -a^ beta/4 + a^3/72`, so `J_fail(k) = 3k^2 E[H]` and `J_2(k) = 3k^2 E[H_2]` with +`k`-free integrands. Exactly, `E[2T] = 6 E[u^2] = 6 (6k + 1/(4k) + 5/(216 k^3))`. Also `E[2(B_-)^3] = 4 integral_0^inf +phi_(sqrt(2/k))(a) g(a^2/12) da` with `g(m) = E[(m - beta)_+^3]`. + +## 3. Monotonicity and the `(B, w)`-calculus (Lemma 2) + +Let `r = B/x`, so `r in [-1, 2]` on the big root and `r in [-2, -1]` on the middle root. Let `w = |u| = (T/3)^(1/2)`. +Implicit differentiation of `P(x) = T` gives: +- `x_T = 1/(6x(x + B))` (positive on the big root, negative on the middle root); +- `x_B = -(x - B)/(2x)`; +- `r_B = (2x - B)(x + B)/(2x^3)`; +- `H_T = (4 - r^2)/(2(1 + r))`; +- `V_T = -r(r + 4)/(2(1 + r))`; +- `H_B / x^2 = -(3/2)(r - 2)^2 (r + 1)` for `B > 0` and `-(3/2)(r + 2)(r^2 - r + 2)` for `B < 0`; +- `V_B = -(3/2) x^2 (r - 2)^2 (r + 1)`. + +Therefore: +- `H` is non-increasing in `B` and non-decreasing in `T`. +- `V` is non-increasing in `B`. In `T` it is non-increasing for `B > 0` and non-decreasing for `B < 0`. +- `H_2` is non-increasing in both. +- On the big root, `x` increases in `T` and, for fixed `T`, is unimodal in `B` with minimum `x = B` at `B = (T/4)^(1/3)`. + `r` is non-decreasing in `B`, increasing in `T` for `B < 0` and decreasing for `B > 0`. +- On the middle root, `x` and `r` are non-increasing in both `B` and `T`. + +In `(B, w)` with `sigma = +1` on the big root and `-1` on the middle root: + + V_w = -sigma sqrt3 x^(3/2) r (4 + r)(2 - r)^(1/2) + V_BB = x p(r), p(r) = 6 + 3r - 9r^2/2 - 3r^3/4 + 3r^4/4 + V_Bw = -sigma (sqrt3/2) x^(1/2) e(r), e(r) = (2 - r)^(3/2)(4 + 2r + r^2) (decreasing on r <= 2) + V_ww = g(r) = -r(4 - r + r^2) (strictly decreasing) + +- The `T`-derivatives blow up as `T -> 0` for `B < 0` (`V_TT ~ T^(-3/2)`). The `(B, w)`-Hessian is bounded: its entries are + `x p(r)`, `x^(1/2) e(r)` and `g(r)` with `r` in a bounded range. At the triple point `(0, 0)` it is bounded but + discontinuous, and `V` is `C^(1,1)` there, which is all the integral remainder needs. +- `V` depends on `u` only through `w`, with a kink at `u = 0` when `B < 0` (`V_w(B, 0) = 9|B|^(3/2)`). +- On the support boundary `T = |B|^3/2` the middle root is `|B|/2` (`r = -2`). There `H_2 = 0`, + `H_2,B = -(3/2) x^2 (r + 2)(r^2 - r + 2) = 0` and `H_2,T = (4 - r^2)/(2(1 + r)) = 0`. So the zero extension of `H_2` is + `C^1` and piecewise `C^2`. + +All formulas were checked against finite differences of the root map (relative agreement `6e-6` at step `1e-4`) and +symbolically. + +## 4. Box enclosures (Lemma 3) + +For a box `Q` of jet coordinates, let `[B_l, B_h] x [T_l, T_h]` contain its image. The ranges of `B` and `d` are exact +up to outward rounding: `d` is linear in `beta`, and its `a`-extremes are at the ends or at `a = +-sqrt(6 beta)`. The +`u`-range is the `c`-range plus the `d`-range. Two enclosures of `integral_Q F dmu` (`F = V`, `H`, `H_2`) are formed and intersected; an empty +intersection would raise. + +- **First order.** `I_1 = mu(Q) [min F, max F]`. The extremes sit at rectangle corners by Lemma 2, so the corner roots + suffice (four big roots, or two middle roots). +- **Second order** (only if the `u`-range of `Q` avoids `0`, so `w = s u` with a fixed sign `s`). Take a float point + `(B*, w*)` of the rectangle and Taylor's formula with integral remainder along the segment to `(B, w)` (the rectangle + is convex): + + F = F* + F_B* dB + F_w* dw + 1/2 (h_BB dB^2 + 2 h_Bw dB dw + h_ww dw^2). + + - Each averaged Hessian entry `h` lies in the interval range `S` of the second derivative over the rectangle, + computed from the corner ranges of `x` and `r`, with `p` evaluated in mean-value form. For `H_2` on boxes meeting + the support boundary, the ranges are hulled with `0`. + - Hence `integral_Q F dmu` lies in + + F* m + F_B* E[dB] + s F_w* E[du] + 1/2 (S_BB E[dB^2] + S_ww E[du^2]) + s c E[dB du] +- delta (E[dB^2] E[du^2])^(1/2), + + where `c` and `delta` are the centre and radius of `S_Bw`. + - The six moments are exact Gaussian polynomial moments (orders `a <= 6`, `beta, c <= 2`). The second-order part + has width `O(diam(S) h^2 mu(Q))`, so the scheme is third order. +- **Partial moments.** `M_n = (n - 1) s^2 M_(n-2) + s^2 (x_0^(n-1) phi(x_0) - x_1^(n-1) phi(x_1))`. Cell masses are + differences of upper tails `Q(t)` on the side away from `0`, so they keep relative precision far out. +- **Roots.** The roots are bracketed by Newton iteration from a Cardano / trigonometric start. The bracket is then + verified: `P(a) - T < 0 < P(b) - T` for the big root, reversed for the middle root. The sign test is a binary64 + evaluation of `(x + B)^2 (2x - B) - T` with an a priori rounding bound. +- **Refinement.** An initial `8^3` grid covers `[-8 sigma, 8 sigma]^3`. Boxes are bisected depth-first along their + standardized widest side until the enclosure width is at most `10^-8`. Partial sums are rounded outward, and the + initial cells are dealt to processes in a fixed pattern. + +## 5. Exact and one-dimensional parts; tails + +- `E[2T]` is the exact rational above. +- `E[2(B_-)^3]` is a one-dimensional integral over `4000` cells, each by second-order Taylor expansion about the + midpoint. + - `G(a) = g(a^2/12)` has `G'' = g''(m) a^2/36 + g'(m)/6`. + - `g' = 3[(m^2 + 2) Phi(m/sqrt2) + sqrt2 m phi(m/sqrt2)]` and `g'' = 6[m Phi(m/sqrt2) + sqrt2 phi(m/sqrt2)]` are increasing, + so `G''` is bounded by its values at the cell ends. + - The tail `a > 12 sigma` uses `g(m) <= m^3 + 6m + sqrt2 (m^2 + 4)/sqrt(2 pi)`. +- **Outside the box.** Lemma 1 and `x_b <= |B| + (T/2)^(1/3)` (since `P(|B| + t) >= 2t^3`) give + `|V|, H, H_2 <= 17|B|^3 + 8T`. This is a polynomial in `|a|, |beta|, |c|`, integrated over the complement by the union + bound with exact tail moments. + +These are expectations over the scaled law. `J_fail = 3k^2 (E[2(B_-)^3] + E[2T] + E[V])` and `J_2 = 3k^2 E[H_2]`. The +tail bound enters `E[V]` with both signs and `E[H_2]` from above. + +| `k` | `E[2(B_-)^3]` | `E[2T]` exact | `E[V]` boxes | `E[H_2]` boxes | tail bound | +|---|---|---|---|---|---| +| `1/2` | `[8.3574312330, 8.3574329380]` | `199/9` (`22.11111111`) | `[18.5755374, 18.5805207]` | `[1.4748651, 1.4792648]` | `1.5E-9` | +| `1` | `[5.8760207620, 5.8760219460]` | `1355/36` (`37.63888888`) | `[14.1005477, 14.1051535]` | `[1.0843739, 1.0874697]` | `3.3E-10` | +| `2` | `[5.0934267220, 5.0934277450]` | `20957/288` (`72.76736111`) | `[14.8044089, 14.8084356]` | `[0.73542405, 0.73771834]` | `2.1E-10` | + +## 6. The `k`-integral + +`I_k` is a four-dimensional box integral in the original coordinates `(a, beta, c, k)`, with weight +`W(k) = 3k^(1/3) exp(-12k^2)` (that is, `3k^2 exp(-12k^2) k^(-5/3)`): + + B = beta - a^2/(12k), u = sqrt(k) c - a beta/(4 sqrt k) + a^3/(72 k^(3/2)). + +- Each 4D box maps into a box of scaled coordinates, so the rectangle machinery of §4 applies unchanged. +- The Taylor moments are polynomials in `(a, beta, c)` times `k^(p/2)`, `p = -6..2`. Their `k`-factors + `K_p = integral W(k) k^(p/2) dk` are enclosed on the `4096` elementary dyadic cells of `[1/2, 2]` by midpoint (below) + and trapezoid (above) rules on `24` sub-cells. This is valid because every `f(k) = k^q exp(-12k^2)` with + `q in [-8/3, 4/3]` is convex on `[1/2, 2]`: `f''/f = (24k - q/k)^2 - q/k^2 - 24 >= 87 - 16/3 - 24` for `q >= 0` and + `>= 120` for `q < 0`. +- `k`-cells are dyadic and never finer than elementary. +- The tail is bounded uniformly in `k` by `17|B|^3 + 8T <= 68|beta|^3 + (68/216 + 576/5184) a^6 + 144 c^2 + 9 a^2 beta^2`. + +Values: +- `int_K W(k) dk` is in `[0.00907557755, 0.00907557942]`; +- `const` is in `[0.008250963314, 0.008250964965]`; +- `I_b = sqrt(pi)(Phi(sqrt2) - 1/2)` is in `[0.7468240581, 0.7468242075]`; +- `I_k` (boxes plus tail) is in `[0.4421366, 0.4432234]`. + +`I_b^(3)` is integrated on `100000` cells of `[0, 1]`. On each cell `exp(-b^2)` decreases and `g` increases (both positive), +so the cell integral lies between `h exp(-b_1^2) g(b_0)` and `h exp(-b_0^2) g(b_1)`. + +## 6a. Lemma T: transfer to the torus + +[LP] and [S] work on the torus `R^2/(L Z^2)` with +`K_L(z) = sum_n exp(-|z + Ln|^2/2) / sum_n exp(-|Ln|^2/2)` ([LP] §2, by Poisson summation). Everything certified in the +plane is a function of **one-site** jet covariances of order at most `6`: +- the odd conditional law of `(a, beta, c)` given `(f_x, f_z, f_xxx) = (0, 0, 12k)`; +- the even conditional law of `A = f_zz` given `(f, f_xx, f_xz) = (b, 0, 0)`, which gives `p_b(0)` and `z_0`; +- the factors of `pi_0` ([LP] §15). + +`K_L` is even, so odd and even jets stay exactly uncorrelated. + +**Lemma T.** For every side `L >= 10` and every frame: +1. **Covariances.** Every one-site covariance of order `<= 6` differs from its reference value by at most + `delta <= 2.86E-14`. In a frame `R`, + `d^g K_L(0) = sum_n (d^g G)(R^T L n)/Theta_L`, `|R^T L n| = L|n|`, and `|He_g1(y1) He_g2(y2)| <= Hh_g1(|y|) Hh_g2(|y|)` with + `Hh_j(t) = sum_k j!/(k!(j-2k)! 2^k) t^(j-2k)`. The shell `|n|_inf = j` has `8j` points with `|n| >= j`. Each product + `Hh_a Hh_c e^(-t^2/2)` decreases for `t^2 > a + c`, so the bound at `L = 10` holds for every `L >= 10`. +2. **Regressions.** Interval regressions with every entry widened by `delta` give the torus conditional laws: + - odd: `N(12k v, S_c)` with `||diag(2,2,6)^(-1/2)(S_c - diag(2,2,6))diag(2,2,6)^(-1/2)|| <= eps = 1.30e-12`; + - even: `N(c_1 b, s^2)`; + - the ratio `(pi_0 p_b(0))^L/(pi_0 p_b(0))^ref` lies in `1 +- 8.8e-12`, uniformly in `(b, k)` in the band; + - `z_0` cancels in `A_0 a_fail` on the torus as well. +3. **Expectations.** Let `S_c = S_0^(1/2)(I + E)S_0^(1/2)` with eigenvalues `e_i` of `E`, and `mu~ = S_0^(-1/2) mu`. + Coordinate-free, `1 + chi^2 = det(I - E^2)^(-1/2) exp(mu~^T (I - E)^(-1) mu~)`. In an orthonormal eigenbasis of `E`, + with `mu~_i` the coordinates of `mu~` in that basis (not the original jet coordinates), this reads + + 1 + chi^2 = prod (1 - e_i^2)^(-1/2) exp(sum mu~_i^2/(1 - e_i)). + + For `eps^2 <= 0.006` we have `(1 - eps^2)^(-3/2) <= 1.01` and `(1 - eps^2)^(-3/2) - 1 <= 2 eps^2`. With + `e^m - 1 <= m e^m`, this gives `chi^2 <= 1.01 m e^m + 2 eps^2`, where `m = |mu~|^2/(1 - eps)`. The certificate checks + `eps^2 <= 0.006` and obtains `chi^2 <= 2.88E-23`. By Cauchy–Schwarz, `|J^L - J| <= (E[I^2] chi^2)^(1/2)`, using + `I = 3k^2 H <= 12 Pi_4` (the tail polynomial of §6, `k <= 2`), with `E[I^2] <= 5.405e+08`. This bounds `|J^L - J|` by + `1.25E-7` for every `k` in `[1/2, 2]`, for `J_fail` + and for `J_2` (`0 <= I_2 <= I_fail`). +4. **Constants.** Hence `C_fail^L` lies in `(1 +- 8.8e-12)(C_fail^ref +- const I_b K_(-4)/3 . |J^L - J|)`, with + `K_(-4)/3 = integral e^(-12k^2) k^(-5/3) dk`. The torus `alpha_j`, `a_fail` at the grid points use the torus + `p_b(0)/z_0`. + +`torus.py` computes all of this in interval arithmetic. `RESULTS.json` lists the torus enclosures under "torus model". +Lemma T is planar: Corollary 2 stays a reference-kernel statement. + +## 7. Arithmetic + +Every `+, -, *, /, sqrt` is IEEE binary64 round-to-nearest, widened outward by one ulp. No libm transcendental is +trusted: +- `pi` comes from Machin's series with a rational remainder, and `ln 2` from `sum 1/(j 2^j)`. +- `exp(-y)`: `y = n ln2 + r` with `r in [0, 0.75]`, a degree-18 Taylor polynomial with the a priori bound + `(56u + 10^-19) e^0.75`, then exact scaling by `2^-n`. The observed relative width is below `4e-13` for `y <= 700`. +- `Phi(t)` for `t < 2` is `1/2 + phi(t) sum t^(2n+1)/(2n+1)!!`, with an a priori summation bound and a geometric tail. +- `Q(t) = 1 - Phi(t)` for `t >= 2` is `phi(t) R(t)`, where `R` is Laplace's continued fraction + `1/(t + 1/(t + 2/(t + ...)))`. It has positive elements, so consecutive approximants bracket `R`. + +libm `pow`, `acos` and `cos` only seed iterations whose results are verified. + +## 8. Verification + +- **`--check`** (CI, under `-B -S` and `-B -O -S`) runs two stages: + - a self-test: every box enclosure, and its first- and second-order parts separately, on `13` fixed 3D boxes and + `2` 4D boxes contains a pure-Python Gauss–Legendre value (`12`–`16` points per axis; relative slack `1e-9`, or + `1e-6` on the `C^1` support boundary of `H_2`), and the `k`-moments contain Gauss–Legendre values; + - the full replay. +- The replayed enclosures must lie inside the published ones and be at most `2 x 10^-7` relatively (plus one unit of + the 10th digit) tighter. The leaf counts must agree: + - `V`: `1,209,983`, `1,096,831`, `951,074` at `k = 1/2, 1, 2`; + - `H_2`: `1,131,985`, `799,321`, `583,149`; + - `I_k` (4D): `997,266`. +- A full replay takes about 10 minutes on four cores. +- **Mutants** exit 1 at the self-test: + - `vw-sign` flips `V_w`; + - `ww-drop` drops `V_ww`; + - `corner-swap` swaps the monotonicity corners of `V`; + - `mid-branch` gives `H_2` the big-root sign of `V_w`; + - `k-weight` drops `k^(1/3)` from `W`. +- **Development checks** (not in the packet): + - `1,200` random 3D boxes and `250` random 4D boxes, `0` failures against `numpy` Gauss–Legendre; + - the full-`H` integrand, boxed directly, gives an enclosure of `J_fail(1/2)` that intersects the split one + (`controls.py`). + +## 9. Comparisons (floating controls, `controls.py`) + +Floating values of [NUM] (`RESULTS.json`, GH40 / GH60 Gauss–Hermite, Monte Carlo with `10^6` samples) against the +certified enclosures. "in" marks a value inside the enclosure; the last column is the Monte Carlo deviation from +the enclosure midpoint in Monte Carlo standard errors. + +| `k` | | certified | GH40 | GH60 | MC `+-` s.e. | (MC - mid)/s.e. | +|---|---|---|---|---|---|---| +| `1/2` | `J_1` | `[35.67361, 35.68065]` | `35.652` (-0.069%) | `35.676` in (-0.004%) | `35.5595 +- 0.0856` | `-1.37` | +| `1/2` | `J_2` | `[1.106149, 1.109449]` | `1.1211` (+1.205%) | `1.1085` in (+0.065%) | `1.11467 +- 0.00619` | `+1.11` | +| `1` | `J_1` | `[169.584, 169.6071]` | `169.63` (+0.023%) | `169.61` (+0.009%) | `169.374 +- 0.283` | `-0.78` | +| `1` | `J_2` | `[3.253122, 3.262409]` | `3.2494` (-0.256%) | `3.2554` in (-0.073%) | `3.28057 +- 0.0222` | `+1.03` | +| `2` | `J_1` | `[1103.13, 1103.206]` | `1104.2` (+0.097%) | `1103.7` (+0.046%) | `1102.63 +- 1.69` | `-0.32` | +| `2` | `J_2` | `[8.825089, 8.85262]` | `8.4211` (-4.726%) | `8.6412` (-2.236%) | `8.87486 +- 0.0727` | `+0.50` | + +Every Monte Carlo value is within `1.4` standard errors of the certified value. GH60 is inside the enclosure for +`J_1(1/2)`, `J_2(1/2)` and `J_2(1)`, within `0.05%` of it for `J_1(1)` and `J_1(2)`, and `2.2%` low for `J_2(2)`. + +Other controls (`controls.py`): +- the floating Simpson value of `E[2(B_-)^3]` lies inside its enclosure at all three `k`; +- the direct box enclosure of `3k^2 E[H]` at `k = 1/2`, `[36.778523711, 36.792370494]` (`444127` leaves, tolerance + `1e-7`), intersects the split enclosure of `J_fail(1/2)`; +- with the displayed numerators, `C_{B,K}(1/4096) / C_fail^{B,K}` lies in `[1.473, 1.478] x 10^5` and + `Gamma_2 / C_fail^{B,K}` in `[1.125, 1.129] x 10^5`. These are rounded outward from the certified `C_fail` interval; + `Gamma_2 = 307.51` is itself a floating value, so these are comparisons, not certified ratios. +- `d = 3`: + - likewise `C_{B,K}^(3)(1/4096) / C_fail^(3),{B,K}` lies in `[2.729, 2.737] x 10^5`, and + `Gamma_3 / C_fail^(3),{B,K}` in `[1.224, 1.228] x 10^5` (outward; `Gamma_3 = 540.97` is floating); + - Math-#184's floating `a_fail^(3)(0, k)` (`4.8558` at `k = 1/2`, `5.705` at `k = 1`) agree with the certified values to + `0.003%` and `0.008%`. +- Lemma T: direct floating lattice sums of the torus deviation of six one-site covariance entries at `L = 10`, in three + frames, are at most `3.3e-16`, below the certified `delta = 2.86e-14`. + +## 10. What this does not do + +- **Torus sides `L >= 10` only.** Lemma T covers every side `L >= 10` and every frame; `L < 10` is not covered. Corollary 2 + (`d = 3`) is for the reference kernel only. +- **Conditional identification.** + - The certified numbers are the Gaussian integrals that [NUM] (from [CL] Theorem N, [CUB] Theorems C and F) and + [S] §9 define. + - That `a_fail` is the selector-failure coefficient and `C_fail^{B,K}` the lifetime constant are statements of + [S]: AUTHOR_SIDE / HOLD, consumed, not reviewed. + - `nu(1) = alpha_1 + k integral_X Lambda` also contains the remote part of [RM], which is not evaluated here. + - Corollary 2 is conditional in addition on [FIB] (Math-#175, open, AUTHOR_SIDE / HOLD) and [HD] (Math-#184, + open). Its pointwise values are given only at `b = 0`; `R_3(1)` would need `m_(3,1)`, which is not certified here. +- **Coverage.** + - Pointwise coefficients only at `k = 1/2, 1, 2` and `b = 0, 1`, plus the band integral `C_fail^{B,K}`. + - No other band and no `k -> 0`. `d = 3` is covered only through Corollary 2, and there is nothing for `d >= 4`. +- No rate, no finite-`r` statement. +- Not a review of [NUM], [S], [CL], [CUB], [LP] or Math-#215. + +## 11. Provenance + +- **Sources:** `SOURCE_MAP.json` pins [NUM], [S], [LP], [CL], [CUB] and the [NUM] results on `main` `3e0a91b`. + Math-#215 (`d748c2a`, open) is pinned for comparison constants only. Math-#175 (`701dff4`) and Math-#184 (`d37ff5d`), + both open, are pinned for Corollary 2. +- **Packet files:** + - `certificate.py` (drivers, exact parts, tails, assembly, self-test, replay), `boxes.py` (roots, Taylor boxes), + `torus.py` (Lemma T) and `ia.py` (interval arithmetic), all standard-library Python; + - `controls.py` (floating controls); + - `RESULTS.json`, the certificate's output. +- **Manifest:** `SOURCE_FILES.json`. diff --git a/frontiers/c8_elder_failure_coefficient_20261001/README.md b/frontiers/c8_elder_failure_coefficient_20261001/README.md new file mode 100644 index 000000000..68ac15054 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/README.md @@ -0,0 +1,30 @@ +# Certified planar near cluster coefficients and the elder-failure constant + +Certified interval enclosures, for the reference planar kernel, of the near cluster coefficients `alpha_1`, `alpha_2` +(`nu(2) = alpha_2`) of [NUM] at `k = 1/2, 1, 2`, `b = 0, 1`, and of the constant `C_fail^{B,K}` of [S] §9 on +`[0, 1] x [1/2, 2]`. Details and the method are in [NOTE.md](NOTE.md). + +- `J_fail(1/2)` is in `[36.783058, 36.786801]`, `J_fail(1)` in `[172.84636, 172.86021]` and `J_fail(2)` in + `[1111.9822, 1112.0308]`. +- `J_2(1/2)` is in `[1.1061, 1.1095]`, `J_2(1)` in `[3.2531, 3.2625]` and `J_2(2)` in + `[8.8250, 8.8527]`. The last settles the Gauss–Hermite / Monte Carlo disagreement that [NUM] reported at + `k = 2`. +- `C_fail^{B,K}` is in `[0.0027244, 0.0027312]`, about `1.5 x 10^5` below the cap-route constant + `C_{B,K}(1/4096) = 402.53` of Math-#215. +- In `d = 3`, conditionally on Math-#175 and Math-#184, `C_fail^(3),{B,K}` is in `[0.0044058, 0.0044168]`. + +The planar values hold for the reference kernel and, by a torus-transfer lemma, for the torus model of [LP] and [S] with +every side `L >= 10` and every frame. + +**Method.** A third-order box scheme: +- second-order Taylor expansion in `(B, w = |u|)`, where the Hessian of the failure integrand is bounded and explicit; +- interval Hessian ranges times exact Gaussian moments, intersected with first-order monotone bounds; +- `k` as a fourth box coordinate for `C_fail`. + +**Arithmetic.** Standard-library Python and binary64 interval arithmetic; no transcendental library value is trusted. + +**Replay.** `python3 -B -S certificate.py --check` (CI in both interpreter modes, with five mutants). Floating controls: +`python3 -B certificate.py --controls`. + +**Scientific effect:** NONE. Author-side; conditional on the sources' identifications ([NUM] from [CL]/[CUB], [S]). +Same GitHub account as every lane; zero organizational-independence credit. diff --git a/frontiers/c8_elder_failure_coefficient_20261001/RESULTS.json b/frontiers/c8_elder_failure_coefficient_20261001/RESULTS.json new file mode 100644 index 000000000..670297d08 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/RESULTS.json @@ -0,0 +1,533 @@ +{ + "object": "CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1", + "scientific_effect": "NONE", + "certified": true, + "arithmetic": "binary64 interval arithmetic, outward one-ulp widening after every operation; published intervals widened by 1e-07 relatively and rounded outward to 10 significant digits", + "parameters": { + "k": [ + 0.5, + 1.0, + 2.0 + ], + "tol3": { + "V": 1e-08, + "H2": 1e-08 + }, + "tol4": 2.5e-09, + "box_half_width_sigma": 8.0, + "initial_cells": 8, + "initial_k_cells": 8, + "k_split_scale": 0.8, + "k_elementary_level": 12 + }, + "leaves": { + "V k=0.5": 1209983, + "H2 k=0.5": 1131985, + "V k=1": 1096831, + "H2 k=1": 799321, + "V k=2": 951074, + "H2 k=2": 583149, + "I_k 4D": 997266 + }, + "results": { + "k=0.5": { + "J_fail": [ + "36.78305814", + "36.78680021" + ], + "J_1": [ + "35.67360980", + "35.68065110" + ], + "J_2": [ + "1.106148894", + "1.109448563" + ], + "share_J2_over_Jfail": [ + "0.03006918240", + "0.03016194161" + ], + "near_palm_excess": [ + "0.05837775723", + "0.05856276965" + ], + "E[2T] exact": "199/9", + "parts": { + "E[2(B_-)^3]": [ + "8.357431233", + "8.357432938" + ], + "E[2T]": [ + "22.11110889", + "22.11111333" + ], + "E[V] boxes": [ + "18.57553740", + "18.58052069" + ], + "E[H_2] boxes": [ + "1.474865193", + "1.479264750" + ], + "tail bound": [ + "0", + "1.475969601E-9" + ] + }, + "coefficients": { + "b=0": { + "p_b(0)/z_0": [ + "0.03134386261", + "0.03134386889" + ], + "alpha_1": [ + "1.118148836", + "1.118369538" + ], + "alpha_2": [ + "0.03467098246", + "0.03477440682" + ], + "a_fail = alpha_1 + alpha_2": [ + "1.152923236", + "1.153040528" + ] + }, + "b=1": { + "p_b(0)/z_0": [ + "0.008974028845", + "0.008974030641" + ], + "alpha_1": [ + "0.3201360354", + "0.3201992242" + ], + "alpha_2": [ + "0.009926613082", + "0.009956224403" + ], + "a_fail = alpha_1 + alpha_2": [ + "0.3300922578", + "0.3301258393" + ] + } + } + }, + "k=1": { + "J_fail": [ + "172.8463609", + "172.8602041" + ], + "J_1": [ + "169.5839525", + "169.6070816" + ], + "J_2": [ + "3.253121786", + "3.262409093" + ], + "share_J2_over_Jfail": [ + "0.01881938142", + "0.01887461645" + ], + "near_palm_excess": [ + "0.03694156081", + "0.03705188358" + ], + "E[2T] exact": "1355/36", + "parts": { + "E[2(B_-)^3]": [ + "5.876020762", + "5.876021946" + ], + "E[2T]": [ + "37.63888512", + "37.63889266" + ], + "E[V] boxes": [ + "14.10054776", + "14.10515341" + ], + "E[H_2] boxes": [ + "1.084373928", + "1.087469698" + ], + "tail bound": [ + "0", + "3.232088888E-10" + ] + }, + "coefficients": { + "b=0": { + "p_b(0)/z_0": [ + "0.007835965654", + "0.007835967222" + ], + "alpha_1": [ + "1.328854160", + "1.329035399" + ], + "alpha_2": [ + "0.02549135313", + "0.02556412816" + ], + "a_fail = alpha_1 + alpha_2": [ + "1.354418283", + "1.354526758" + ] + }, + "b=1": { + "p_b(0)/z_0": [ + "0.002243507211", + "0.002243507661" + ], + "alpha_1": [ + "0.3804628584", + "0.3805147487" + ], + "alpha_2": [ + "0.007298402918", + "0.007319239058" + ], + "a_fail = alpha_1 + alpha_2": [ + "0.3877820960", + "0.3878131531" + ] + } + } + }, + "k=2": { + "J_fail": [ + "1111.982273", + "1112.030781" + ], + "J_1": [ + "1103.129655", + "1103.205691" + ], + "J_2": [ + "8.825088654", + "8.852620061" + ], + "share_J2_over_Jfail": [ + "0.007936011924", + "0.007961115375" + ], + "near_palm_excess": [ + "0.01574666826", + "0.01579686122" + ], + "E[2T] exact": "20957/288", + "parts": { + "E[2(B_-)^3]": [ + "5.093426722", + "5.093427745" + ], + "E[2T]": [ + "72.76735383", + "72.76736839" + ], + "E[V] boxes": [ + "14.80440892", + "14.80843560" + ], + "E[H_2] boxes": [ + "0.7354240545", + "0.7377183383" + ], + "tail bound": [ + "0", + "2.004587614E-10" + ] + }, + "coefficients": { + "b=0": { + "p_b(0)/z_0": [ + "0.001958991413", + "0.001958991806" + ], + "alpha_1": [ + "2.161021739", + "2.161170692" + ], + "alpha_2": [ + "0.01728827462", + "0.01734220843" + ], + "a_fail = alpha_1 + alpha_2": [ + "2.178363944", + "2.178458970" + ] + }, + "b=1": { + "p_b(0)/z_0": [ + "0.0005608768028", + "0.0005608769151" + ], + "alpha_1": [ + "0.6187198961", + "0.6187625424" + ], + "alpha_2": [ + "0.004949788004", + "0.004965229734" + ], + "a_fail = alpha_1 + alpha_2": [ + "0.6236851249", + "0.6237123314" + ] + } + } + }, + "d=3 (conditional on Math-#175 (L1) and Math-#184 Theorem 1)": { + "R_3(0) = (32 + 28 sqrt2)/17": [ + "4.211645446", + "4.211646289" + ], + "const_3": [ + "0.006583316240", + "0.006583317557" + ], + "I_b^(3)": [ + "1.513649011", + "1.513683082" + ], + "C_fail^(3),{B,K}": [ + "0.004405817806", + "0.004416745237" + ], + "k=0.5, b=0": { + "alpha_1^(3)": [ + "4.709246927", + "4.710176443" + ], + "alpha_2^(3) = nu^(3)(2)": [ + "0.1460219000", + "0.1464574868" + ], + "a_fail^(3)": [ + "4.855704384", + "4.856198372" + ] + }, + "k=1, b=0": { + "alpha_1^(3)": [ + "5.596663133", + "5.597426446" + ], + "alpha_2^(3) = nu^(3)(2)": [ + "0.1073605521", + "0.1076670548" + ], + "a_fail^(3)": [ + "5.704330166", + "5.704787020" + ] + }, + "k=2, b=0": { + "alpha_1^(3)": [ + "9.101458278", + "9.102085610" + ], + "alpha_2^(3) = nu^(3)(2)": [ + "0.07281209039", + "0.07303924044" + ], + "a_fail^(3)": [ + "9.174497503", + "9.174897715" + ] + } + }, + "torus model, every L >= 10, every frame (Lemma T)": { + "jet covariance deviation delta (order <= 6)": [ + "0", + "2.850266693E-14" + ], + "chi^2 bound": [ + "0", + "2.878353120E-23" + ], + "|J^L - J| bound (all k, J_fail and J_2)": [ + "0", + "1.247240771E-7" + ], + "(pi_0 p_b(0))^L / (pi_0 p_b(0))^ref": [ + "0.9999998999", + "1.000000101" + ], + "k=0.5": { + "J_fail": [ + "36.78305802", + "36.78680034" + ], + "J_2": [ + "1.106148770", + "1.109448688" + ], + "b=0": { + "p_b(0)/z_0": [ + "0.03134386261", + "0.03134386889" + ], + "alpha_1": [ + "1.118148828", + "1.118369546" + ], + "alpha_2": [ + "0.03467097855", + "0.03477441073" + ], + "a_fail": [ + "1.152923232", + "1.153040532" + ] + }, + "b=1": { + "p_b(0)/z_0": [ + "0.008974028845", + "0.008974030641" + ], + "alpha_1": [ + "0.3201360332", + "0.3201992265" + ], + "alpha_2": [ + "0.009926611963", + "0.009956225523" + ], + "a_fail": [ + "0.3300922567", + "0.3301258404" + ] + } + }, + "k=1": { + "J_fail": [ + "172.8463608", + "172.8602042" + ], + "J_2": [ + "3.253121661", + "3.262409218" + ], + "b=0": { + "p_b(0)/z_0": [ + "0.007835965654", + "0.007835967222" + ], + "alpha_1": [ + "1.328854158", + "1.329035401" + ], + "alpha_2": [ + "0.02549135216", + "0.02556412914" + ], + "a_fail": [ + "1.354418282", + "1.354526759" + ] + }, + "b=1": { + "p_b(0)/z_0": [ + "0.002243507211", + "0.002243507661" + ], + "alpha_1": [ + "0.3804628579", + "0.3805147493" + ], + "alpha_2": [ + "0.007298402638", + "0.007319239338" + ], + "a_fail": [ + "0.3877820958", + "0.3878131534" + ] + } + }, + "k=2": { + "J_fail": [ + "1111.982273", + "1112.030781" + ], + "J_2": [ + "8.825088529", + "8.852620186" + ], + "b=0": { + "p_b(0)/z_0": [ + "0.001958991413", + "0.001958991806" + ], + "alpha_1": [ + "2.161021738", + "2.161170692" + ], + "alpha_2": [ + "0.01728827438", + "0.01734220867" + ], + "a_fail": [ + "2.178363944", + "2.178458970" + ] + }, + "b=1": { + "p_b(0)/z_0": [ + "0.0005608768028", + "0.0005608769151" + ], + "alpha_1": [ + "0.6187198960", + "0.6187625426" + ], + "alpha_2": [ + "0.004949787934", + "0.004965229804" + ], + "a_fail": [ + "0.6236851248", + "0.6237123315" + ] + } + }, + "C_fail^{B,K}": [ + "0.002724454687", + "0.002731151043" + ] + }, + "C_fail": { + "I_k": [ + "0.4421366818", + "0.4432233938" + ], + "I_k boxes": [ + "0.4421366818", + "0.4432233938" + ], + "I_k tail bound": [ + "0", + "1.296653887E-11" + ], + "int_K W(k) dk": [ + "0.009075577556", + "0.009075579412" + ], + "const": [ + "0.008250963314", + "0.008250964965" + ], + "I_b": [ + "0.7468240581", + "0.7468242075" + ], + "C_fail^{B,K}": [ + "0.002724454695", + "0.002731151035" + ] + } + } +} diff --git a/frontiers/c8_elder_failure_coefficient_20261001/SOURCE_FILES.json b/frontiers/c8_elder_failure_coefficient_20261001/SOURCE_FILES.json new file mode 100644 index 000000000..ce8e98434 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/SOURCE_FILES.json @@ -0,0 +1,63 @@ +{ + "schema": 1, + "scientific_effect": "NONE", + "certified": true, + "review_status": "AUTHOR_SIDE_CERTIFIED_NUMERICS", + "object": "CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1", + "files": [ + { + "path": "NOTE.md", + "bytes": 27438, + "sha256": "de50468afe0cae0bd2147a79e1b5b45a9f57ade5a81a9a99fd95aaaccc826c7c", + "git_blob": "9903dd6f6b2514069bef78399b3d90201280911e" + }, + { + "path": "README.md", + "bytes": 1906, + "sha256": "63313f05a155a5b93db70c82c16eca7ab98c4f1311ddc770dffc5ccb2db85031", + "git_blob": "68ac15054b26ab4e0b9021859727be3385b6dbbd" + }, + { + "path": "RESULTS.json", + "bytes": 9040, + "sha256": "73bd816e0b48d9740c02d58f32fda99a110e4634222760e17c9c98df64b4d870", + "git_blob": "670297d08e516b5c3456fef278abb8ce454cfd79" + }, + { + "path": "SOURCE_MAP.json", + "bytes": 6451, + "sha256": "ab8a34556de2c309d5a379bd17b047159622dc111bf78ecf818f6dd9e339c4c9", + "git_blob": "f48177918014d3b09fe40e810695c986b1a30af5" + }, + { + "path": "boxes.py", + "bytes": 28396, + "sha256": "5d6a9cfca33ed5140065bfd02a630a77596adf8c0b45ed6edb0032005754f9db", + "git_blob": "f6b42e2a0bf0021d85ded2bd19a77205455a5aa0" + }, + { + "path": "certificate.py", + "bytes": 30759, + "sha256": "a0027c937e74c41e746a225170faca3bc646358e0d1873be2cff103ec4cedd21", + "git_blob": "21318a7bba41a16b92ae7f72a34d0bbcbad62c9a" + }, + { + "path": "controls.py", + "bytes": 5593, + "sha256": "d275ee5a120c3246e0c1a6dc5cd8bc8fbe738069c64448894efe3feebb7dd1f1", + "git_blob": "39613832228583a9eb81d3eb2a1ee12e2b791d39" + }, + { + "path": "ia.py", + "bytes": 8930, + "sha256": "6d048d58c6ab18f21c0609fc888ca73497a60e234217bb4f7a309a1502fccc9e", + "git_blob": "cf4bbfa16f577fe65f4ab5fb3621943fa63cf9fa" + }, + { + "path": "torus.py", + "bytes": 9778, + "sha256": "c5e3496d62dd7dcb544a72722982221654b9e6de100d5cf36640ca4abd67d279", + "git_blob": "4a448b2f30e8f1157b758bc9ae5c9c258e112db8" + } + ] +} diff --git a/frontiers/c8_elder_failure_coefficient_20261001/SOURCE_MAP.json b/frontiers/c8_elder_failure_coefficient_20261001/SOURCE_MAP.json new file mode 100644 index 000000000..f48177918 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/SOURCE_MAP.json @@ -0,0 +1,111 @@ +{ + "schema": 1, + "object": "CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1", + "scientific_effect": "NONE", + "source_repository": "d6g8k5htny-coder/Math-", + "main_at_creation": "3e0a91b69892b9717564c592e7a966492a81566a", + "note": "Certified interval enclosures of explicit Gaussian integrals that [NUM] (from [CL]/[CUB]) and [S] define, for the reference (continuum) kernel. The identifications of alpha_1, alpha_2, a_fail and C_fail^{B,K} with model quantities are those of the sources and are not re-proved or accepted here. [LDC] is an open pull request, used only for comparison constants. The d = 3 corollary is conditional in addition on [FIB] (Math-#175) and [HD] (Math-#184), both open.", + "sources": [ + { + "path": "frontiers/c6_cluster_coefficients_numerics_20260930/NOTE.md", + "bytes": 11018, + "sha256": "f20c1ac96b00edc071e428bd2e5dc4c91b8546a3ee0693367b783109b26aa518", + "git_blob": "845881d40520d4acbe15ee196a1ba0d5412eafa1", + "id": "NUM", + "commit": "e27b4c80a4f2b11227a7139f819b3e9ba4de8493", + "role": "section 2: the jet covariances, the even/odd regression, A ~ N(-b, 2), p_b(0) = phi(b/sqrt2)/sqrt2, z_0 = 36 k^2 m_(2,b), (a, beta, c) ~ N(0, diag(2, 2, 6)) and alpha_j = [phi(b/sqrt2)/(sqrt2 36 k^2 m_(2,b))] J_j(k); section 3: B, D, the typed domain s < -|B|/2, the weight 9k^2 (4s^2 - B^2), the classifier (n = 2 iff s > B and g(s) > 0; n = 1 iff g(s) < 0) and I_j; the floating values compared in section 9", + "current_required": true, + "on_main": true + }, + { + "path": "frontiers/local_elder_geometry_20260930/PROOF.md", + "bytes": 39722, + "sha256": "f68038be79b46124b0f9b31205aa6e3682b34b6f5ef81f4a5697ae81f07cc46b", + "git_blob": "ef2aa57959ea9f721bbf2316ce94cf616c1c9113", + "id": "S", + "commit": "79f18f0cb318fceccd5bdf4e6dad29815165bc7d", + "role": "Theorem S, (S2): r^-3 (1 - p_r) -> a_fail = alpha_1 + alpha_2; section 9: C_fail^{B,K} = int_{B x K x S^1} A_0 a_fail/(3k^(5/3)), nu_rej ~ C_fail ell^(2/3), E N_rej(0, t] ~ (3/5) C_fail t^(5/3). AUTHOR_SIDE / HOLD in its own README; consumed as a statement, not reviewed", + "current_required": true, + "on_main": true + }, + { + "path": "imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md", + "bytes": 40261, + "sha256": "9350ad6eaba6626b93c3dedeef9e2ff816e5cdf1c8318e85fb27499141c84bc7", + "git_blob": "dfed3b8d318a3ab1950957f393307733a4bef3f2", + "id": "LP", + "commit": "b6282ecc44cf63664a989efd055e473aa1abc50b", + "role": "section 2: the kernel; (10.2)-(10.3): A_r = 12 pi_r Z_r / r^2 -> A_0 = 12 pi_0 z_0; section 15: pi_0 = p_(f,V_u)(b, 0) p_G(0) phi_tau_u(12k), z_0 = (6k)^2 E[(det A_u)^2 1{A_u < 0} | f = b, V_u = 0]", + "current_required": true, + "on_main": true + }, + { + "path": "frontiers/c6_cluster_law_20260929/PROOF.md", + "bytes": 70131, + "sha256": "96b8e2d59defee6059bf151a896b7fc7c91fc529b385a3db8a87574901197b38", + "git_blob": "ba492c8e62e58bfc055fc8254c790e893d346ba5", + "id": "CL", + "commit": "4807bfffb880ce9e4ae3b04a3faf0dd777a9945c", + "role": "Theorem N, (1.3) and (3.11): the near coefficients alpha_1, alpha_2 and nu(2) = alpha_2 (through [NUM]); nonauthor-accepted (OpenAI review 5360192822) and merged", + "current_required": true, + "on_main": true + }, + { + "path": "frontiers/planar_cubic_cluster_20260929/PROOF.md", + "bytes": 19889, + "sha256": "117e9299a71e6139270266889eb772ed97518e402d9af029bbdb67caf9d4a0f0", + "git_blob": "bb446d08db8a944537a743ad550b88c1c2ad5758", + "id": "CUB", + "commit": "325c5ebdc2ca0165b483cddce8972ee2d0252c79", + "role": "Theorem C (the classifier of the sheared cubic), Theorem F (alpha_1, alpha_2 finite and positive), Theorem B (parity factorization), through [NUM]", + "current_required": true, + "on_main": true + }, + { + "path": "frontiers/c6_cluster_coefficients_numerics_20260930/RESULTS.json", + "bytes": 4095, + "sha256": "9c88ea85b6ba08199996c9b5ff74ecd06afc0f458e08c28b4c3af6191f437fcb", + "git_blob": "2648314f8c11f0637ce301ba76360142a8309e8f", + "id": "NUM-R", + "commit": "0b35e76ff05676e6ad5c89171d65807668c8de15", + "role": "floating Gauss-Hermite and Monte Carlo values of J_1, J_2; read only by controls.py for the comparison table (not part of the certificate)", + "current_required": true, + "on_main": true + }, + { + "path": "frontiers/c8_lifetime_difference_constant_20261001/NOTE.md", + "bytes": 30864, + "sha256": "d77c9ab1792f1cba55c512828ad0601216068167b56bc4c912f9e8e7d985ae3e", + "git_blob": "0f6c7fb79988b52e26c1e2f1456ab196f075ecf0", + "id": "LDC", + "commit": "d748c2abee2cbfbad6504b5700c70bbd83ba5785", + "role": "Math-#215 (open, not on main): C_{B,K}(1/4096) = 402.5297, Gamma_2 = 307.51 (floating) and Theorem L (eps_2 <= 3.96e-10); comparison and scope remarks only", + "current_required": false, + "on_main": false + }, + { + "path": "frontiers/concave_fibre_elder_20260930/PROOF.md", + "bytes": 23416, + "sha256": "5681896a6808f48fb7508914d0860ec508af7646e9e706fc260d935dd951c75e", + "git_blob": "923d3236b5291a29bd162937f172190db84dcef5", + "id": "FIB", + "commit": "701dff49119897932ebe1891657126eafe395711", + "role": "Math-#175 (OpenAI/Codex, open, AUTHOR_SIDE/HOLD): Theorem F (S3), (1 - p_r)/r^3 -> a1 + a2 in fixed d >= 3, and (L1), C_fail^{B,K} = int_{B x K x S^(d-1)} A0 a_fail/(3k^(5/3)); consumed as statements for the d = 3 corollary only", + "current_required": false, + "on_main": false + }, + { + "path": "frontiers/c6_hard_direction_factorization_20260930/PROOF.md", + "bytes": 22456, + "sha256": "26743f82bf4b86d4ae4c7dfb085d4d15bf36014aa627dee5edfc86a75ca64c05", + "git_blob": "91d22fce0fd7d4630bfe2e2c93eab289babc97f2", + "id": "HD", + "commit": "d37ff5d0315a4f71453d37f0cae2da1f40fdef3c", + "role": "Math-#184 (Anthropic/Claude, open; OpenAI nonauthor ACCEPTs 5366497205 (Slices A/B) and 5366261768 (Slice C)): Theorem 1 and Corollary 1, a_j^(d) = R_d(b) a_j^(2), R_d = N_d m_(2,b)/m_(d,b), N_3(b) = pi (2 pi)^(-1/2) int_0^inf h^3 phi_2(h - b) dh, R_3(0) = (32 + 28 sqrt2)/17; consumed for the d = 3 corollary only", + "current_required": false, + "on_main": false + } + ], + "nonauthor_acceptance": false, + "governing_status_changed": false +} diff --git a/frontiers/c8_elder_failure_coefficient_20261001/boxes.py b/frontiers/c8_elder_failure_coefficient_20261001/boxes.py new file mode 100644 index 000000000..f6b42e2a0 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/boxes.py @@ -0,0 +1,734 @@ +"""Box enclosures of the failure integrands (certified engine; standard library only). + +Jet coordinates. (a, beta, c) ~ N(0, diag(2, 2, 6)) are the odd jets of [NUM] section 2. For a gap mark k the scaled +coordinates a^ = a / sqrt k ~ N(0, 2/k), beta ~ N(0, 2), c^ = sqrt(k) c ~ N(0, 6k) give + B = beta - a^2/12, u = c^ + d, d = -a^ beta/4 + a^3/72, w = |u|, T = 3 w^2 ( = 12 k D^2 of [NUM] section 3). +Roots. P(x) = (x + B)^2 (2x - B). The big root x(B, T) >= xm = max(-B, B/2) solves P(x) = T (P increasing and convex on +[xm, inf), P(xm) = 0). For B < 0 and T <= |B|^3/2 the middle root solves P(x) = T in [|B|/2, |B|] (P decreasing there). +Integrands. With r = B/x: + H = |B|^3 + 4x^3 - 3B^2 x = 2(B_-)^3 + 2T + V, V = -3B(2x - B)(x + B) (big root) + H_2 = |B|^3 + 4x^3 - 3B^2 x = 2|B|^3 + 2T + V (middle root, inside its support) + V_B = x^2 (-6 + 9r^2/2 - 3r^3/2), V_w = -sigma sqrt3 x^(3/2) r (4 + r)(2 - r)^(1/2) + V_BB = x p(r), p(r) = 6 + 3r - 9r^2/2 - 3r^3/4 + 3r^4/4 + V_Bw = -sigma (sqrt3/2) x^(1/2) e(r), e(r) = (2 - r)^(3/2)(4 + 2r + r^2) + V_ww = g(r) = -r(4 - r + r^2) +with sigma = +1 on the big root and -1 on the middle root (derivatives in (B, w) at fixed w, resp. fixed B). +Per box: I1 = Gaussian mass x pointwise range (monotonicity at the corners of the image rectangle), and, on one side of +u = 0, I2 = second-order Taylor expansion in (B, w) about a float point (B*, w*) of the rectangle with exact Gaussian +polynomial moments and the integral-form remainder, whose averaged Hessian entries lie in the interval ranges of the +second derivatives over the rectangle. The box integral lies in I1 and in I2; the intersection is returned (an empty +intersection raises).""" +import math +from ia import (dn, up, add, sub, neg, mul, div, scal, divc, sqr, pw, isqrt, exp_neg, Q_point, SQRT2PI) + +MUTANT = None # set by certificate.py --mutant (the certificate must then fail) +SQRT3 = isqrt((3.0, 3.0)) +HSQRT3 = scal(0.5, SQRT3) +_U = 2.0 ** -53 + + +# ================================================================ 1D Gaussian partial moments +def phis(x, s, s2): + """phi_sigma(x) for a float x, sigma interval s, sigma^2 interval s2.""" + return div(exp_neg(div(sqr((x, x)), scal(2.0, s2))), mul(s, SQRT2PI)) + + +def _Qiv(t): + """Q = 1 - Phi over an interval t >= 0 (decreasing).""" + return (Q_point(t[1])[0], Q_point(t[0])[1]) + + +def cell_mass(x0, x1, s): + """Gaussian mass of [x0, x1], computed on the side away from 0 as a difference of upper tails.""" + t0, t1 = div((x0, x0), s), div((x1, x1), s) + if x0 >= 0.0: + return sub(_Qiv(t0), _Qiv(t1)) + if x1 <= 0.0: + return sub(_Qiv(neg(t1)), _Qiv(neg(t0))) + return sub(sub((1.0, 1.0), _Qiv(t1)), _Qiv(neg(t0))) + + +def moments(x0, x1, s, s2, nmax): + """M_n = int_{x0}^{x1} x^n phi_sigma(x) dx, n = 0..nmax: M_1 = s^2 (phi(x0) - phi(x1)), + M_n = (n - 1) s^2 M_{n-2} + s^2 (x0^(n-1) phi(x0) - x1^(n-1) phi(x1)).""" + p0, p1 = phis(x0, s, s2), phis(x1, s, s2) + m0 = cell_mass(x0, x1, s) + M = [(max(0.0, m0[0]), m0[1])] + if nmax >= 1: + M.append(mul(s2, sub(p0, p1))) + x0i, x1i = (x0, x0), (x1, x1) + for n in range(2, nmax + 1): + t = sub(mul(pw(x0i, n - 1), p0), mul(pw(x1i, n - 1), p1)) + M.append(add(mul(scal(float(n - 1), s2), M[n - 2]), mul(s2, t))) + return M + + +class Moments: + """Cache of the partial moments of one coordinate.""" + def __init__(self, sig, sig2, nmax): + self.sig, self.sig2, self.nmax, self.c = sig, sig2, nmax, {} + + def get(self, x0, x1): + key = (x0, x1) + r = self.c.get(key) + if r is None: + r = moments(x0, x1, self.sig, self.sig2, self.nmax) + self.c[key] = r + return r + + +# ================================================================ verified roots +def _sgn_f(x, B, T): + """Rigorous sign of P(x) - T at floats x, B, T >= 0 (+1, -1, or 0 if undecided). The binary64 error of + p = x + B, q = 2x - B (2x exact), p*p*q and the final subtraction is below 5u (p^2 |q|) + u |f|; the test uses + 8u (p^2 |q| + T) + 4u |f| + 1e-300.""" + p = x + B + q = 2.0 * x - B + r = p * p * q + f = r - T + e = 8.0 * _U * (abs(r) + T) + 4.0 * _U * abs(f) + 1e-300 + if f > e: + return 1 + if f < -e: + return -1 + return 0 + + +def _guess_big(B, T): + aB = abs(B) + if B < 0 and T < aB ** 3: + c = -1.0 + 2.0 * T / aB ** 3 + ph = math.acos(max(-1.0, min(1.0, c))) + return aB * (math.cos(ph / 3.0) + 0.5) + h = B ** 3 / 8.0 + T / 4.0 + sd = math.sqrt(max(0.0, T * (B ** 3 + T) / 16.0)) + A, C = h + sd, h - sd + y = math.copysign(abs(A) ** (1.0 / 3.0), A) + math.copysign(abs(C) ** (1.0 / 3.0), C) + return y - B / 2.0 + + +def _guess_mid(B, T): + aB = -B + c = -1.0 + 2.0 * T / aB ** 3 + ph = math.acos(max(-1.0, min(1.0, c))) + return aB * (math.cos((2.0 * math.pi - ph) / 3.0) + 0.5) + + +def _newton_big(B, T, xm): + lo = xm + hi = xm + 2.0 * (T / 2.0) ** (1.0 / 3.0) + abs(B) + 1.0 # starting bracket only (the result is verified) + x = _guess_big(B, T) + if not (lo < x < hi): + x = 0.5 * (lo + hi) + for _ in range(60): + fx = (x + B) * (x + B) * (2 * x - B) - T + d = 6 * x * (x + B) + if fx > 0: + hi = min(hi, x) + else: + lo = max(lo, x) + nx = x - fx / d if d > 0 else 0.5 * (lo + hi) + if not (lo < nx < hi): + nx = 0.5 * (lo + hi) + if abs(nx - x) <= 2e-15 * max(1.0, abs(x)): + return nx + x = nx + return x + + +def xbig(B, Tlo, Thi=None): + """Interval containing the big root x(B, T) for every T in [Tlo, Thi]: the lower end a has a = xm or P(a) < Tlo, the + upper end b has P(b) > Thi (P increasing on [xm, inf)).""" + if Thi is None: + Thi = Tlo + xm = -B if B < 0 else B / 2.0 + if Thi <= 0.0: + return (xm, xm) + x = _newton_big(B, 0.5 * (Tlo + Thi), xm) + delta = 8e-16 * max(1.0, abs(x)) + for _ in range(60): + a, b = max(xm, x - delta), x + delta + if (a <= xm or _sgn_f(a, B, Tlo) == -1) and _sgn_f(b, B, Thi) == 1: + return (a, b) + delta *= 4.0 + raise RuntimeError('big-root bracket failed: B=%r T=[%r, %r]' % (B, Tlo, Thi)) + + +def xmid(B, Tlo, Thi=None): + """Interval containing the middle root x(B, T) in [|B|/2, |B|] for every T in [Tlo, Thi] (B < 0, + Thi <= |B|^3/2): lower end a = |B|/2 or P(a) > Thi, upper end b = |B| or P(b) < Tlo (P decreasing there).""" + if Thi is None: + Thi = Tlo + aB = -B + lo0, hi0 = aB / 2.0, aB + if Thi <= 0.0: + return (hi0, hi0) + T = 0.5 * (Tlo + Thi) + lo, hi = lo0, hi0 + x = _guess_mid(B, T) + if not (lo < x < hi): + x = 0.5 * (lo + hi) + for _ in range(80): + fx = (x + B) * (x + B) * (2 * x - B) - T + if fx > 0: + lo = x + else: + hi = x + d = 6 * x * (x + B) + nx = x - fx / d if d < 0 else 0.5 * (lo + hi) + if not (lo < nx < hi): + nx = 0.5 * (lo + hi) + if abs(nx - x) <= 2e-15 * max(1.0, abs(x)): + x = nx + break + x = nx + delta = 8e-16 * max(1.0, abs(x)) + for _ in range(60): + a, b = max(lo0, x - delta), min(hi0, x + delta) + if (a <= lo0 or _sgn_f(a, B, Thi) == 1) and (b >= hi0 or _sgn_f(b, B, Tlo) == -1): + return (a, b) + delta *= 4.0 + raise RuntimeError('middle-root bracket failed: B=%r T=[%r, %r]' % (B, Tlo, Thi)) + + +def _cbrt_lo(y): + """A float c >= 0 with c^3 <= y (y >= 0).""" + if y <= 0: + return 0.0 + c = y ** (1.0 / 3.0) + while pw((c, c), 3)[1] > y: + c = dn(c) + return c + + +# ================================================================ pointwise pieces +def _nonneg(a): + return (max(0.0, a[0]), max(0.0, a[1])) + + +def _hull0(a): + return (min(0.0, a[0]), max(0.0, a[1])) + + +def V_iv(Bi, xi): + """V = -3B(2x - B)(x + B).""" + return mul(scal(-3.0, Bi), mul(sub(scal(2.0, xi), Bi), add(xi, Bi))) + + +def VB_iv(Bi, xi): + """V_B = -6x^2 + 9B^2/2 - 3B^3/(2x) (either root, x > 0).""" + return add(add(scal(-6.0, sqr(xi)), scal(4.5, sqr(Bi))), neg(div(scal(1.5, pw(Bi, 3)), xi))) + + +def Vw_big_iv(Bi, xi): + """V_w on the big root: -sqrt3 B (4x + B)(2x - B)^(1/2)/x (the middle root has the opposite sign).""" + q = _nonneg(sub(scal(2.0, xi), Bi)) + v = neg(div(mul(mul(SQRT3, mul(Bi, add(scal(4.0, xi), Bi))), isqrt(q)), xi)) + return neg(v) if MUTANT == 'vw-sign' else v + + +_P = (6.0, 3.0, -4.5, -0.75, 0.75) # p(r), increasing degree +_DP = (3.0, -9.0, -2.25, 3.0) # p'(r) + + +def _horner(coef, R): + acc = (coef[-1], coef[-1]) + for c in reversed(coef[:-1]): + acc = add(mul(acc, R), (c, c)) + return acc + + +def p_range(R): + """Mean-value enclosure of p over the interval R.""" + rc = 0.5 * (R[0] + R[1]) + return add(_horner(_P, (rc, rc)), mul(_horner(_DP, R), sub(R, (rc, rc)))) + + +def _g_pt(r): + ri = (r, r) + return neg(mul(ri, add(sub((4.0, 4.0), ri), sqr(ri)))) + + +def g_range(R): + """g(r) = -r(4 - r + r^2) is strictly decreasing (g' = -4 + 2r - 3r^2 < 0).""" + if MUTANT == 'ww-drop': + return (0.0, 0.0) + return (_g_pt(R[1])[0], _g_pt(R[0])[1]) + + +def _e_pt(r): + t = _nonneg(sub((2.0, 2.0), (r, r))) + ri = (r, r) + return mul(mul(t, isqrt(t)), add(add((4.0, 4.0), scal(2.0, ri)), sqr(ri))) + + +def e_range(R): + """e(r) = (2 - r)^(3/2)(4 + 2r + r^2) is decreasing on r <= 2 (e' = (2 - r)^(1/2)(-2 - r - 7r^2/2) < 0).""" + return (_e_pt(R[1])[0], _e_pt(R[0])[1]) + + +# ================================================================ image rectangles +def d_iv(ai, bi): + return add(divc(neg(mul(ai, bi)), 4.0), divc(pw(ai, 3), 72.0)) + + +def d_range(a0, a1, b0, b1): + """Range of d = -a beta/4 + a^3/72 over a rectangle: linear in beta (extremes at beta = b0, b1); for fixed beta + the a-extremes are at a0, a1 or at a = +-sqrt(6 beta) (beta > 0), where d = -+ beta^(3/2)/sqrt6.""" + cands = [] + for b in (b0, b1): + bi = (b, b) + for a in (a0, a1): + cands.append(d_iv((a, a), bi)) + if b > 0: + r = isqrt(scal(6.0, bi)) + val = div(mul(bi, isqrt(bi)), isqrt((6.0, 6.0))) + if r[0] < a1 and a0 < r[1]: + cands.append(neg(val)) + if -r[1] < a1 and a0 < -r[0]: + cands.append(val) + return (min(c[0] for c in cands), max(c[1] for c in cands)) + + +def box_geometry(box): + """(Bl, Bh, ul, uh): ranges of B and u over a box of scaled coordinates.""" + (a0, a1), (b0, b1), (c0, c1) = box + A2 = sqr((a0, a1)) + Bl = dn(b0 - up(A2[1] / 12.0)) + Bh = up(b1 - dn(A2[0] / 12.0)) + dl, dh = d_range(a0, a1, b0, b1) + return Bl, Bh, dn(c0 + dl), up(c1 + dh) + + +def wT(ul, uh): + """(s, wl, wh, Tl, Th): s = sign of u on the box (0 if the u-range meets 0), ranges of w = |u| and T = 3w^2.""" + if ul > 0.0: + s, wl, wh = 1, ul, uh + elif uh < 0.0: + s, wl, wh = -1, -uh, -ul + else: + s, wl, wh = 0, 0.0, max(-ul, uh) + Tl = max(0.0, dn(3.0 * max(0.0, dn(wl * wl)))) if wl > 0.0 else 0.0 + Th = up(3.0 * up(wh * wh)) + return s, wl, wh, Tl, Th + + +def half_cube(B): + """|B|^3/2 for a float B < 0.""" + return scal(0.5, pw((-B, -B), 3)) + + +# ================================================================ second-order pieces +def taylor_moments3(A, Bm, C, Bs, us): + """(m, int dB, int du, int dB^2, int du^2, int dB du) over a 3D box, dB = B - Bs, du = u - us (exact polynomial + moments; A to order 6, Bm and C to order 2).""" + ab = lambda i, j: mul(A[i], Bm[j]) + ab00, ab01, ab02 = ab(0, 0), ab(0, 1), ab(0, 2) + ab11, ab12 = ab(1, 1), ab(1, 2) + ab20, ab21, ab22 = ab(2, 0), ab(2, 1), ab(2, 2) + ab30, ab31 = ab(3, 0), ab(3, 1) + ab40, ab41 = ab(4, 0), ab(4, 1) + ab50, ab60 = ab(5, 0), ab(6, 0) + Bsi, usi = (Bs, Bs), (us, us) + gB0 = sub(sub(ab01, divc(ab20, 12.0)), mul(Bsi, ab00)) # int (B - Bs) over (a, beta) + gBB = add(add(add(ab02, divc(ab40, 144.0)), mul(sqr(Bsi), ab00)), + add(sub(divc(mul(Bsi, ab20), 6.0), divc(ab21, 6.0)), scal(-2.0, mul(Bsi, ab01)))) + Id1 = add(divc(neg(ab11), 4.0), divc(ab30, 72.0)) # int d + Id2 = add(add(divc(ab22, 16.0), divc(neg(ab41), 144.0)), divc(ab60, 5184.0)) # int d^2 + IBd = sub(add(add(divc(neg(ab12), 4.0), divc(scal(5.0, ab31), 144.0)), divc(neg(ab50), 864.0)), mul(Bsi, Id1)) + C0, C1, C2 = C[0], C[1], C[2] + Cs1 = sub(C1, mul(usi, C0)) # int (c - us) + Cs2 = add(sub(C2, scal(2.0, mul(usi, C1))), mul(sqr(usi), C0)) # int (c - us)^2 + m = mul(ab00, C0) + E1B = mul(gB0, C0) + E1u = add(mul(ab00, Cs1), mul(Id1, C0)) + QB = mul(gBB, C0) + Qu = add(add(mul(ab00, Cs2), scal(2.0, mul(Id1, Cs1))), mul(Id2, C0)) + QBu = add(mul(gB0, Cs1), mul(IBd, C0)) + return m, E1B, E1u, _nonneg(QB), _nonneg(Qu), QBu + + +def _quad(SBB, SBw, Sww, QB, Qu, QBu, s): + """Remainder integral 1/2 int (h_BB dB^2 + 2 h_Bw dB dw + h_ww dw^2) with averaged Hessian entries h in the ranges + S: in 1/2 (SBB QB + Sww Qu) + s c QBu +- d sqrt(QB Qu), c = centre and d = radius of SBw (dw = s du).""" + c = 0.5 * (SBw[0] + SBw[1]) + dBw = up(max(up(SBw[1] - c), up(c - SBw[0]))) + cross = scal(float(s), mul((c, c), QBu)) + rad = up(dBw * up(math.sqrt(up(QB[1] * Qu[1]))) * (1.0 + 4 * _U)) + q = add(add(scal(0.5, mul(SBB, QB)), scal(0.5, mul(Sww, Qu))), cross) + return (dn(q[0] - rad), up(q[1] + rad)) + + +def _expansion_point(ac, bc, cc, Bl, Bh, s, wl, wh): + """Float point (B*, w*) of the image rectangle near the image of the box centre (a^, beta, c^) = (ac, bc, cc).""" + Bs = min(max(bc - ac * ac / 12.0, Bl), Bh) + ws = min(max(s * (cc - ac * bc / 4.0 + ac ** 3 / 72.0), wl), wh) + return Bs, ws + + +def _intersect(I1, I2, what): + lo, hi = max(I1[0], I2[0]), min(I1[1], I2[1]) + if lo > hi: + raise RuntimeError('inconsistent %s enclosures: %r vs %r' % (what, I1, I2)) + return (lo, hi) + + +class Rect: + """Image rectangle of a box with its corner roots (big root).""" + def __init__(self, Bl, Bh, ul, uh): + self.Bl, self.Bh, self.ul, self.uh = Bl, Bh, ul, uh + self.s, self.wl, self.wh, self.Tl, self.Th = wT(ul, uh) + self._c = None + + def corners(self): + if self._c is None: + Bl, Bh, Tl, Th = self.Bl, self.Bh, self.Tl, self.Th + self._c = (xbig(Bl, Tl), xbig(Bl, Th), xbig(Bh, Tl), xbig(Bh, Th)) + return self._c + + +def V_range(R): + """[min V, max V] over the rectangle: V non-increasing in B; in T non-increasing for B > 0, non-decreasing for + B < 0; V(0, T) = 0.""" + xLL, xLH, xHL, xHH = R.corners() + Bl, Bh = R.Bl, R.Bh + if MUTANT == 'corner-swap': + xHH, xHL = xHL, xHH + vmin = V_iv((Bh, Bh), xHH)[0] if Bh > 0 else (V_iv((Bh, Bh), xHL)[0] if Bh < 0 else 0.0) + vmax = V_iv((Bl, Bl), xLH)[1] if Bl < 0 else (V_iv((Bl, Bl), xLL)[1] if Bl > 0 else 0.0) + return (vmin, vmax) + + +def _Bneg3(Bi): + lo = pw((max(0.0, -Bi[1]),) * 2, 3)[0] + hi = pw((max(0.0, -Bi[0]),) * 2, 3)[1] + return (max(0.0, lo), hi) + + +def _Bneg(Bi): + return (max(0.0, -Bi[1]), max(0.0, -Bi[0])) + + +def _H_at(B, T, xi): + return add(add(scal(2.0, _Bneg3((B, B))), (dn(2.0 * T), up(2.0 * T))), V_iv((B, B), xi)) + + +def H_range(R): + """[min H, max H]: H is non-increasing in B and non-decreasing in T.""" + xLL, xLH, xHL, xHH = R.corners() + return (max(0.0, _H_at(R.Bh, R.Tl, xHL)[0]), _H_at(R.Bl, R.Th, xLH)[1]) + + +def big_hessian(R, with_H): + """Interval ranges (S_BB, S_Bw, S_ww) of the (B, w)-Hessian of V (or of H if with_H) over the rectangle (one side of + u = 0). x is increasing in T and, for fixed T, unimodal in B with minimum x = B at B = (T/4)^(1/3); r = B/x is + non-decreasing in B, increasing in T for B < 0 and decreasing for B > 0; r in [-1, 2].""" + xLL, xLH, xHL, xHH = R.corners() + Bl, Bh, Tl = R.Bl, R.Bh, R.Tl + XL, XH = min(xLL[0], xHL[0]), max(xLH[1], xHH[1]) + if Bh > 0: + XL = min(XL, _cbrt_lo(max(0.0, dn(Tl / 4.0)))) + if XL <= 0.0: + Rr = (-1.0, 2.0) + XL = 0.0 + else: + rmin = div((Bl, Bl), xLL)[0] if Bl < 0 else (div((Bl, Bl), xLH)[0] if Bl > 0 else 0.0) + rmax = div((Bh, Bh), xHL)[1] if Bh > 0 else (div((Bh, Bh), xHH)[1] if Bh < 0 else 0.0) + Rr = (max(-1.0, rmin), min(2.0, rmax)) + X = (XL, XH) + SBB = mul(X, p_range(Rr)) + if with_H: + SBB = add(SBB, scal(12.0, _Bneg((Bl, Bh)))) + SBw = neg(mul(HSQRT3, mul(isqrt(X), e_range(Rr)))) + Sww = g_range(Rr) + if with_H: + Sww = add(Sww, (12.0, 12.0)) + return SBB, SBw, Sww + + +def big_center(Bs, ws, with_H): + """(value, d/dB, d/dw) of V (or H) at the float point (B*, w*), w* >= 0.""" + w2 = (max(0.0, dn(ws * ws)), up(ws * ws)) + T2 = (dn(3.0 * w2[0]), up(3.0 * w2[1])) + xs = xbig(Bs, T2[0], T2[1]) + if xs[0] <= 0.0: + return None + Bsi = (Bs, Bs) + v0, vB, vw = V_iv(Bsi, xs), VB_iv(Bsi, xs), Vw_big_iv(Bsi, xs) + if with_H: + v0 = add(add(scal(2.0, _Bneg3(Bsi)), scal(2.0, T2)), v0) + vB = sub(vB, scal(6.0, sqr(_Bneg(Bsi)))) + vw = add(vw, scal(12.0, (ws, ws))) + return v0, vB, vw + + +def second_order(R, Bs, ws, mom, with_H): + """Second-order enclosure (None if the expansion point has x = 0).""" + c = big_center(Bs, ws, with_H) + if c is None: + return None + v0, vB, vw = c + SBB, SBw, Sww = big_hessian(R, with_H) + m, E1B, E1u, QB, Qu, QBu = mom + lin = add(add(mul(v0, m), mul(vB, E1B)), scal(float(R.s), mul(vw, E1u))) + return add(lin, _quad(SBB, SBw, Sww, QB, Qu, QBu, R.s)) + + +# ================================================================ 3D box integrals (scaled coordinates, fixed k) +def _box3_common(box, Ma, Mb, Mc): + (a0, a1), (b0, b1), (c0, c1) = box + A, Bm, C = Ma.get(a0, a1), Mb.get(b0, b1), Mc.get(c0, c1) + R = Rect(*box_geometry(box)) + return A, Bm, C, R + + +def big_box3(box, Ma, Mb, Mc, with_H=False, second=True, parts=False): + """Enclosure of int_box V dmu (with_H: of H) over a 3D box of scaled coordinates.""" + A, Bm, C, R = _box3_common(box, Ma, Mb, Mc) + m0 = mul(mul(A[0], Bm[0]), C[0]) + I1 = mul(m0, H_range(R) if with_H else V_range(R)) + I2 = None + if second and R.s != 0: + (a0, a1), (b0, b1), (c0, c1) = box + Bs, ws = _expansion_point(0.5 * (a0 + a1), 0.5 * (b0 + b1), 0.5 * (c0 + c1), R.Bl, R.Bh, R.s, R.wl, R.wh) + I2 = second_order(R, Bs, ws, taylor_moments3(A, Bm, C, Bs, R.s * ws), with_H) + I = I1 if I2 is None else _intersect(I1, I2, 'H' if with_H else 'V') + return (I, I1, I2) if parts else I + + +def H2_point(B, T): + """H_2(B, T) = 2|B|^3 + 2T + V(B, T; middle root) on {B < 0, T < |B|^3/2}, 0 outside; near the boundary + 0 <= H_2(B, T) <= H_2(B, T') for T' <= T inside (H_2 is non-increasing in T).""" + if B >= 0.0: + return (0.0, 0.0) + half = half_cube(B) + if T >= half[1]: + return (0.0, 0.0) + if T > half[0]: + return (0.0, H2_point(B, dn(dn(half[0]) * (1.0 - 1e-12)))[1]) + v = add(add(scal(2.0, pw((-B, -B), 3)), (dn(2 * T), up(2 * T))), V_iv((B, B), xmid(B, T))) + return (max(0.0, v[0]), max(0.0, v[1])) + + +def _H2_center(Bs, ws, A, Bm, C, s): + w2 = (max(0.0, dn(ws * ws)), up(ws * ws)) + T2 = (dn(3.0 * w2[0]), up(3.0 * w2[1])) + if Bs < 0.0 and T2[1] < half_cube(Bs)[0]: + xs = xmid(Bs, T2[0], T2[1]) + Bsi = (Bs, Bs) + h0 = add(add(scal(-2.0, pw(Bsi, 3)), scal(2.0, T2)), V_iv(Bsi, xs)) + hB = add(scal(-6.0, sqr(Bsi)), VB_iv(Bsi, xs)) + vw = neg(Vw_big_iv(Bsi, xs)) # middle root: opposite sign + if MUTANT == 'mid-branch': + vw = neg(vw) + hw = add(scal(12.0, (ws, ws)), vw) + return h0, hB, hw + if Bs >= 0.0 or T2[0] > half_cube(Bs)[1]: + return (0.0, 0.0), (0.0, 0.0), (0.0, 0.0) # outside the support: value and gradient 0 + return None + + +def H2_box3(box, Ma, Mb, Mc, second=True, parts=False): + """Enclosure of int_box H_2 dmu. Inside the support the middle root x and r = B/x are non-increasing in B and T, + r in [-2, -1]; on boxes meeting the support boundary the zero extension of H_2 is C^1 with zero value and gradient + on {T = |B|^3/2} and piecewise C^2, so the Taylor remainder holds with Hessian ranges hulled with 0, the ranges being + taken over rectangle and support, where |B|/2 <= x <= x(Bl, Tl) and r in [-2, r(Bl, Tl)].""" + (a0, a1), (b0, b1), (c0, c1) = box + Bl, Bh, ul, uh = box_geometry(box) + zero = ((0.0, 0.0),) * 3 + if Bl >= 0.0: + return zero if parts else (0.0, 0.0) + s, wl, wh, Tl, Th = wT(ul, uh) + if Tl >= half_cube(Bl)[1]: + return zero if parts else (0.0, 0.0) + A, Bm, C = Ma.get(a0, a1), Mb.get(b0, b1), Mc.get(c0, c1) + m0 = mul(mul(A[0], Bm[0]), C[0]) + inside = Bh < 0.0 and Th < half_cube(Bh)[0] + I2 = None + if inside: + xLL, xHH = xmid(Bl, Tl), xmid(Bh, Th) + hmin = add(add(scal(2.0, pw((-Bh, -Bh), 3)), (dn(2 * Th), up(2 * Th))), V_iv((Bh, Bh), xHH))[0] + hmax = add(add(scal(2.0, pw((-Bl, -Bl), 3)), (dn(2 * Tl), up(2 * Tl))), V_iv((Bl, Bl), xLL))[1] + I1 = mul(m0, (max(0.0, hmin), max(0.0, hmax))) + X = (xHH[0], xLL[1]) + Rr = (max(-2.0, div((Bh, Bh), xHH)[0]), min(-1.0, div((Bl, Bl), xLL)[1])) + Bhe = Bh + else: + I1 = mul(m0, (H2_point(Bh, Th)[0], H2_point(Bl, Tl)[1])) + if Tl < half_cube(Bl)[0]: + Bhe = min(Bh, 0.0) + xLL = xmid(Bl, Tl) + X = (max(0.0, -Bhe / 2.0), xLL[1]) + Rr = (-2.0, min(-1.0, div((Bl, Bl), xLL)[1])) + else: + X = None + if second and s != 0 and X is not None: + SBB = add(scal(-12.0, (Bl, Bhe)), mul(X, p_range(Rr))) + SBw = mul(HSQRT3, mul(isqrt(X), e_range(Rr))) + Sww = add((12.0, 12.0), g_range(Rr)) + if not inside: + SBB, SBw, Sww = _hull0(SBB), _hull0(SBw), _hull0(Sww) + Bs, ws = _expansion_point(0.5 * (a0 + a1), 0.5 * (b0 + b1), 0.5 * (c0 + c1), Bl, Bh, s, wl, wh) + cen = _H2_center(Bs, ws, A, Bm, C, s) + if cen is not None: + h0, hB, hw = cen + m, E1B, E1u, QB, Qu, QBu = taylor_moments3(A, Bm, C, Bs, s * ws) + lin = add(add(mul(h0, m), mul(hB, E1B)), scal(float(s), mul(hw, E1u))) + I2 = add(lin, _quad(SBB, SBw, Sww, QB, Qu, QBu, s)) + I = I1 if I2 is None else _intersect(I1, I2, 'H2') + return (I, I1, I2) if parts else I + + +# ================================================================ 4D: the gap mark k as a box coordinate +K_LO, K_HI = 0.5, 2.0 +PMIN, PMAX = -6, 2 +KLEVEL = 12 # elementary dyadic cells of [1/2, 2]: width 3 * 2^-13 + + +def _cbrt_iv(x): + """Enclosure of x^(1/3) for a float x > 0.""" + c = x ** (1.0 / 3.0) + lo, hi = c, c + while pw((lo, lo), 3)[1] > x: + lo = dn(lo) + while pw((hi, hi), 3)[0] < x: + hi = up(hi) + return (lo, hi) + + +def kweights(k): + """f_p(k) = 3 k^(1/3) k^(p/2) e^(-12 k^2), p = PMIN..PMAX, at a float k > 0 (3k^2 x e^(-12k^2) k^(-5/3) x k^(p/2)).""" + third = _cbrt_iv(k) + if MUTANT == 'k-weight': + third = (1.0, 1.0) + base = mul(scal(3.0, third), exp_neg((dn(12.0 * dn(k * k)), up(12.0 * up(k * k))))) + sk = isqrt((k, k)) + isk = div((1.0, 1.0), sk) + out = {} + v = base + for p in range(0, PMAX + 1): + out[p] = v + v = mul(v, sk) + v = mul(base, isk) + for p in range(-1, PMIN - 1, -1): + out[p] = v + v = mul(v, isk) + return out + + +def elementary_kmoments(nsub=24): + """K_p = int f_p over each elementary cell. Every f_p is convex on [1/2, 2]: f''/f = (24k - q/k)^2 - q/k^2 - 24 with + q = 1/3 + p/2 in [-8/3, 4/3] is >= 87 - 16/3 - 24 > 0 (q >= 0) and >= 144 - 24 > 0 (q < 0). Hence on every sub-cell + h f(mid) <= int f <= h (f(l) + f(r))/2. All cell and sub-cell ends are exact dyadic floats.""" + n = 2 ** KLEVEL + width = (K_HI - K_LO) / n + cells = [] + prev = kweights(K_LO) + for i in range(n): + k0 = K_LO + width * i + acc = {p: (0.0, 0.0) for p in range(PMIN, PMAX + 1)} + for j in range(nsub): + l = k0 + width * j / nsub + r = k0 + width * (j + 1) / nsub + h = (dn(r - l), up(r - l)) + mid = kweights(0.5 * (l + r)) + fr = kweights(r) + for p in acc: + lo = mul(h, mid[p])[0] + hi = mul(h, scal(0.5, add(prev[p], fr[p])))[1] + acc[p] = (dn(acc[p][0] + lo), up(acc[p][1] + hi)) + prev = fr + cells.append(acc) + return cells + + +class KMoments: + """K_p on dyadic k-cells (unions of elementary cells), cached.""" + def __init__(self, elem): + self.elem, self.c = elem, {} + self.width = (K_HI - K_LO) / len(elem) + + def get(self, k0, k1): + key = (k0, k1) + r = self.c.get(key) + if r is None: + i0 = round((k0 - K_LO) / self.width) + i1 = round((k1 - K_LO) / self.width) + if K_LO + i0 * self.width != k0 or K_LO + i1 * self.width != k1: + raise ValueError('k-cell is not a union of elementary cells: %r' % (key,)) + r = {} + for p in range(PMIN, PMAX + 1): + lo = hi = 0.0 + for i in range(i0, i1): + e = self.elem[i][p] + lo, hi = dn(lo + e[0]), up(hi + e[1]) + r[p] = (lo, hi) + self.c[key] = r + return r + + +def _integrate4(terms, A, Bm, C, K): + acc = (0.0, 0.0) + for (cf, i, j, l, p) in terms: + acc = add(acc, mul(cf, mul(mul(A[i], Bm[j]), mul(C[l], K[p])))) + return acc + + +def _prod4(t1, t2): + out = {} + for (c1, i1, j1, l1, p1) in t1: + for (c2, i2, j2, l2, p2) in t2: + key = (i1 + i2, j1 + j2, l1 + l2, p1 + p2) + c = mul(c1, c2) + out[key] = add(out[key], c) if key in out else c + return [(c,) + key for key, c in sorted(out.items())] + + +_C12 = div((1.0, 1.0), (12.0, 12.0)) +_C72 = div((1.0, 1.0), (72.0, 72.0)) + + +def taylor_moments4(A, Bm, C, K, Bs, us): + """The six Taylor moments over a 4D box against W(k) dk x Gaussian, in the original coordinates: + dB = beta - a^2 k^-1/12 - Bs, du = c k^(1/2) - a beta k^(-1/2)/4 + a^3 k^(-3/2)/72 - us + (terms: coefficient, powers of a, beta, c, and p for k^(p/2)).""" + dB = [((1.0, 1.0), 0, 1, 0, 0), (neg(_C12), 2, 0, 0, -2), ((-Bs, -Bs), 0, 0, 0, 0)] + du = [((1.0, 1.0), 0, 0, 1, 1), ((-0.25, -0.25), 1, 1, 0, -1), (_C72, 3, 0, 0, -3), ((-us, -us), 0, 0, 0, 0)] + m = _integrate4([((1.0, 1.0), 0, 0, 0, 0)], A, Bm, C, K) + E1B = _integrate4(dB, A, Bm, C, K) + E1u = _integrate4(du, A, Bm, C, K) + QB = _integrate4(_prod4(dB, dB), A, Bm, C, K) + Qu = _integrate4(_prod4(du, du), A, Bm, C, K) + QBu = _integrate4(_prod4(dB, du), A, Bm, C, K) + return m, E1B, E1u, _nonneg(QB), _nonneg(Qu), QBu + + +def _scaled_range(x0, x1, f0, f1): + p = (x0 * f0, x0 * f1, x1 * f0, x1 * f1) + return dn(min(p)), up(max(p)) + + +def scaled_box(box): + """Outer box of the scaled coordinates (a/sqrt k, beta, sqrt(k) c) over a 4D box (a, beta, c, k).""" + (a0, a1), (b0, b1), (c0, c1), (k0, k1) = box + sk0, sk1 = isqrt((k0, k0)), isqrt((k1, k1)) + isk = (div((1.0, 1.0), sk1)[0], div((1.0, 1.0), sk0)[1]) + sk = (sk0[0], sk1[1]) + return (_scaled_range(a0, a1, isk[0], isk[1]), (b0, b1), _scaled_range(c0, c1, sk[0], sk[1])) + + +def H_box4(box, Ma, Mb, Mc, Km, second=True, parts=False): + """Enclosure of int_box W(k) H dmu dk over a 4D box (a, beta, c, k), W(k) = 3 k^(1/3) e^(-12k^2).""" + (a0, a1), (b0, b1), (c0, c1), (k0, k1) = box + A, Bm, C, K = Ma.get(a0, a1), Mb.get(b0, b1), Mc.get(c0, c1), Km.get(k0, k1) + R = Rect(*box_geometry(scaled_box(box))) + m0 = mul(mul(A[0], Bm[0]), mul(C[0], K[0])) + I1 = mul(m0, H_range(R)) + I2 = None + if second and R.s != 0: + ac, bc, cc, kc = 0.5 * (a0 + a1), 0.5 * (b0 + b1), 0.5 * (c0 + c1), 0.5 * (k0 + k1) + skc = math.sqrt(kc) + Bs, ws = _expansion_point(ac / skc, bc, cc * skc, R.Bl, R.Bh, R.s, R.wl, R.wh) + I2 = second_order(R, Bs, ws, taylor_moments4(A, Bm, C, K, Bs, R.s * ws), True) + I = I1 if I2 is None else _intersect(I1, I2, 'H (4D)') + return (I, I1, I2) if parts else I diff --git a/frontiers/c8_elder_failure_coefficient_20261001/certificate.py b/frontiers/c8_elder_failure_coefficient_20261001/certificate.py new file mode 100644 index 000000000..21318a7bb --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/certificate.py @@ -0,0 +1,730 @@ +#!/usr/bin/env python3 +"""Certificate for CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1 (standard library only). + + python3 -B -S certificate.py --write [--procs N] full run, writes RESULTS.json + python3 -B -S certificate.py --check [--procs N] self-test, then full replay; every replayed enclosure must lie + inside the published one (and the published one must not be + looser than the replay by more than the publication margin), + leaf counts must agree + python3 -B -S certificate.py --check --mutant NAME must exit 1 + python3 -B certificate.py --controls floating-point cross-checks (not part of the certificate) + +Certified quantities (planar near cluster coefficients of [NUM], reference kernel): + J_fail(k) = J_1 + J_2 = 3k^2 (E[2(B_-)^3] + E[2T] + E[V]), J_2(k) = 3k^2 E[H_2] (k = 1/2, 1, 2) + I_k = int_{1/2}^{2} e^(-12k^2) k^(-5/3) J_fail(k) dk, C_fail^{B,K} = const * I_b * I_k, + const = 4 / (sqrt2 (2 pi)^(5/2) sqrt12), I_b = int_0^1 e^(-b^2) db = sqrt(pi) (Phi(sqrt2) - 1/2). +d = 3 (conditional on Math-#175 (L1) and Math-#184 Theorem 1): a_j^(3) = R_3(b) alpha_j, R_3(0) = (32 + 28 sqrt2)/17, and + C_fail^(3) = const_3 I_b^(3) I_k, const_3 = 4 pi (4/sqrt2) (2 pi)^-4 12^(-1/2), I_b^(3) = (1/2) int_0^1 e^(-b^2) g(b) db, + g(b) = E[(b + sqrt2 Z)_+^3] = (b^3 + 6b) Phi(b/sqrt2) + sqrt2 (b^2 + 4) phi(b/sqrt2).""" +import sys, os, json, time, math, argparse +from decimal import Decimal, ROUND_FLOOR, ROUND_CEILING +from fractions import Fraction as Fr +from multiprocessing import Pool + +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +import ia +from ia import (dn, up, add, sub, neg, mul, div, scal, divc, sqr, pw, isqrt, exp_neg, Phi, Phi_point, SQRT2PI, + from_frac, PI) +import boxes as BX +import torus as TOR + +PACKET = 'CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1' +KS = (0.5, 1.0, 2.0) +TOL3 = {'V': 1e-8, 'H2': 1e-8} # per-box width tolerance (unscaled expectations) +TOL4 = 2.5e-9 +L_BOX = 8.0 # boxes cover [-8 sigma, 8 sigma] in every jet coordinate +N0, NK0 = 8, 8 # initial grid: 8 cells per jet coordinate, 8 dyadic k-cells +SK = 0.8 # k-scale used only to choose the split direction +MAXDEPTH3, MAXDEPTH4 = 60, 80 +MARGIN = 1e-7 # publication: widen relatively by MARGIN, then round outward to 10 digits +DIGITS = 10 +SQRT2 = isqrt((2.0, 2.0)) + + +# ================================================================ drivers +def sig3(k): + """Scaled-law standard deviations for k in KS (all variances exact floats).""" + v = ((2.0 / k), 2.0, 6.0 * k) + if v[0] * k != 2.0 or v[2] / k != 6.0: + raise ValueError('k must make 2/k and 6k exact') + return [(isqrt((x, x)), (x, x)) for x in v] + + +_ST = {} + + +def _init3(k): + S = sig3(k) + _ST['M'] = (BX.Moments(S[0][0], S[0][1], 6), BX.Moments(S[1][0], S[1][1], 2), BX.Moments(S[2][0], S[2][1], 2)) + _ST['sd'] = tuple(s[0][1] for s in S) + + +def _work3(args): + which, tol, cells, mutant = args + BX.MUTANT = mutant + Ma, Mb, Mc = _ST['M'] + sd = _ST['sd'] + fn = BX.H2_box3 if which == 'H2' else (lambda b, x, y, z: BX.big_box3(b, x, y, z, with_H=(which == 'H'))) + lo = hi = 0.0 + nleaf = 0 + for cell in cells: + stack = [(cell, 0)] + while stack: + box, dep = stack.pop() + I = fn(box, Ma, Mb, Mc) + if I[1] - I[0] > tol and dep < MAXDEPTH3: + w = [(box[i][1] - box[i][0]) / sd[i] for i in range(3)] + d = w.index(max(w)) + l_, h_ = box[d] + mid = 0.5 * (l_ + h_) + for half in ((l_, mid), (mid, h_)): + nb = list(box) + nb[d] = half + stack.append((tuple(nb), dep + 1)) + else: + lo, hi = dn(lo + I[0]), up(hi + I[1]) + nleaf += 1 + return lo, hi, nleaf + + +def run3(k, which, tol, procs, mutant=None): + """Sum of box enclosures of E[which] over [-L sigma, L sigma]^3 (which in V, H2, H).""" + S = sig3(k) + edges = [[-L_BOX * s[0][1] + 2 * L_BOX * s[0][1] * i / N0 for i in range(N0 + 1)] for s in S] + cells = [((edges[0][i], edges[0][i + 1]), (edges[1][j], edges[1][j + 1]), (edges[2][l], edges[2][l + 1])) + for i in range(N0) for j in range(N0) for l in range(N0)] + chunks = [(which, tol, cells[p::procs], mutant) for p in range(procs)] + with Pool(procs, initializer=_init3, initargs=(k,)) as pool: + res = pool.map(_work3, chunks) + lo = hi = 0.0 + n = 0 + for r in res: + lo, hi, n = dn(lo + r[0]), up(hi + r[1]), n + r[2] + return (lo, hi), n, edges + + +def _init4(elem): + S = sig3(1.0) # original jet law N(0, diag(2, 2, 6)) + _ST['M'] = (BX.Moments(S[0][0], S[0][1], 6), BX.Moments(S[1][0], S[1][1], 2), BX.Moments(S[2][0], S[2][1], 2)) + _ST['K'] = BX.KMoments(elem) + _ST['sd'] = tuple(s[0][1] for s in S) + (SK,) + + +def _work4(args): + tol, cells, mutant = args + BX.MUTANT = mutant + Ma, Mb, Mc = _ST['M'] + Km = _ST['K'] + sd = _ST['sd'] + lo = hi = 0.0 + nleaf = 0 + for cell in cells: + stack = [(cell, 0)] + while stack: + box, dep = stack.pop() + I = BX.H_box4(box, Ma, Mb, Mc, Km) + if I[1] - I[0] > tol and dep < MAXDEPTH4: + w = [(box[i][1] - box[i][0]) / sd[i] for i in range(4)] + if box[3][1] - box[3][0] <= Km.width: + w[3] = -1.0 + d = w.index(max(w)) + l_, h_ = box[d] + mid = 0.5 * (l_ + h_) + for half in ((l_, mid), (mid, h_)): + nb = list(box) + nb[d] = half + stack.append((tuple(nb), dep + 1)) + else: + lo, hi = dn(lo + I[0]), up(hi + I[1]) + nleaf += 1 + return lo, hi, nleaf + + +def run4(tol, procs, elem, mutant=None): + S = sig3(1.0) + edges = [[-L_BOX * s[0][1] + 2 * L_BOX * s[0][1] * i / N0 for i in range(N0 + 1)] for s in S] + kw = (BX.K_HI - BX.K_LO) / NK0 + kedges = [BX.K_LO + kw * j for j in range(NK0 + 1)] + cells = [((edges[0][i], edges[0][i + 1]), (edges[1][j], edges[1][j + 1]), (edges[2][l], edges[2][l + 1]), + (kedges[q], kedges[q + 1])) for q in range(NK0) for i in range(N0) for j in range(N0) for l in range(N0)] + chunks = [(tol, cells[p::procs], mutant) for p in range(procs)] + with Pool(procs, initializer=_init4, initargs=(elem,)) as pool: + res = pool.map(_work4, chunks) + lo = hi = 0.0 + n = 0 + for r in res: + lo, hi, n = dn(lo + r[0]), up(hi + r[1]), n + r[2] + return (lo, hi), n + + +# ================================================================ exact and one-dimensional parts, tails +def E2T(k): + """E[2T] = 6 E[u^2] = 6 (6k + 1/(4k) + 5/(216 k^3)) (exact rational).""" + k = Fr(k) + return 6 * (6 * k + Fr(1) / (4 * k) + Fr(5) / (216 * k ** 3)) + + +def _gs(mi): + """g = E[(m - beta)_+^3], g' = 3 E[(m - beta)_+^2], g'' = 6 E[(m - beta)_+] for beta ~ N(0, 2), t = m/sqrt2: + g = (m^3 + 6m) Phi(t) + sqrt2 (m^2 + 4) phi(t), g' = 3 [(m^2 + 2) Phi(t) + sqrt2 m phi(t)], + g'' = 6 [m Phi(t) + sqrt2 phi(t)]; all increasing in m.""" + t = div(mi, SQRT2) + P = Phi(t) + ph = div(exp_neg(divc(sqr(mi), 4.0)), SQRT2PI) + g = add(mul(add(pw(mi, 3), scal(6.0, mi)), P), mul(mul(SQRT2, add(sqr(mi), (4.0, 4.0))), ph)) + g1 = scal(3.0, add(mul(add(sqr(mi), (2.0, 2.0)), P), mul(mul(SQRT2, mi), ph))) + g2 = scal(6.0, add(mul(mi, P), mul(SQRT2, ph))) + return g, g1, g2 + + +def tail_moments(L, s, s2, nmax): + """int_L^inf x^n phi_s(x) dx, n = 0..nmax, L >= 0.""" + pL = BX.phis(L, s, s2) + q = BX._Qiv(div((L, L), s)) + M = [q, mul(s2, pL)] + for n in range(2, nmax + 1): + M.append(add(mul(scal(float(n - 1), s2), M[n - 2]), mul(s2, mul(pw((L, L), n - 1), pL)))) + return M + + +def E2B(k, L=12.0, N=4000): + """E[2(B_-)^3] = 4 int_0^inf phi_sa(a) G(a) da with G(a) = g(a^2/12), sa^2 = 2/k: second-order Taylor expansion on each + of N cells about the midpoint (G'' = g''(m) a^2/36 + g'(m)/6 is increasing in a >= 0, bounded by its values at the + cell ends), plus the tail a > L sa with g(m) <= m^3 + 6m + sqrt2 (m^2 + 4)/sqrt(2 pi).""" + s, s2 = sig3(k)[0] + top = L * s[1] + acc = (0.0, 0.0) + for i in range(N): + a0, a1 = top * i / N, top * (i + 1) / N + M = BX.moments(a0, a1, s, s2, 2) + ac = 0.5 * (a0 + a1) + aci = (ac, ac) + g, g1, _ = _gs(divc(sqr(aci), 12.0)) + _, g1l, g2l = _gs(divc(sqr((a0, a0)), 12.0)) + _, g1h, g2h = _gs(divc(sqr((a1, a1)), 12.0)) + G2lo = add(divc(mul(g2l, sqr((a0, a0))), 36.0), divc(g1l, 6.0))[0] + G2hi = add(divc(mul(g2h, sqr((a1, a1))), 36.0), divc(g1h, 6.0))[1] + M1c = sub(M[1], mul(aci, M[0])) + M2c = BX._nonneg(add(sub(M[2], scal(2.0, mul(aci, M[1]))), mul(sqr(aci), M[0]))) + cell = add(add(mul(g, M[0]), mul(mul(g1, divc(aci, 6.0)), M1c)), scal(0.5, mul((G2lo, G2hi), M2c))) + acc = add(acc, cell) + T = tail_moments(top, s, s2, 6) + tail = add(add(divc(T[6], 1728.0), divc(T[2], 2.0)), + mul(div(SQRT2, SQRT2PI), add(divc(T[4], 144.0), scal(4.0, T[0])))) + return (scal(4.0, acc)[0], scal(4.0, add(acc, (0.0, tail[1])))[1]) + + +def abs_moments(s, s2, nmax): + """E|X|^n, X ~ N(0, s^2): 1, 2 s^2 phi_s(0), then (n - 1) s^2 E|X|^(n-2).""" + M = [(1.0, 1.0), scal(2.0, mul(s2, BX.phis(0.0, s, s2)))] + for n in range(2, nmax + 1): + M.append(mul(scal(float(n - 1), s2), M[n - 2])) + return M + + +def _tail_union(S, L, terms): + """Upper bound of E[sum_terms c |X1|^i |X2|^j |X3|^l 1{some |X_m| > L sigma_m}] (union bound over m; independence).""" + full = [abs_moments(S[m][0], S[m][1], 6) for m in range(3)] + tl = [[scal(2.0, t) for t in tail_moments(L * S[m][0][1], S[m][0], S[m][1], 6)] for m in range(3)] + tot = (0.0, 0.0) + for m in range(3): + for cf, e in terms: + v = cf + for j in range(3): + v = mul(v, tl[j][e[j]] if j == m else full[j][e[j]]) + tot = add(tot, v) + return tot[1] + + +def tail3(k): + """Outside the box, |V|, H and H_2 are at most 17|B|^3 + 8T <= 68|beta|^3 + (68/1728 + 96/5184) a^6 + 48 c^2 + + 6 a^2 beta^2 (scaled coordinates; x <= |B| + (T/2)^(1/3), |B|^3 <= 4(|beta|^3 + a^6/1728), T <= 6c^2 + 6d^2, + d^2 <= 2(a^2 beta^2/16 + a^6/5184)).""" + c6 = add(divc((68.0, 68.0), 1728.0), divc((96.0, 96.0), 5184.0)) + terms = [((68.0, 68.0), (0, 3, 0)), (c6, (6, 0, 0)), ((48.0, 48.0), (0, 0, 2)), ((6.0, 6.0), (2, 2, 0))] + return _tail_union(sig3(k), L_BOX, terms) + + +def tail4(Ktot): + """Outside the jet box, for every k in [1/2, 2] (original coordinates): H <= 17|B|^3 + 8T <= 68|beta|^3 + 68 a^6/216 + + 576 a^6/5184 + 144 c^2 + 9 a^2 beta^2 (|B| <= |beta| + a^2/6; u^2 <= 6c^2 + 3a^2 beta^2/8 + 24 a^6/5184); times + int W(k) dk = K_0([1/2, 2]).""" + c6 = add(divc((68.0, 68.0), 216.0), divc((576.0, 576.0), 5184.0)) + terms = [((68.0, 68.0), (0, 3, 0)), (c6, (6, 0, 0)), ((144.0, 144.0), (0, 0, 2)), ((9.0, 9.0), (2, 2, 0))] + return up(_tail_union(sig3(1.0), L_BOX, terms) * Ktot[1] * (1.0 + 4e-16)) + + +# ================================================================ assembly +def m2(b): + """m_(2,b) = (b^2 + 2) Phi(b/sqrt2) + sqrt2 b phi(b/sqrt2) (b a float >= 0).""" + bi = (b, b) + t = div(bi, SQRT2) + ph = div(exp_neg(divc(sqr(t), 2.0)), SQRT2PI) + return add(mul(add(sqr(bi), (2.0, 2.0)), Phi(t)), mul(mul(SQRT2, bi), ph)) + + +def prefactor(b, k): + """p_b(0)/z_0 = phi(b/sqrt2) / (sqrt2 36 k^2 m_(2,b)).""" + t = div((b, b), SQRT2) + ph = div(exp_neg(divc(sqr(t), 2.0)), SQRT2PI) + return div(ph, mul(mul(SQRT2, (36.0 * k * k, 36.0 * k * k)), m2(b))) + + +def constants(): + sqrt_pi = isqrt(PI) + two_pi = scal(2.0, PI) + const = div((4.0, 4.0), mul(mul(SQRT2, mul(sqr(two_pi), isqrt(two_pi))), isqrt((12.0, 12.0)))) + Ib = mul(sqrt_pi, sub(Phi(SQRT2), (0.5, 0.5))) + return const, Ib + + +def constants3(): + """d = 3: const_3 = 4 pi (4/sqrt2) (2 pi)^-4 12^(-1/2) and R_3(0) = (32 + 28 sqrt2)/17; also checks + N_3(0) = pi (2 pi)^(-1/2) g(0) = 2 sqrt2, the identity behind I_b^(3).""" + two_pi = scal(2.0, PI) + const3 = mul(mul(scal(4.0, PI), div((4.0, 4.0), SQRT2)), div((1.0, 1.0), mul(pw(two_pi, 4), isqrt((12.0, 12.0))))) + R30 = divc(add((32.0, 32.0), scal(28.0, SQRT2)), 17.0) + N30 = mul(div(PI, isqrt(two_pi)), _gs((0.0, 0.0))[0]) + two_sqrt2 = scal(2.0, SQRT2) + if N30[1] < two_sqrt2[0] or two_sqrt2[1] < N30[0]: + raise RuntimeError('N_3(0) = 2 sqrt2 check failed: %r' % (N30,)) + return const3, R30 + + +def Ib3(N=100000): + """I_b^(3) = (1/2) int_0^1 e^(-b^2) g(b) db: on each cell e^(-b^2) is decreasing and g increasing (both positive). + (The prefactor e^(-3b^2/4) phi(b/sqrt2) pi (2 pi)^(-1/2) of N_3 equals e^(-b^2)/2.)""" + lo = hi = 0.0 + prev_g = _gs((0.0, 0.0))[0] + prev_e = (1.0, 1.0) + for i in range(1, N + 1): + b0, b1 = (i - 1) / N, i / N + h = (dn(b1 - b0), up(b1 - b0)) + g1 = _gs((b1, b1))[0] + e1 = exp_neg(sqr((b1, b1))) + lo = dn(lo + mul(h, mul(e1, prev_g))[0]) + hi = up(hi + mul(h, mul(prev_e, g1))[1]) + prev_g, prev_e = g1, e1 + return (dn(0.5 * lo), up(0.5 * hi)) + + +def iv_frac(q): + return from_frac(q) + + +def certify(procs, mutant=None, log=print, only=None): + """Run everything; returns the raw results (floats) and leaf counts.""" + out = {'J': {}, 'leaves': {}} + t0 = time.time() + for k in KS: + if only is not None and k not in only: + continue + s = 3.0 * k * k + EV, nV, _ = run3(k, 'V', TOL3['V'], procs, mutant) + EH2, nH2, _ = run3(k, 'H2', TOL3['H2'], procs, mutant) + e2b = E2B(k) + e2t = iv_frac(E2T(k)) + tl = tail3(k) + EVt = (dn(EV[0] - tl), up(EV[1] + tl)) + EH2t = (EH2[0], up(EH2[1] + tl)) + Jf = scal(s, add(add(e2b, e2t), EVt)) + J2 = scal(s, EH2t) + out['J'][k] = {'J_fail': Jf, 'J_2': J2, 'E2B': e2b, 'E2T': e2t, 'EV_boxes': EV, 'EH2_boxes': EH2, 'tail': tl} + out['leaves']['V k=%g' % k] = nV + out['leaves']['H2 k=%g' % k] = nH2 + log('k=%g J_fail in [%.10f, %.10f] J_2 in [%.10f, %.10f] (leaves %d, %d; %.0fs)' % + (k, Jf[0], Jf[1], J2[0], J2[1], nV, nH2, time.time() - t0)) + if only is None: + elem = BX.elementary_kmoments() + Km = BX.KMoments(elem) + Ktot = Km.get(BX.K_LO, BX.K_HI)[0] + Ik_boxes, n4 = run4(TOL4, procs, elem, mutant) + t4 = tail4(Ktot) + Ik = (Ik_boxes[0], up(Ik_boxes[1] + t4)) + const, Ib = constants() + C = mul(mul(const, Ib), Ik) + out['Ik'] = {'Ik': Ik, 'boxes': Ik_boxes, 'tail': t4, 'const': const, 'Ib': Ib, 'C_fail': C, 'Wtot': Ktot, + 'K_m4': Km.get(BX.K_LO, BX.K_HI)[-4]} + out['torus'] = TOR.transfer(TOR.E_I2_upper()) + const3, R30 = constants3() + Ib3_ = Ib3() + out['d3'] = {'const3': const3, 'R30': R30, 'Ib3': Ib3_, 'C_fail3': mul(mul(const3, Ib3_), Ik)} + out['leaves']['I_k 4D'] = n4 + log('I_k in [%.10f, %.10f] C_fail in [%.12f, %.12f] (leaves %d; %.0fs)' % (Ik[0], Ik[1], C[0], C[1], n4, + time.time() - t0)) + return out + + +# ================================================================ publication rounding and the result document +def pub(iv_): + """Outward publication of an enclosure: widen by MARGIN relatively, round outward to DIGITS significant digits.""" + lo, hi = Decimal(iv_[0]), Decimal(iv_[1]) + lo = lo - abs(lo) * Decimal(MARGIN) + hi = hi + abs(hi) * Decimal(MARGIN) + + def rnd(x, mode): + if x == 0: + return Decimal(0) + e = x.adjusted() - (DIGITS - 1) + return x.quantize(Decimal(1).scaleb(e), rounding=mode) + return [str(rnd(lo, ROUND_FLOOR)), str(rnd(hi, ROUND_CEILING))] + + +def contains_pub(pub_iv, iv_): + return Decimal(pub_iv[0]) <= Decimal(iv_[0]) and Decimal(iv_[1]) <= Decimal(pub_iv[1]) + + +def not_looser(pub_iv, iv_): + """The published interval is within 2 MARGIN + one unit of the last digit of the replayed one.""" + lo, hi = Decimal(iv_[0]), Decimal(iv_[1]) + unit = lambda x: Decimal(1).scaleb(x.adjusted() - (DIGITS - 1)) if x != 0 else Decimal(0) + lo_min = lo - abs(lo) * Decimal(2 * MARGIN) - unit(lo) + hi_max = hi + abs(hi) * Decimal(2 * MARGIN) + unit(hi) + return Decimal(pub_iv[0]) >= lo_min and Decimal(pub_iv[1]) <= hi_max + + +def derived(raw): + """Published quantities derived from the raw enclosures.""" + res = {} + for k in KS: + J = raw['J'][k] + Jf, J2 = J['J_fail'], J['J_2'] + J1 = sub(Jf, J2) + share = div(J2, Jf) + palm = div(scal(2.0, J2), add(Jf, J2)) + row = {'J_fail': Jf, 'J_1': J1, 'J_2': J2, 'share_J2_over_Jfail': share, 'near_palm_excess': palm, + 'E[2T] exact': str(E2T(k)), + 'parts': {'E[2(B_-)^3]': J['E2B'], 'E[2T]': J['E2T'], 'E[V] boxes': J['EV_boxes'], + 'E[H_2] boxes': J['EH2_boxes'], 'tail bound': (0.0, J['tail'])}} + coeffs = {} + for b in (0.0, 1.0): + P = prefactor(b, k) + coeffs['b=%g' % b] = {'p_b(0)/z_0': P, 'alpha_1': mul(P, J1), 'alpha_2': mul(P, J2), + 'a_fail = alpha_1 + alpha_2': mul(P, Jf)} + row['coefficients'] = coeffs + res['k=%g' % k] = row + if 'd3' in raw: + D = raw['d3'] + sec = {'R_3(0) = (32 + 28 sqrt2)/17': D['R30'], 'const_3': D['const3'], 'I_b^(3)': D['Ib3'], + 'C_fail^(3),{B,K}': D['C_fail3']} + for k in KS: + c = res['k=%g' % k]['coefficients']['b=0'] + sec['k=%g, b=0' % k] = {'alpha_1^(3)': mul(D['R30'], c['alpha_1']), 'alpha_2^(3) = nu^(3)(2)': mul(D['R30'], c['alpha_2']), + 'a_fail^(3)': mul(D['R30'], c['a_fail = alpha_1 + alpha_2'])} + res['d=3 (conditional on Math-#175 (L1) and Math-#184 Theorem 1)'] = sec + if 'torus' in raw and 'Ik' in raw: + T = raw['torus'] + dJ = (-T['dJ'], T['dJ']) + sec = {'jet covariance deviation delta (order <= 6)': (0.0, T['delta']), 'chi^2 bound': (0.0, T['chi2']), + '|J^L - J| bound (all k, J_fail and J_2)': (0.0, T['dJ']), + '(pi_0 p_b(0))^L / (pi_0 p_b(0))^ref': T['pi0p_ratio']} + for k in KS: + J = raw['J'][k] + JfL, J2L = add(J['J_fail'], dJ), add(J['J_2'], dJ) + row = {'J_fail': JfL, 'J_2': J2L} + for b in (0.0, 1.0): + mu = mul(T['mean_coeff_A'], (b, b)) + s = isqrt(T['var_A']) + t = div(mu, s) + ph = div(exp_neg(divc(sqr(t), 2.0)), SQRT2PI) + p = div(ph, s) + m2L = sub(mul(add(sqr(mu), T['var_A']), sub((1.0, 1.0), Phi(t))), mul(mul(mu, s), ph)) + P = div(p, mul((36.0 * k * k, 36.0 * k * k), m2L)) + row['b=%g' % b] = {'p_b(0)/z_0': P, 'alpha_1': mul(P, sub(JfL, J2L)), 'alpha_2': mul(P, J2L), + 'a_fail': mul(P, JfL)} + sec['k=%g' % k] = row + I = raw['Ik'] + dC = mul(mul(mul(I['const'], I['Ib']), divc(I['K_m4'], 3.0)), (T['dJ'], T['dJ'])) + CL = mul(T['pi0p_ratio'], add(I['C_fail'], (-dC[1], dC[1]))) + sec['C_fail^{B,K}'] = CL + res['torus model, every L >= 10, every frame (Lemma T)'] = sec + if 'Ik' in raw: + I = raw['Ik'] + res['C_fail'] = {'I_k': I['Ik'], 'I_k boxes': I['boxes'], 'I_k tail bound': (0.0, I['tail']), + 'int_K W(k) dk': I['Wtot'], 'const': I['const'], 'I_b': I['Ib'], 'C_fail^{B,K}': I['C_fail']} + return res + + +def _walk(d, f): + if isinstance(d, dict): + return {k: _walk(v, f) for k, v in d.items()} + if isinstance(d, str): + return d + return f(d) + + +def document(raw): + der = derived(raw) + return { + 'object': PACKET, + 'scientific_effect': 'NONE', + 'certified': True, + 'arithmetic': 'binary64 interval arithmetic, outward one-ulp widening after every operation; published ' + 'intervals widened by %g relatively and rounded outward to %d significant digits' % (MARGIN, DIGITS), + 'parameters': {'k': list(KS), 'tol3': TOL3, 'tol4': TOL4, 'box_half_width_sigma': L_BOX, 'initial_cells': N0, + 'initial_k_cells': NK0, 'k_split_scale': SK, 'k_elementary_level': BX.KLEVEL}, + 'leaves': raw['leaves'], + 'results': _walk(der, pub), + } + + +def compare(doc, raw): + """Replay check: every replayed interval inside the published one, not looser by more than the margin.""" + der = derived(raw) + bad = [] + + def rec(p, r, path): + if isinstance(r, dict): + if sorted(p) != sorted(r): + bad.append(('/'.join(path), 'keys', sorted(p), sorted(r))) + return + for key in r: + rec(p[key], r[key], path + [key]) + elif isinstance(r, str): + if p != r: + bad.append(('/'.join(path), p, r)) + else: + if not contains_pub(p, r) or not not_looser(p, r): + bad.append(('/'.join(path), p, r)) + rec(doc['results'], der, []) + if doc['leaves'] != raw['leaves']: + bad.append(('leaves', doc['leaves'], raw['leaves'])) + return bad + + +# ================================================================ self-test (pure-Python Gauss-Legendre references) +def _gl(n): + xs, ws = [], [] + for i in range(1, n + 1): + x = math.cos(math.pi * (i - 0.25) / (n + 0.5)) + for _ in range(100): + p0, p1 = 1.0, x + for j in range(2, n + 1): + p0, p1 = p1, ((2 * j - 1) * x * p1 - (j - 1) * p0) / j + dp = n * (x * p1 - p0) / (x * x - 1) + dx = p1 / dp + x -= dx + if abs(dx) < 1e-16: + break + xs.append(x) + ws.append(2.0 / ((1 - x * x) * dp * dp)) + return xs, ws + + +def _xb_f(B, T): + xm = -B if B < 0 else B / 2.0 + return xm if T <= 0 else BX._newton_big(B, T, xm) + + +def _xm_f(B, T): + b = BX.xmid(B, T) + return 0.5 * (b[0] + b[1]) + + +def _H_f(B, T): + x = _xb_f(B, T) + return abs(B) ** 3 + 4 * x ** 3 - 3 * B * B * x + + +def _V_f(B, T): + x = _xb_f(B, T) + return -3 * B * (2 * x - B) * (x + B) + + +def _H2_f(B, T): + if B >= 0 or T >= abs(B) ** 3 / 2: + return 0.0 + x = _xm_f(B, T) + return abs(B) ** 3 + 4 * x ** 3 - 3 * B * B * x + + +def _ref3(box, k, f, n=12): + xs, ws = _gl(n) + var = (2.0 / k, 2.0, 6.0 * k) + pts = [] + for (l, h), v in zip(box, var): + pts.append([(0.5 * (h - l) * x + 0.5 * (h + l), + 0.5 * (h - l) * w * math.exp(-(0.5 * (h - l) * x + 0.5 * (h + l)) ** 2 / (2 * v)) / math.sqrt(2 * math.pi * v)) + for x, w in zip(xs, ws)]) + tot = 0.0 + for a, wa in pts[0]: + for be, wb in pts[1]: + B = be - a * a / 12.0 + d = -a * be / 4.0 + a ** 3 / 72.0 + for c, wc in pts[2]: + u = c + d + tot += wa * wb * wc * f(B, 3 * u * u) + return tot + + +def _ref4(box, n=10): + xs, ws = _gl(n) + var = (2.0, 2.0, 6.0) + pts = [] + for (l, h), v in zip(box[:3], var): + pts.append([(0.5 * (h - l) * x + 0.5 * (h + l), + 0.5 * (h - l) * w * math.exp(-(0.5 * (h - l) * x + 0.5 * (h + l)) ** 2 / (2 * v)) / math.sqrt(2 * math.pi * v)) + for x, w in zip(xs, ws)]) + l, h = box[3] + kp = [(0.5 * (h - l) * x + 0.5 * (h + l), 0.5 * (h - l) * w) for x, w in zip(xs, ws)] + tot = 0.0 + for k, wk in kp: + Wk = wk * 3 * k ** (1.0 / 3.0) * math.exp(-12 * k * k) + sk = math.sqrt(k) + for a, wa in pts[0]: + for be, wb in pts[1]: + B = be - a * a / (12 * k) + for c, wc in pts[2]: + u = sk * c - a * be / (4 * sk) + a ** 3 / (72 * k * sk) + tot += Wk * wa * wb * wc * _H_f(B, 3 * u * u) + return tot + + +SELFTEST3 = [ # (k, integrand, box in scaled coordinates (a, beta, c)); boxes on one side of u = 0 + (0.5, 'V', ((0.30, 0.55), (-1.20, -0.95), (0.80, 1.10))), + (0.5, 'V', ((-2.10, -1.80), (0.40, 0.70), (-1.60, -1.30))), + (0.5, 'V', ((0.05, 0.25), (0.90, 1.20), (2.00, 2.30))), + (0.5, 'H', ((1.00, 1.30), (-0.60, -0.30), (0.30, 0.60))), + (1.0, 'V', ((0.10, 0.35), (-2.00, -1.70), (0.20, 0.45))), + (1.0, 'H', ((-0.80, -0.55), (0.10, 0.40), (-2.60, -2.30))), + (2.0, 'V', ((0.40, 0.60), (-0.90, -0.60), (3.00, 3.40))), + (0.5, 'H2', ((0.20, 0.40), (-2.40, -2.15), (0.60, 0.80))), # inside the support + (1.0, 'H2', ((-0.30, -0.10), (-3.00, -2.70), (-1.20, -0.95))), # inside the support + (0.5, 'H2', ((0.10, 0.30), (-1.30, -1.10), (0.60, 0.85))), # meets the support boundary + (0.5, 'V', ((0.30, 0.32), (-1.00, -0.98), (0.50, 1.50))), # elongated in c: the w-curvature dominates + (1.0, 'H', ((0.30, 0.32), (0.40, 0.42), (1.00, 2.00))), + (0.5, 'H2', ((0.20, 0.22), (-2.40, -2.38), (0.20, 0.90))), +] +SELFTEST4 = [ + ((0.40, 0.70), (-1.20, -0.90), (1.00, 1.40), (0.5, 0.5 + 3 * 2.0 ** -7)), + ((-1.50, -1.20), (0.50, 0.80), (-2.00, -1.60), (0.5 + 3 * 2.0 ** -6, 0.5 + 3 * 2.0 ** -5)), +] + + +def selftest(mutant=None, log=print): + """Every box enclosure (and its first- and second-order parts separately) must contain a high-order floating-point + Gauss-Legendre value; the k-moments must contain Gauss-Legendre values. Returns a list of failures.""" + BX.MUTANT = mutant + fails = [] + for k, which, box in SELFTEST3: + S = sig3(k) + Ma, Mb, Mc = BX.Moments(S[0][0], S[0][1], 6), BX.Moments(S[1][0], S[1][1], 2), BX.Moments(S[2][0], S[2][1], 2) + if which == 'H2': + I, I1, I2 = BX.H2_box3(box, Ma, Mb, Mc, parts=True) + ref = _ref3(box, k, _H2_f, 16) + slack = 1e-6 * abs(ref) + 1e-15 # the zero extension is only C^1 at the support boundary + else: + try: + I, I1, I2 = BX.big_box3(box, Ma, Mb, Mc, with_H=(which == 'H'), parts=True) + except RuntimeError as e: + fails.append((which, k, box, 'raised: %s' % e)) + continue + ref = _ref3(box, k, _H_f if which == 'H' else _V_f) + slack = 1e-9 * abs(ref) + 1e-15 + for name, J in (('I', I), ('I1', I1), ('I2', I2)): + if J is not None and not (J[0] - slack <= ref <= J[1] + slack): + fails.append((which, k, box, name, J, ref)) + if I2 is None: + fails.append((which, k, box, 'no second-order enclosure')) + S = sig3(1.0) + Ma, Mb, Mc = BX.Moments(S[0][0], S[0][1], 6), BX.Moments(S[1][0], S[1][1], 2), BX.Moments(S[2][0], S[2][1], 2) + elem = {} + for box in SELFTEST4: + k0, k1 = box[3] + w = (BX.K_HI - BX.K_LO) / 2 ** BX.KLEVEL + i0, i1 = round((k0 - BX.K_LO) / w), round((k1 - BX.K_LO) / w) + Km = _partial_kmoments(i0, i1) # elementary k-moments of the needed cells only + try: + I, I1, I2 = BX.H_box4(box, Ma, Mb, Mc, Km, parts=True) + except RuntimeError as e: + fails.append(('H4', box, 'raised: %s' % e)) + continue + ref = _ref4(box) + slack = 1e-9 * abs(ref) + 1e-15 + for name, J in (('I', I), ('I1', I1), ('I2', I2)): + if J is not None and not (J[0] - slack <= ref <= J[1] + slack): + fails.append(('H4', box, name, J, ref)) + if I2 is None: + fails.append(('H4', box, 'no second-order enclosure')) + xs, ws = _gl(20) + K = Km.get(k0, k1) + for p in (BX.PMIN, 0, BX.PMAX): + refk = sum(0.5 * (k1 - k0) * wq * 3 * kk ** (1 / 3 + p / 2) * math.exp(-12 * kk * kk) + for kk, wq in ((0.5 * (k1 - k0) * x + 0.5 * (k1 + k0), wq) for x, wq in zip(xs, ws))) + if not (K[p][0] * (1 - 1e-12) <= refk <= K[p][1] * (1 + 1e-12)): + fails.append(('K', box[3], p, K[p], refk)) + log('self-test: %d boxes, %d failures' % (len(SELFTEST3) + len(SELFTEST4), len(fails))) + return fails + + +def _partial_kmoments(i0, i1, nsub=24): + """KMoments object whose elementary cells i0..i1-1 are computed (others are never touched).""" + n = 2 ** BX.KLEVEL + width = (BX.K_HI - BX.K_LO) / n + elem = [None] * n + for i in range(i0, i1): + k0 = BX.K_LO + width * i + acc = {p: (0.0, 0.0) for p in range(BX.PMIN, BX.PMAX + 1)} + prev = BX.kweights(k0) + for j in range(nsub): + l = k0 + width * j / nsub + r = k0 + width * (j + 1) / nsub + h = (dn(r - l), up(r - l)) + mid = BX.kweights(0.5 * (l + r)) + fr = BX.kweights(r) + for p in acc: + lo = mul(h, mid[p])[0] + hi = mul(h, scal(0.5, add(prev[p], fr[p])))[1] + acc[p] = (dn(acc[p][0] + lo), up(acc[p][1] + hi)) + prev = fr + elem[i] = acc + return BX.KMoments(elem) + + +# ================================================================ main +def main(): + ap = argparse.ArgumentParser() + ap.add_argument('--write', action='store_true') + ap.add_argument('--check', action='store_true') + ap.add_argument('--controls', action='store_true') + ap.add_argument('--mutant', default=None) + ap.add_argument('--procs', type=int, default=os.cpu_count() or 2) + a = ap.parse_args() + path = os.path.join(HERE, 'RESULTS.json') + if a.controls: + import controls + controls.main() + return 0 + fails = selftest(a.mutant) + if fails: + for f in fails[:10]: + print('SELF-TEST FAILURE', f) + return 1 + try: + raw = certify(a.procs, a.mutant) + except RuntimeError as e: + print('CERTIFICATE FAILURE', e) + return 1 + if a.write: + doc = document(raw) + with open(path, 'w') as fh: + json.dump(doc, fh, indent=1, sort_keys=False) + fh.write('\n') + print('wrote', path) + return 0 + if a.check: + with open(path) as fh: + doc = json.load(fh) + bad = compare(doc, raw) + if bad: + for b in bad[:20]: + print('MISMATCH', b) + return 1 + print(json.dumps({'object': PACKET, 'check': 'passed', 'scientific_effect': 'NONE'})) + return 0 + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/frontiers/c8_elder_failure_coefficient_20261001/controls.py b/frontiers/c8_elder_failure_coefficient_20261001/controls.py new file mode 100644 index 000000000..396138322 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/controls.py @@ -0,0 +1,106 @@ +"""Floating-point controls for CL-C8-ELDER-FAILURE-COEFFICIENT-20261001-v1 (not part of the certificate). + +1. E[2(B_-)^3] by composite Simpson (float) against the certified enclosure. +2. Consistency of the split J_fail = 3k^2 (E[2(B_-)^3] + E[2T] + E[V]) with a direct box enclosure of E[H] at k = 1/2 + (same engine, different integrand: the two enclosures must intersect). +3. The certified J_1, J_2 against the Gauss-Hermite and Monte Carlo values of [NUM] (RESULTS.json on this tree). +4. C_fail^{B,K} against the cap-route constants of Math-#215 (C_{B,K}(1/4096) = 402.5297, Gamma_2 = 307.51), and the + d = 3 constant against C_{B,K}^(3)(1/4096) = 1205.7089, Gamma_3 = 540.97. +5. The d = 3 near coefficients at b = 0 against Math-#184's floating values. +6. Lemma T: floating lattice sums of the torus jet-covariance deviation at L = 10 against the certified delta.""" +import os, sys, json, math +from decimal import Decimal +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +import certificate as C +from ia import dn, up + +NUM = os.path.join(HERE, '..', 'c6_cluster_coefficients_numerics_20260930', 'RESULTS.json') + + +def e2b_float(k, N=200000): + sa = math.sqrt(2 / k) + def g(m): + return (m ** 3 + 6 * m) * 0.5 * (1 + math.erf(m / 2)) + math.sqrt(2) * (m * m + 4) * math.exp(-m * m / 4) / math.sqrt(2 * math.pi) + L = 12 * sa + h = L / N + f = lambda a: math.exp(-a * a / (2 * sa * sa)) / (sa * math.sqrt(2 * math.pi)) * g(a * a / 12) + S = f(0) + f(L) + sum((4 if i % 2 else 2) * f(i * h) for i in range(1, N)) + return 4 * S * h / 3 + + +def main(): + doc = json.load(open(os.path.join(HERE, 'RESULTS.json'))) + R = doc['results'] + print('1. E[2(B_-)^3]: float Simpson vs certified') + for k in C.KS: + e = e2b_float(k) + iv = R['k=%g' % k]['parts']['E[2(B_-)^3]'] + print(' k=%g float %.12f certified [%s, %s] inside: %s' % (k, e, iv[0], iv[1], Decimal(iv[0]) <= Decimal(e) <= Decimal(iv[1]))) + print('2. direct E[H] boxes at k = 1/2 (tol 1e-7) vs the split') + (lo, hi), n, _ = C.run3(0.5, 'H', 1e-7, os.cpu_count() or 2) + tl = C.tail3(0.5) + s = 0.75 + J = R['k=0.5']['J_fail'] + direct = (s * lo, s * (hi + tl)) + print(' direct 3k^2 E[H] in [%.9f, %.9f] (%d leaves); split J_fail [%s, %s]; intersect: %s' % + (direct[0], direct[1], n, J[0], J[1], Decimal(direct[0]) <= Decimal(J[1]) and Decimal(J[0]) <= Decimal(direct[1]))) + print('3. [NUM] floating values against the certified enclosures') + num = json.load(open(NUM))['per_k'] + for k, key in ((0.5, '1/2'), (1.0, '1'), (2.0, '2')): + r = R['k=%g' % k] + for j in ('1', '2'): + iv = [float(x) for x in r['J_%s' % j]] + mid = 0.5 * (iv[0] + iv[1]) + v = num[key] + line = ' k=%g J_%s certified [%.6f, %.6f]' % (k, j, iv[0], iv[1]) + for lab in ('gh40', 'gh60', 'mc'): + x = v['J%s_%s' % (j, lab)] + inside = iv[0] <= x <= iv[1] + line += ' %s %.5f (%+.3f%%%s)' % (lab, x, 100 * (x - mid) / mid, ', inside' if inside else '') + line += ' mc s.e. %.4f -> (mc - mid)/se = %+.2f' % (v['J%s_mc_se' % j], (v['J%s_mc' % j] - mid) / v['J%s_mc_se' % j]) + print(line) + print('4. C_fail^{B,K} against the cap-route constants') + Cf = [float(x) for x in R['C_fail']['C_fail^{B,K}']] + for name, val in (('C_{B,K}(1/4096), Math-#215', 402.5297), ('Gamma_2 (floating), Math-#215', 307.51)): + print(' %s = %g : ratio to C_fail in [%.4g, %.4g]' % (name, val, val / Cf[1], val / Cf[0])) + D3 = R['d=3 (conditional on Math-#175 (L1) and Math-#184 Theorem 1)'] + Cf3 = [float(x) for x in D3['C_fail^(3),{B,K}']] + for name, val in (('C_{B,K}^(3)(1/4096), Math-#215', 1205.7089), ('Gamma_3 (floating), Math-#215', 540.97)): + print(' d = 3: %s = %g : ratio to C_fail^(3) in [%.4g, %.4g]' % (name, val, val / Cf3[1], val / Cf3[0])) + print('5. d = 3 near coefficients at b = 0 against Math-#184 (floating, from [NUM] GH60 times R_3(0))') + hd = {0.5: 4.8558, 1.0: 5.705} + for k, v in hd.items(): + iv = [float(x) for x in D3['k=%g, b=0' % k]['a_fail^(3)']] + print(' k=%g a_fail^(3) certified [%.6f, %.6f] Math-#184 %.4f (%+.3f%%)' % (k, iv[0], iv[1], v, 100 * (v / (0.5 * (iv[0] + iv[1])) - 1))) + print('6. Lemma T: direct floating lattice sums of the torus deviation at L = 10 (three frames) against delta') + import torus as TOR + dlt = TOR.delta() + def He(j, t): + h0, h1 = 1.0, t + if j == 0: + return h0 + for i in range(1, j): + h0, h1 = h1, t * h1 - i * h0 + return h1 + worst = 0.0 + for th in (0.0, 0.3, math.pi / 4): + c, s_ = math.cos(th), math.sin(th) + for g in ((6, 0), (3, 3), (4, 2), (2, 2), (5, 1), (0, 0)): + tot = 0.0 + theta = 0.0 + for n1 in range(-3, 4): + for n2 in range(-3, 4): + if n1 == 0 and n2 == 0: + continue + y1, y2 = 10.0 * (c * n1 + s_ * n2), 10.0 * (-s_ * n1 + c * n2) + e = math.exp(-(y1 * y1 + y2 * y2) / 2) + tot += (-1) ** (g[0] + g[1]) * He(g[0], y1) * He(g[1], y2) * e + theta += e + d0 = TOR._D(g) + dev = abs((tot - d0 * theta) / (1.0 + theta)) # (d0 + tot)/(1 + theta) - d0 without cancellation + worst = max(worst, dev) + print(' max floating deviation %.3e <= delta %.3e: %s' % (worst, dlt, worst <= dlt)) + +if __name__ == '__main__': + main() diff --git a/frontiers/c8_elder_failure_coefficient_20261001/ia.py b/frontiers/c8_elder_failure_coefficient_20261001/ia.py new file mode 100644 index 000000000..cf4bbfa16 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/ia.py @@ -0,0 +1,279 @@ +"""Interval arithmetic on IEEE-754 binary64: every +, -, *, / and sqrt is correctly rounded to nearest and the result is +widened outward by one ulp (math.nextafter), so each interval contains the exact real result. No libm transcendental is +trusted: pi and ln 2 come from exact rational series, exp(-y) from a reduced Taylor polynomial with an a priori error +bound, Phi from its positive power series and Q = 1 - Phi (t >= 2) from Laplace's continued fraction for the Mills ratio. +libm values (pow, acos, cos) are used only as starting points of iterations whose results are then verified. +An interval is a pair (lo, hi) of floats with lo <= hi.""" +import math +from fractions import Fraction as Fr + +INF = math.inf +_na = math.nextafter + + +def dn(x): + return _na(x, -INF) + + +def up(x): + return _na(x, INF) + + +def iv(x): + """Exact float point as an interval.""" + return (x, x) + + +def from_frac(q): + """Tightest float interval containing the rational q.""" + f = float(q) + fq = Fr(f) + if fq == q: + return (f, f) + return (f, up(f)) if fq < q else (dn(f), f) + + +def add(a, b): + return (dn(a[0] + b[0]), up(a[1] + b[1])) + + +def sub(a, b): + return (dn(a[0] - b[1]), up(a[1] - b[0])) + + +def neg(a): + return (-a[1], -a[0]) + + +def mul(a, b): + p = (a[0] * b[0], a[0] * b[1], a[1] * b[0], a[1] * b[1]) + return (dn(min(p)), up(max(p))) + + +def div(a, b): + if b[0] <= 0.0 <= b[1]: + raise ZeroDivisionError('interval division by an interval containing 0') + p = (a[0] / b[0], a[0] / b[1], a[1] / b[0], a[1] / b[1]) + return (dn(min(p)), up(max(p))) + + +def scal(c, a): + """c * a for an exact float c.""" + if c >= 0: + return (dn(c * a[0]), up(c * a[1])) + return (dn(c * a[1]), up(c * a[0])) + + +def divc(a, c): + """a / c for an exact float c > 0.""" + return (dn(a[0] / c), up(a[1] / c)) + + +def sqr(a): + if a[0] >= 0: + return (max(0.0, dn(a[0] * a[0])), up(a[1] * a[1])) + if a[1] <= 0: + return (max(0.0, dn(a[1] * a[1])), up(a[0] * a[0])) + m = max(-a[0], a[1]) + return (0.0, up(m * m)) + + +def pw(a, n): + """a ** n for an integer n >= 0 (monotone pieces).""" + if n == 0: + return (1.0, 1.0) + if n == 1: + return a + if n % 2 == 1: + # odd powers are increasing: bound each endpoint by repeated outward multiplication + return (_pw_pt(a[0], n, False), _pw_pt(a[1], n, True)) + r = pw(sqr(a), n // 2) + return (max(0.0, r[0]), r[1]) + + +def _pw_pt(x, n, upper): + """Outward bound of x ** n at a point (odd n): upper -> rounded up, else rounded down.""" + r = (x, x) + for _ in range(n - 1): + r = mul(r, (x, x)) + return r[1] if upper else r[0] + + +def hull(*ivs): + return (min(x[0] for x in ivs), max(x[1] for x in ivs)) + + +def absmax(a): + return max(-a[0], a[1]) + + +def isqrt(a): + if a[0] < 0: + raise ValueError('sqrt of a negative interval') + lo = math.sqrt(a[0]) + hi = math.sqrt(a[1]) + return (max(0.0, dn(lo)), up(hi)) + + +def contains(a, x): + return a[0] <= x <= a[1] + + +# ---------------------------------------------------------------- constants (exact rationals, then outward floats) +def _pi_frac(digits=40): + """pi by Machin's formula with an explicit rational enclosure.""" + def arctan_inv(n, terms): + s = Fr(0) + for k in range(terms): + s += Fr((-1) ** k, (2 * k + 1) * n ** (2 * k + 1)) + return s + T = 40 + lo = 16 * arctan_inv(5, T) - 4 * arctan_inv(239, T) + err = Fr(16, (2 * T + 1) * 5 ** (2 * T + 1)) + Fr(4, (2 * T + 1) * 239 ** (2 * T + 1)) + return lo - err, lo + err + + +_PI_LO, _PI_HI = _pi_frac() +PI = (from_frac(_PI_LO)[0], from_frac(_PI_HI)[1]) +SQRT2PI = isqrt(scal(2.0, PI)) +SQRT2 = isqrt((2.0, 2.0)) +INV_SQRT2PI = div((1.0, 1.0), SQRT2PI) + + +# ---------------------------------------------------------------- exp(-y) for y >= 0 +_INVFACT = [1.0] +for _i in range(1, 19): + _INVFACT.append(_INVFACT[-1] / _i) +_U = 2.0 ** -53 + + +def _ln2_frac(): + """ln 2 = sum_{j>=1} 1/(j 2^j), with the tail after J terms below 2/(J 2^J).""" + J = 80 + s = Fr(0) + for j in range(1, J + 1): + s += Fr(1, j * 2 ** j) + return s, s + Fr(2, J * 2 ** J) + + +_L2LO, _L2HI = _ln2_frac() +LN2 = (from_frac(_L2LO)[0], from_frac(_L2HI)[1]) + + +def _poly_exp_neg(z): + """Binary64 Horner value p of the degree-18 Taylor polynomial of exp(-z) at a float z with |z| <= 0.75, and a bound + of |p - exp(-z)|: Horner and coefficient rounding (coefficients 1/i! carry relative error <= 18u) give at most + (2*18 + 18 + 2) u e^|z| (1 + 1e-12), truncation at most |z|^19/19! e^|z| < 1e-19 e^|z|; exp(-z) >= e^-|z| >= 0.47.""" + p = _INVFACT[18] + for i in range(17, -1, -1): + p = _INVFACT[i] - z * p + ez = 2.12 # >= e^0.75 + err = (56.0 * _U * (1.0 + 1e-12) + 1e-19) * ez + return p, err + + +def _exp_neg_point(y): + """Enclosure of exp(-y) for a float y >= 0: y = n ln2 + r with n = floor(y / ln2) and r in an interval inside + [-1e-9, 0.75] (from the rational enclosure of ln 2), exp(-y) = 2^-n exp(-r) (exact power-of-two scaling).""" + if y == 0.0: + return (1.0, 1.0) + if y > 700.0: + return (0.0, math.ldexp(1.0, -1000)) # exp(-700) < 2^-1000 + n = int(y / 0.6931471805599453) + nL = mul((float(n), float(n)), LN2) + r = sub((y, y), nL) + if not (-1e-9 <= r[0] <= r[1] <= 0.75): + raise RuntimeError('exp range reduction failed at y=%r' % y) + plo, elo = _poly_exp_neg(r[1]) # exp(-r) is decreasing: lower end at r[1] + phi_, ehi = _poly_exp_neg(r[0]) + lo = dn(plo - elo) + hi = up(phi_ + ehi) + return (max(0.0, math.ldexp(lo, -n)), math.ldexp(hi, -n)) + + +def exp_neg(a): + """Enclosure of exp(-y) for y in the interval a, a[0] >= 0 (decreasing).""" + if a[0] < 0: + # a lower bound that is negative only by underflow-level rounding (y >= -1e-300): exp(-y) <= exp(1e-300) < up(1) + if a[0] < -1e-300: + raise ValueError('exp_neg needs y >= 0') + return (_exp_neg_point(max(0.0, a[1]))[0], up(1.0)) + return (_exp_neg_point(a[1])[0], _exp_neg_point(a[0])[1]) + + +def phi_std(t): + """Standard normal density at the interval t.""" + return mul(exp_neg(divc(sqr(t), 2.0)), INV_SQRT2PI) + + +def _Phi_pos_point(t): + """Enclosure of Phi(t) for a float 0 <= t < 2: Phi(t) = 1/2 + phi(t) S(t), S(t) = sum t^(2n+1)/(2n+1)!! (positive + terms), summed in binary64. On 0 < t < 2 the loop stops by n = 24 (there term_n/S <= 4^n/(2n+1)!! < 5e-18), so the + accumulated relative error is at most (3 * 24 + 3) u < 1e-14, inside the 1e-13 allowance; the cap n <= 300 only + guards termination. The geometric tail bound applies once the term ratio t^2/(2n+3) is below 1/2.""" + if t == 0.0: + return (0.5, 0.5) + if t < 1e-100: + return (0.5, up(0.5 + t)) + if not t < 2.0: + raise ValueError('series branch is for t < 2') + t2 = t * t + term = t + S = t + n = 0 + while True: + n += 1 + if n > 300: + raise RuntimeError('Phi series did not terminate at t=%r' % t) + term = term * t2 / (2 * n + 1) + S += term + ratio = t2 * (1.0 + 4 * _U) / (2 * n + 3) + if ratio < 0.5 and term < 1e-17 * S: + tail = term * ratio / (1.0 - ratio) * (1.0 + 1e-12) + break + Slo, Shi = dn(S * (1.0 - 1e-13)), up((S + tail) * (1.0 + 1e-13)) + r = add((0.5, 0.5), mul(phi_std((t, t)), (Slo, Shi))) + return (r[0], min(r[1], 1.0)) + + +def _mills(t): + """Enclosure of the Mills ratio R(t) = Q(t)/phi(t), t >= 2 float: Laplace's continued fraction + R(t) = 1/(t + 1/(t + 2/(t + 3/(t + ...)))) has positive partial numerators and denominators, so consecutive + approximants bracket R(t) (Stieltjes); the hull of the depth-n and depth-(n+1) approximants is returned once its + relative width is below 1e-14 (each approximant evaluated in outward interval arithmetic).""" + ti = (t, t) + n = 16 + while True: + res = [] + for depth in (n, n + 1): + v = ti + for j in range(depth, 0, -1): + v = add(ti, div((float(j), float(j)), v)) + res.append(div((1.0, 1.0), v)) + h = (min(res[0][0], res[1][0]), max(res[0][1], res[1][1])) + if h[1] - h[0] <= 1e-14 * h[0] or n > 4000: + return h + n *= 2 + + +def Q_point(t): + """Enclosure of the upper tail Q(t) = 1 - Phi(t) for a float t >= 0, with relative precision for large t.""" + if t >= 2.0: + return mul(phi_std((t, t)), _mills(t)) + r = _Phi_pos_point(t) + return (max(0.0, dn(1.0 - r[1])), up(1.0 - r[0])) + + +def Phi_point(t): + if t >= 0: + if t >= 2.0: + q = Q_point(t) + return (dn(1.0 - q[1]), min(1.0, up(1.0 - q[0]))) + return _Phi_pos_point(t) + q = Q_point(-t) + return q + + +def Phi(a): + """Normal CDF over an interval argument (increasing).""" + return (Phi_point(a[0])[0], Phi_point(a[1])[1]) diff --git a/frontiers/c8_elder_failure_coefficient_20261001/torus.py b/frontiers/c8_elder_failure_coefficient_20261001/torus.py new file mode 100644 index 000000000..4a448b2f3 --- /dev/null +++ b/frontiers/c8_elder_failure_coefficient_20261001/torus.py @@ -0,0 +1,219 @@ +"""Lemma T (torus transfer) for CL-C8-ELDER-FAILURE-COEFFICIENT-20261001: the planar coefficients of the torus model of +[LP] / [S] (side L >= 10, any frame) against the reference-kernel ones. + +One-site jet covariances. Cov(d^a f, d^b f) = (-1)^|b| d^(a+b) K(0). For the torus kernel +K_L(z) = sum_n G(z + Ln) / Theta_L, G = exp(-|z|^2/2), Theta_L = sum_n exp(-L^2 |n|^2/2) ([LP] section 2), in a rotated +frame d^g K_L(0) = sum_n (d^g G)(R^T L n) / Theta_L with |R^T L n| = L|n|, and |d^g G(y)| = |He_g1(y1) He_g2(y2)| G(y) +<= Hh_g1(|y|) Hh_g2(|y|) G(y), Hh_j(t) = sum_k j!/(k!(j-2k)! 2^k) t^(j-2k). Grouping n != 0 by |n|_inf = j (8j points, +|n| >= j) and using that t -> Hh_a(t) Hh_c(t) e^(-t^2/2) decreases for t^2 > a + c, every entry of order |a| + |b| <= 6 +differs from the reference by at most delta(10) for every L >= 10 and every frame. + +Propagation. Exact interval regressions with every entry widened by delta give: the odd conditional law +N(12k v, S_c) of (f_xxz, f_xzz, f_zzz) given (f_x, f_z, f_xxx) = (0, 0, 12k) (reference: N(0, diag(2, 2, 6))); the even +conditional law N(c1 b, s^2) of f_zz given (f, f_xx, f_xz) = (b, 0, 0) (reference N(-b, 2)); and the factors of +pi_0 = p_(f,V_u)(b, 0) p_G(0) phi_tau(12k) ([LP] section 15). Parity (K_L even) keeps odd and even jets uncorrelated. + +Expectations. With S_c = S0^(1/2) (I + E) S0^(1/2), ||E|| <= eps, mu~ = S0^(-1/2) mu: + 1 + chi^2(N(mu, S_c) || N(0, S0)) = prod (1 - e_i^2)^(-1/2) exp(sum mu~_i^2/(1 - e_i)) <= (1 - eps^2)^(-3/2) exp(|mu~|^2/(1 - eps)), +and for eps^2 <= 0.006, where (1 - eps^2)^(-3/2) <= 1.01 and (1 - eps^2)^(-3/2) - 1 <= 2 eps^2, with e^m - 1 <= m e^m, + chi^2 <= 1.01 m e^m + 2 eps^2 (m = |mu~|^2/(1 - eps)), +and |E_L[I] - E_ref[I]| <= (E_ref[I^2] chi^2)^(1/2) +(Cauchy-Schwarz), with I = 3k^2 H <= 12 Pi_4 for k <= 2 (Pi_4 the tail polynomial of certificate.tail4).""" +import math +from fractions import Fraction as Fr +from functools import reduce +from ia import (dn, up, add, sub, neg, mul, div, scal, divc, sqr, pw, isqrt, exp_neg, from_frac, absmax, Phi, SQRT2PI) + +L_MIN = 10.0 +ORDER = 6 + + +def _Hh(j): + return {j - 2 * k: Fr(math.factorial(j), math.factorial(k) * math.factorial(j - 2 * k) * 2 ** k) for k in range(j // 2 + 1)} + + +def _poly(coef, t): + acc = (0.0, 0.0) + for p, c in coef.items(): + acc = add(acc, mul(from_frac(c), pw(t, p))) + return acc + + +def _M(t): + """max over a + c <= ORDER of Hh_a(t) Hh_c(t) (upper end), t >= 1.""" + best = 0.0 + for m in range(ORDER + 1): + for a in range(m + 1): + best = max(best, mul(_poly(_Hh(a), t), _poly(_Hh(m - a), t))[1]) + return best + + +def delta(L=L_MIN, J=6): + """Bound on |c_L - c_ref| for every one-site covariance entry of order <= ORDER, every side >= L and every frame.""" + S = 0.0 + term = 0.0 + for j in range(1, J + 1): + t = (L * j, L * j) + term = up(8.0 * j * up(_M(t) * exp_neg(divc(sqr(t), 2.0))[1])) + S = up(S + term) + S = up(S + 2.0 * term) # shells j > J: consecutive bounds decrease by far more than 1/2 + theta_minus_1 = S # sum_{n != 0} e^(-L^2 |n|^2/2) <= S since Hh products are >= 1 + return up(S + up(15.0 * theta_minus_1)) # |d^g G(0)| <= 15 for |g| <= 6; |1/Theta - 1| <= Theta - 1 + + +def _dfact(n): + r = 1 + while n > 1: + r, n = r * n, n - 2 + return r + + +def _D(g): + if g[0] % 2 or g[1] % 2: + return 0 + return (-1) ** ((g[0] + g[1]) // 2) * _dfact(g[0] - 1) * _dfact(g[1] - 1) + + +def cov_ref(a, b): + """Reference covariance of the jets d^a f and d^b f at one site (multi-indices (x-order, z-order)).""" + return (-1) ** (b[0] + b[1]) * _D((a[0] + b[0], a[1] + b[1])) + + +def _icov(rows, cols, d): + return [[(dn(cov_ref(a, b) - d), up(cov_ref(a, b) + d)) for b in cols] for a in rows] + + +def _det2(M): + return sub(mul(M[0][0], M[1][1]), mul(M[0][1], M[1][0])) + + +def _det3(M): + return sub(add(mul(M[0][0], sub(mul(M[1][1], M[2][2]), mul(M[1][2], M[2][1]))), + mul(M[0][2], sub(mul(M[1][0], M[2][1]), mul(M[1][1], M[2][0])))), + mul(M[0][1], sub(mul(M[1][0], M[2][2]), mul(M[1][2], M[2][0])))) + + +def _inv3(M): + d = _det3(M) + adj = [[None] * 3 for _ in range(3)] + for i in range(3): + for j in range(3): + r = [x for x in range(3) if x != i] + c = [x for x in range(3) if x != j] + m = sub(mul(M[r[0]][c[0]], M[r[1]][c[1]]), mul(M[r[0]][c[1]], M[r[1]][c[0]])) + adj[j][i] = m if (i + j) % 2 == 0 else neg(m) + return [[div(adj[i][j], d) for j in range(3)] for i in range(3)] + + +def _inv2(M): + d = _det2(M) + return [[div(M[1][1], d), div(neg(M[0][1]), d)], [div(neg(M[1][0]), d), div(M[0][0], d)]] + + +def _mm(A, B): + return [[reduce(add, [mul(A[i][t], B[t][j]) for t in range(len(B))]) for j in range(len(B[0]))] for i in range(len(A))] + + +def _T(A): + return [list(r) for r in zip(*A)] + + +def _expiv(x): + """exp over a small interval x (|x| < 1) from the verified exp_neg of ia: exp(t) = exp_neg(-t) for t <= 0 and + 1/exp_neg(t) for t > 0.""" + def e(t): + return exp_neg((-t, -t)) if t <= 0 else div((1.0, 1.0), exp_neg((t, t))) + return (e(x[0])[0], e(x[1])[1]) + + +def transfer(E_I2_upper): + """Lemma T bounds. E_I2_upper: an upper bound of E_ref[I^2], I = 3k^2 H, uniform over k in [1/2, 2]. + Returns a dict of bounds (all floats, upper bounds unless stated).""" + d = delta() + # ---- odd block: Y = (f_x, f_z, f_xxx) pinned at (0, 0, 12k); X = (f_xxz, f_xzz, f_zzz) + Y = [(1, 0), (0, 1), (3, 0)] + X = [(2, 1), (1, 2), (0, 3)] + SYY, SXY, SXX = _icov(Y, Y, d), _icov(X, Y, d), _icov(X, X, d) + Bc = _mm(SXY, _inv3(SYY)) + BS = _mm(Bc, _T(SXY)) + Sc = [[sub(SXX[i][j], BS[i][j]) for j in range(3)] for i in range(3)] + S0 = (2.0, 2.0, 6.0) + eps2 = (0.0, 0.0) # Frobenius bound of E = S0^(-1/2) (S_c - S0) S0^(-1/2) + for i in range(3): + for j in range(3): + dev = sub(Sc[i][j], ((S0[i] if i == j else 0.0),) * 2) + e = div((absmax(dev),) * 2, isqrt((S0[i] * S0[j],) * 2)) + eps2 = add(eps2, sqr(e)) + eps2 = eps2[1] + eps = isqrt((eps2, eps2))[1] + mu2 = (0.0, 0.0) # |mu~|^2, mu = 12k Bc[:, 2] with k <= 2 + for i in range(3): + m_i = scal(24.0, (absmax(Bc[i][2]),) * 2) + mu2 = add(mu2, divc(sqr(m_i), S0[i])) + mu2 = mu2[1] + if not eps2 <= 0.006: # the range of the chi^2 bound (module docstring) + raise RuntimeError('torus transfer: perturbation too large') + m = div((mu2, mu2), sub((1.0, 1.0), (eps, eps)))[1] + chi2 = add(mul(from_frac(Fr(101, 100)), mul((m, m), _expiv((m, m)))), scal(2.0, (eps2, eps2)))[1] + dJ = isqrt(mul((E_I2_upper, E_I2_upper), (chi2, chi2)))[1] + # ---- even block: W = (f, f_xx, f_xz) pinned at (b, 0, 0); A = f_zz + W = [(0, 0), (2, 0), (1, 1)] + A = [(0, 2)] + SWW, SAW, SAA = _icov(W, W, d), _icov(A, W, d), _icov(A, A, d) + Wi = _inv3(SWW) + cA = _mm(SAW, Wi)[0] # regression row: mean of A = cA[0] b + sA2 = sub(SAA[0][0], _mm(_mm(SAW, Wi), _T(SAW))[0][0]) + # p_b(0)^L / p_b(0)^ref = (sqrt2/s_L) exp(-b^2 (cA0^2/(2 s_L^2) - 1/4)), b in [0, 1] + x = sub(div(sqr(cA[0]), scal(2.0, sA2)), (0.25, 0.25)) + ax = absmax(x) + rp = mul(div(isqrt((2.0, 2.0)), isqrt(sA2)), _expiv((-ax, ax))) + # pi_0 factors: p_W(b,0,0) ratio, p_G(0) ratio, phi_tau(12k) ratio + det_ref_W = (2.0, 2.0) # det [[1,-1,0],[-1,3,0],[0,0,1]] + y = sub(Wi[0][0], (1.5, 1.5)) + ay = absmax(y) + rW = mul(isqrt(div(det_ref_W, _det3(SWW))), _expiv((-0.5 * ay, 0.5 * ay))) + G = [(1, 0), (0, 1)] + SGG = _icov(G, G, d) + rG = isqrt(div((1.0, 1.0), _det2(SGG))) + T3 = [(3, 0)] + S3G, S33 = _icov(T3, G, d), _icov(T3, T3, d) + tau2 = sub(S33[0][0], _mm(_mm(S3G, _inv2(SGG)), _T(S3G))[0][0]) + z = sub(div((1.0, 1.0), tau2), div((1.0, 1.0), (6.0, 6.0))) + az = up(288.0 * absmax(z)) # 72 k^2 |1/tau_L^2 - 1/6|, k <= 2 + rT = mul(isqrt(div((6.0, 6.0), tau2)), _expiv((-az, az))) + rho = mul(mul(mul(rW, rG), rT), rp) # (pi_0 p_b(0))^L / (pi_0 p_b(0))^ref, uniform in (b, k) + return {'delta': d, 'eps': eps, 'mu2': mu2, 'chi2': chi2, 'dJ': dJ, 'mean_coeff_A': cA[0], 'var_A': sA2, + 'p_ratio': rp, 'pi0p_ratio': rho, 'tau2': tau2, 'S_c': Sc} + + + +def abs_moment(s2, n): + """E|X|^n for X ~ N(0, s2) with s2 an exact float: n even (n-1)!! s2^(n/2); n odd 2^(n/2) Gamma((n+1)/2)/sqrt(pi) s^n, + i.e. (n-1)!! s^n sqrt(2/pi) (enclosed).""" + if n % 2 == 0: + return pw(scal(float(_dfact(n - 1)), (1.0, 1.0)), 1) if n == 0 else mul((float(_dfact(n - 1)),) * 2, pw((s2, s2), n // 2)) + from ia import PI + s = isqrt((s2, s2)) + return mul(mul((float(_dfact(n - 1)),) * 2, pw(s, n)), isqrt(div((2.0, 2.0), PI))) + + +def E_I2_upper(): + """E_ref[I^2] <= 144 E[Pi_4^2], I = 3k^2 H <= 12 Pi_4 for k in [1/2, 2], (a, beta, c) ~ N(0, diag(2, 2, 6)), + Pi_4 = 68|beta|^3 + c6 a^6 + 144 c^2 + 9 a^2 beta^2, c6 = 68/216 + 576/5184.""" + c6 = add(divc((68.0, 68.0), 216.0), divc((576.0, 576.0), 5184.0)) + terms = [((68.0, 68.0), (0, 3, 0)), (c6, (6, 0, 0)), ((144.0, 144.0), (0, 0, 2)), ((9.0, 9.0), (2, 2, 0))] + var = (2.0, 2.0, 6.0) + tot = (0.0, 0.0) + for c1, e1 in terms: + for c2, e2 in terms: + v = mul(c1, c2) + for j in range(3): + v = mul(v, abs_moment(var[j], e1[j] + e2[j])) + tot = add(tot, v) + return up(144.0 * tot[1]) + + +if __name__ == '__main__': + t = transfer(E_I2_upper()) + for kk in ('delta', 'eps', 'mu2', 'chi2', 'dJ', 'mean_coeff_A', 'var_A', 'p_ratio', 'pi0p_ratio', 'tau2'): + print(kk, t[kk])