diff --git a/.github/workflows/cusp-coefficient-certified.yml b/.github/workflows/cusp-coefficient-certified.yml new file mode 100644 index 00000000..26db78d2 --- /dev/null +++ b/.github/workflows/cusp-coefficient-certified.yml @@ -0,0 +1,80 @@ +name: Cusp coefficient c1 (certified replay) + +'on': + pull_request: + paths: + - 'frontiers/cusp_coefficient_certified_20261001/**' + - '.github/workflows/cusp-coefficient-certified.yml' + - 'coefficients/side24_v1/PROOF.md' + - 'coefficients/side24_v1/ENCLOSURE.json' + - 'imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md' + - 'frontiers/cusp_second_order_20261001/PROOF.md' + workflow_dispatch: + +permissions: + contents: read + +jobs: + certified-replay: + runs-on: ubuntu-latest + timeout-minutes: 90 + 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/cusp_coefficient_certified_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/cusp_coefficient_certified_20261001 + python $FLAGS certificate.py --check --procs "$(nproc)" + for m in kernel-A gl-weight bound-scale phase-sign tail-drop image-order6; do + set +e + python $FLAGS certificate.py --mutant "$m" > /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/cusp_coefficient_certified_20261001/NOTE.md b/frontiers/cusp_coefficient_certified_20261001/NOTE.md new file mode 100644 index 00000000..b2866871 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/NOTE.md @@ -0,0 +1,391 @@ +# Certified cusp coefficient `c1` of Theorem CU (Math-#207): Gaussian kernel, `d = 1, 2, 3`; [P] torus field, `L >= 24` + +**Object:** CL-CU-CUSP-COEFFICIENT-20261001-v1. +**Author:** Anthropic Claude (Claude Code session `session_015wNj8LPTKXsaT68G3DgPPh`), 1 October 2026. +**Disposition:** certified numerical enclosures of the explicit Gaussian integral (CU.2) of Math-#207, for the Gaussian +kernel (v1) and, by Lemma S, for the [P] torus field of every side `L >= 24`, including SIDE24 (v1.1). Author-side. Nonauthor +reads: OpenAI accepts R.1 (review 5378866830 on `357d0a1`) and R.2–R.6 with Lemma Q (review 5384432735 on `7a04873`). Scope claim: Math-#207 comment +5929273376. +**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. Nothing in Math-#207, #214, #216 or #218 is touched. +**Delivered under:** Dylan Roy's explicit instructions in this session, given after the 2026-09-27 owner stop (quoted in +`SOURCE_MAP.json`, `delivered_under`). `OWNER_STOP.md` addresses Cursor agents and automations and is untouched. + +## 0. Statement + +Math-#207 ([CU], `PROOF.md` blob `f6df5a73`, merged at `566b1a1`; merging is not a status promotion) proves, at candidate +status, the second-order term of the +elder lifetime law, `nu_eld(ell) = c ell^(-1/3) + c1 ell^(1/4) + o(ell^(1/4))`, with + + c1 = -(192/7) 2^(1/4) integral_(S^(d-1)) integral_R pi_0(u; v_0(b,0)) E_0[ |Y|^(7/4) |Delta|^(1/4) 1{A < 0} | b ] db dsigma(u), (CU.2) + +where `Y = (f4/12) Delta - gamma^T adj(A) gamma/4` and `Delta = det A`, built from the zero-gap contact jets `f4 = d_u^4 f(0)`, +`gamma = grad_Theta d_u^2 f(0)` and `A = D_Theta^2 f(0)` ([CU] (0.1)). [CU] evaluates (CU.2) by floating quadrature and lists +certified numerics as not claimed. This packet certifies it. + +**Theorem (certified enclosures).** For the Gaussian kernel `K(z) = exp(-|z|^2/2)` on `R^d`, the integral (CU.2) and the +candidate coefficient `I^cand = (3^(1/4)/2) c1` of [CU] (CU'.1) satisfy + +| `d` | `c1` | `I^cand` | `I^cand - c1` | +|---|---|---|---| +| `1` | `[-0.22760635877558, -0.22760635877504]` | `[-0.14977340698364, -0.14977340698328]` | `[0.077832951791707, 0.077832951791985]` | +| `2` | `[-0.269398825674, -0.269398825672]` | `[-0.177274396795, -0.177274396793]` | `[0.0921244288786, 0.0921244288791]` | +| `3` | `[-0.211848347, -0.211848346]` | `[-0.139404052, -0.139404051]` | `[0.0724442943, 0.0724442950]` | + +With the leading coefficients `c_(d,ref)` of `coefficients/side24_v1` (1) ([SIDE24]), enclosed here from their closed +form, + +| `d` | `c_(d,ref)` | `c1 / c_(d,ref)` | +|---|---|---| +| `2` | `[0.0734069193059, 0.0734069193062]` | `[-3.66993776909, -3.66993776907]` | +| `3` | `[0.0417759318405, 0.0417759318407]` | `[-5.07106215, -5.07106213]` | + +**Corollary S (the [P] torus field; v1.1).** Let `L >= 24`, and let the field on `R^d/(L Z^d)` have covariance +`K_L(z) = sum_(n in Z^d) exp(-|z + L n|^2/2)`, normalized to variance one or not. In particular this covers SIDE24 +(`L = 24`), the manuscript's field ([CU] §0). Then for `d = 2, 3` the coefficient (CU.2) of this field satisfies +`|c1^(L)/c1^(ref) - 1| <= eta_d`, with `eta_2 <= 5.9e-105` and `eta_3 <= 1.4e-104` (Lemma S, §3b). Hence + +| `d` | `c1^(L)` | `I^cand,(L)` | `c1^(24) / c_(d,24)` | +|---|---|---|---| +| `2` | `[-0.269398825674, -0.269398825672]` | `[-0.177274396795, -0.177274396793]` | `[-3.66993776909, -3.66993776907]` | +| `3` | `[-0.211848347, -0.211848346]` | `[-0.139404052, -0.139404051]` | `[-5.07106215, -5.07106213]` | + +Here `c_(d,24)` is taken from `coefficients/side24_v1/ENCLOSURE.json` ([SIDE24-E], consumed): `c_(2,24)` in +`[0.07340691930603427103, 0.07340691930603427104]` and `c_(3,24)` in `[0.04177593184059834334, 0.04177593184059834335]`. + +`RESULTS.json` holds the certificate's own intervals (16 digits, after widening by `1e-12` relative); the tables round them +further outward. + +**Consequences.** +- **[CU]'s floating values are confirmed to all their digits.** [CU] §8 (remark 3) reports `c1 = -0.26939883` (`d = 2`) + and `-0.21184835` (`d = 3`), `c = 0.0734069193` and `0.0417759318`, `c1/c = -3.669938` and `-5.071062`, and `I^cand`, + `I^cand - c1` to seven digits. Each is the correct rounding of the certified midpoint (`controls.py`, §7). The `d = 3` + referee Monte Carlo `-0.211848 +- 0.000024` is `0.01` standard errors from it. +- **The manuscript's correction ratio is certified.** For the manuscript's field (`d = 3`, `L = 24`) the relative + correction of [CU] (CU.1) is `(c1/c) ell^(7/12)`, with `c1/c_(3,24)` in `[-5.07106215, -5.07106213]` (Corollary S). [CU] §0 + states `c1/c ~ -5.071` for this field. +- **`d = 1` is a closed form.** `c1^(1) = -(384/7) 2^(1/4) 12^(-9/4) (2 pi)^-2 sqrt(4 pi/3) 30^(7/8) 2^(7/8) + Gamma(11/8)/sqrt(pi)`. It agrees with Math-#214's `C_1` (#216's table: `-0.227606`). As [CU]'s author notes (#207 comment + 5925314397), (CU.2) at `d = 1` is #214's coefficient. + +## 1. What is evaluated + +The integral (CU.2) as [CU] defines it. +- `E_0[. | b]` is the expectation over the contact jets under the zero-gap target `v_0(b, 0)`: `f(0) = b`, + `d_u f(0) = d_u^2 f(0) = d_u^3 f(0) = 0`, `grad_Theta f(0) = grad_Theta d_u f(0) = 0`. +- `pi_0(u; v_0(b, 0))` is the density of those pins. By [LP] §15 and [CU] (CU.2) it factorizes as + `pi_0 = p_(f, V_u)(b, 0) p_G(0) phi_tau(0)`, where `G = grad f(0)` and `V_u = H_f(0) u`. +- `d sigma` is the surface measure on `S^(d-1)`, with `|S^0| = 2`, `|S^1| = 2 pi`, `|S^2| = 4 pi`. + +The kernel is `K(z) = exp(-|z|^2/2)`, which is isotropic, so the `u`-integral is `|S^(d-1)|` times the value at `u = e_1`. +`Cov(d^a f, d^b f) = (-1)^|b| d^(a+b) K(0)`, and `d^g K(0) = prod_i (-1)^(g_i/2) (g_i - 1)!!` for even `g`, 0 otherwise. + +## 2. Reduction (Lemma R) + +**R.1 (the zero-gap contact law).** Given the pins, +- `f4 ~ N(-3b, 24)`; +- `gamma ~ N(0, 2 I_(d-1))`; +- `A = -b I_(d-1) + G`, with `G` symmetric, `G_jj ~ N(0, 2)` and `G_jk ~ N(0, 1)` for `j < k`; +- all of these independent; +- `pi_0 = (2 pi)^(-(d+1)) 12^(-1/2) e^(-3b^2/4)`. + +*Proof.* +- **Odd jets.** `(f_x, f_y, f_xxx, gamma)` are uncorrelated with the even jets (`K` is even). + - Each `gamma_j = f_xxy_j` has variance 3 and covariance `-1` with `f_y_j`. Its covariances with `f_x`, `f_xxx` and the + other `gamma`'s vanish, so given the pins its variance is `3 - 1 = 2`. + - `tau^2 = Var(f_xxx | grad f = 0) = 15 - 9 = 6`, so `phi_tau(0) = (12 pi)^(-1/2)`, and `p_G(0) = (2 pi)^(-d/2)`. +- **Even jets.** `W = (f, f_xx, f_xy_j)` has covariance `[[1, -1], [-1, 3]] (+) I_(d-1)`, with + `W^(-1)_(11) = 3/2` and determinant 2, so `p_W(b, 0) = (2 pi)^(-(d+1)/2) 2^(-1/2) e^(-3b^2/4)`. + - Regression on `W = (b, 0, 0)` gives `E[f4] = (3, -15).(3/2, 1/2) b = -3b` and `Var f4 = 105 - 81 = 24`. + - It gives `E[A_jj] = -b`, `Var A_jj = 3 - 1 = 2`, `Cov(A_jj, A_kk) = 1 - 1 = 0` and `Var A_jk = 1`. + - It gives `Cov(f4, A_jj) = -3 - (-3) = 0`. + - All other covariances vanish by parity of the individual indices. +- Multiplying the three factors gives `pi_0`. ∎ + +**R.2 (Mellin form).** For `0 < p < 2`, `|y|^p = K^(-1) integral_0^oo (1 - cos(omega y)) omega^(-1-p) d omega`, with +`K = pi/(2 Gamma(1+p) sin(pi p/2))`. At `p = 7/4`, `K = pi/(2 Gamma(11/4) sin(7 pi/8))`. + +Tonelli applies to the real nonnegative function `(1 - cos(omega Y)) omega^(-11/4)` times the nonnegative determinant and +cone weight. At each fixed `omega`, the bounded characteristic function is then integrated against the finite Gaussian +weight. The two terms are never split into separate, divergent, outer Mellin integrals. Hence +`E|Y|^(7/4) = K^(-1) integral (1 - Re E e^(i omega Y)) omega^(-11/4) d omega` (wording per 5384432735). + +**R.3 (`d = 1`).** Here `A` is empty and `Y = f4/12`. The `b`-weight `e^(-3b^2/4)` makes `f4` a centered normal of variance +`24 + 9 (2/3) = 30`. With `E|Z|^(7/4) = 2^(7/8) Gamma(11/8)/sqrt(pi)` this gives the closed form of §0. + +**R.4 (`d = 2`).** On `A = -a < 0` we have `Delta = A` and `adj A = 1`. +- Given `(a, b)`, `(f4/12) A ~ N(ab/4, a^2/6)` and `E e^(-i omega gamma^2/4) = (1 + i omega)^(-1/2)`. +- The `b`-weight `pi_0 phi_(-b,2)(-a)` equals `(2 pi)^-3 12^(-1/2) (4 pi)^(-1/2) e^(-3a^2/16) e^(-(b - a/4)^2)`, and + `integral db` gives the factor `sqrt(pi) e^(i omega a^2/16 - omega^2 a^2/64)`. +- Then `integral_0^oo a^(1/4) e^(-beta a^2) da = Gamma(5/8)/(2 beta^(5/8))`, for `Re beta > 0` and the principal power, + with `beta(omega) = 3/16 - i omega/16 + 19 omega^2/192`. +- Hence + c1^(2) = -(192/7) 2^(1/4) (2 pi)^-2 12^(-1/2) Gamma(5/8) K^(-1)/4 . J, + J = integral_0^oo omega^(-11/4) [ beta(0)^(-5/8) - Re((1 + i omega)^(-1/2) beta(omega)^(-5/8)) ] d omega. +- With `omega = u^4`, `beta = (3/16)(1 - i omega/3 + 19 omega^2/36)` and `C0 = (3/16)^(-5/8)`, we get + `J = 4 C0 integral_0^oo g2(u) du`, where `g2(u) = (1 - Re H(u^4))/u^8` and + `H(omega) = (1 + i omega)^(-1/2)(1 - i omega/3 + 19 omega^2/36)^(-5/8)`. + +**R.5 (`d = 3`).** +- `A < 0` has eigenvalues `-a1 <= -a2 < 0`. The laws of `G` and `gamma` are rotation invariant. +- On `2 x 2` symmetric matrices `dA = |lambda_1 - lambda_2| d lambda_1 d lambda_2 d theta` (the Jacobian is exactly + `|lambda_1 - lambda_2|`), with `theta` in `[0, pi)`, so `integral_(A<0) dA = pi integral_(a1>a2>0) (a1 - a2) da1 da2`. +- In the eigenframe `gamma^T adj(A) gamma = -(a2 gamma1^2 + a1 gamma2^2)`. Given `(A, b)`, + E e^(i omega Y) = e^(-i omega b a1 a2/4 - omega^2 (a1 a2)^2/12) (1 - i omega a1)^(-1/2) (1 - i omega a2)^(-1/2). +- The `b`-weight is `e^(-(5/4)(b - S/5)^2 + S^2/20 - (a1^2 + a2^2)/4)`, with `S = a1 + a2`, so `integral db` gives + `sqrt(4 pi/5) e^(S^2/20 - (a1^2 + a2^2)/4) e^(-i omega a1 a2 S/20 - omega^2 (a1 a2)^2/80)`. Hence + c1^(3) = -(192/7) 2^(1/4) 4 pi (2 pi)^(-11/2) 12^(-1/2) (1/2) sqrt(4 pi/5) pi K^(-1) T3, + Phi(omega) = (1 - i omega a1)^(-1/2) (1 - i omega a2)^(-1/2) exp(-i omega a1 a2 (a1 + a2)/20 - 23 omega^2 (a1 a2)^2/240). +- Normalize by `a1 = s^4`, `a2 = a1 v^4` and `omega = t^4/a1`. Then `da1 da2 = 16 s^7 v^3 ds dv` and + `omega^(-11/4) d omega = 4 a1^(7/4) t^(-8) dt`, and + T3 = 64 integral_0^oo ds integral_0^1 dv integral_0^oo dt s^20 v^4 (1 - v^4) e^(-s^8 E(v)) Gt, + E(v) = (2v^8 - v^4 + 2)/10 >= 3/16, Gt = (1 - Re Phi)/t^8. + +**R.6 (safe form).** +- Write `R = Re log Phi` and `I = Im log Phi`. Then `1 - Re Phi = -expm1(R) + e^R 2 sin^2(I/2)`. With + `E1(x) = expm1(x)/x`, `L(x) = log1p(x)/x`, `A(x) = atan(x)/x` and `Sinc(x) = sin(x)/x`, + Rt = -(1/4)(L(t^8) + v^8 L(t^8 v^8)) - 23 a1^2 v^8/240, It = (1/2)(A(t^4) + v^4 A(t^4 v^4)) - a1^2 v^4 (1 + v^4)/20, + Gt = -E1(t^8 Rt) Rt + e^(t^8 Rt) Sinc(t^4 It/2)^2 It^2/2. +- `g2` is the same with `R/omega^2 = -(1/4)L(omega^2) - (5/16) q L(omega^2 q)`, `q = 7/6 + 361 omega^2/1296`, and + `I/omega = -(1/2)A(omega) + (5/8)A(omega/v)/v`, `v = 3 + 19 omega^2/12`. +- No step divides by a small quantity, so the near-zero cancellation of `1 - Re Phi` is exact. +- `cusp.py` writes both integrands once, against real-interval and complex-box operations. + +## 3. Rigorous quadrature + +**Lemma Q (Gauss–Legendre on an ellipse).** Let `f` be analytic in the open Bernstein ellipse `c + h E_rho` of `[a, b]` +(`h = (b - a)/2`) with `|f| <= M` there. Then + |integral_a^b f - Q_n f| <= h . 4 M (1 + 1/(4n^2 - 1)) rho^(-2n) / (1 - rho^(-2)). +*Proof.* +- On `[-1, 1]` the Chebyshev coefficients satisfy `|a_k| <= 2 M rho^(-k)`. +- `Q_n` is exact up to degree `2n - 1`, and exact on odd `T_k` by symmetry. +- For even `k`, `|integral T_k| <= 2/(k^2 - 1)` and `|Q_n T_k| <= sum w = 2`. +- Summing over `k >= 2n` gives the bound. ∎ + +The nodes are certified in exact rational arithmetic (`gauss.py`): +- Newton steps on the grid `2^-200`, then a strict sign change of `P_n` at `r -+ 2^-150`. +- The weights `2(1 - x^2)/(n P_(n-1)(x))^2` are enclosed by an exact midpoint value plus a Lipschitz term. +- Nodes are one ulp wide and weights are tight to `2e-16` relative. + +**Lemma A (certifying analyticity).** +- The bound `M` is taken over complex boxes covering the closed ellipse: slabs along the real axis, and in `d = 3` real + sub-intervals of the other variables. +- There are two representations: + - the safe form of R.6, with the integral-hull enclosures `E1(z) = int e^(tau z)`, `L(z) = int (1 + tau z)^-1`, + `A(z) = int (1 + tau^2 z^2)^-1`, `Sinc(z) = int cos(tau z)`, or their closed forms for `|z| >= 1` (`cx.py`); + - the direct form `(1 - (Phi + Phi*)/2)/t^8`, bounded through `|z^a| = |z|^a`, where every base of a principal power must + avoid `(-oo, 0]`. Here `Phi*(w) = conj(Phi(conj w))` is the analytic coefficient-conjugate, as in `cusp.py`. It is not + a nonholomorphic `Re` of the variable, and `(Phi + Phi*)/2 = Re Phi` holds on the real axis only. +- A representation that evaluates without exception on a box is analytic on a neighbourhood of it: every denominator and + every logarithm or power argument has been checked away from its singular set. +- Every slab meets the real axis in a segment where the representation equals the integrand. By the identity theorem it is + the integrand's analytic continuation on that slab. +- The bound is therefore valid even when different slabs use different representations. + +**Lemma B (`d = 3` boxes).** +- The `(s, t)`-plane up to `s <= 3`, `t <= 15` is cut into 247 boxes. The `s`- and `v`-breakpoints include `c/t` near + the origin, because for large `t` the integrand has branch points at `|v| = 1/t` and grows for `arg s > pi/16`. +- `v` is innermost. On a box `S x T`, + `|I - Q| <= |(I_S - Q_S) I_T I_V f| + |Q_S (I_T - Q_T) I_V f| + |Q_S Q_T (I_V - Q_V) f|`. Each term takes Lemma Q in its + own variable, with the other variables real, and the positive weights summing to `|S|` and `|T|`. +- Since `f3 >= 0`, a box is enclosed by `[0, volume . sup f3]` when that is below `1e-12` (76 boxes). The other + 171 use Gauss–Legendre (1406081 nodes), each to `1e-11`. +- **Tails.** + - For `t > 15`, `integral = wgt (15^-7/7 + theta 15^-9/9)` with `|theta| <= 1`, since `|Phi| <= t^-2` on the real axis. + The `s`-integral `integral s^20 e^(-s^8 E) ds = Gamma(21/8)/(8 E^(21/8))` is exact. + - For `s > 3`, `Gt <= min(E[Y^2]/(2 a1^2), 2 t^-8)` gives less than `1e-280`. +- In `d = 2`, 15 `u`-panels up to `u = 15` are each taken to `5e-16`. The tail is + `15^-7/7 +- 1.5 . 15^-14/14`, since `|Re H(u^4)| <= (36/19)^(5/8) u^-7`. + +## 3b. The [P] torus field (Lemma S) + +**Lemma S.** Let `d in {2, 3}` and `L >= 24`, let `C_L` be the covariance of the `N` jets of (CU.2) for the [P] field in an +orthonormal frame `(u, Theta)`, and let `C_0` be that of the Gaussian kernel. Here `N = 9` in `d = 2` and `N = 14` in `d = 3`, +counting the `2d + 2` pinned coordinates `(f, d_u f, d_u^2 f, d_u^3 f, grad_Theta f, grad_Theta d_u f)` and the free +`(f4, gamma, A)`. Then `(1 - eps) C_0 <= C_L <= (1 + eps) C_0` in every frame, with `eps = N E/mu`, and +`|c1^(L)/c1^(ref) - 1| <= eta = 8 M eps`, `M = N/2 + 1`. + +*Proof.* +- **S.1 (image bound).** Every covariance entry is `(-1)^|b| D^(|a|+|b|) K(0)` contracted with frame vectors, of total + order `q <= 8` (`Cov(f4, f4) = d_u^8 K(0)`). + - For unit `v_i`, `D^q phi(x)[v]` is `phi(x)` times a sum over the partial matchings of `{1..q}` of products of + `-` and `-`. So `|D^q phi(x)[v]| <= T_8 |x|^8 phi(x)` for `|x| >= 1`, with `T_8 = 764` matchings, and + `<= 7!! = 105` at `x = 0`. + - For `d <= 3` the shell `max|n_i| = j` has at most `27 j^3` points, with `|n|^2 <= 3 j^2`. Consecutive shell terms + shrink by at least `1/2`, so `sum_(n != 0) |n|^8 e^(-288|n|^2) <= 4374 e^(-288)` and `S_24 - 1 <= 54 e^(-288)`. + - `e^(288/125) > 10` by the positive Taylor sum through 20 terms, so `e^(-288) < 10^-125`. + - `L^8 e^(-L^2 |n|^2/2)` decreases in `L` once `L|n| >= 3`, so the bound at `L = 24` holds for every `L >= 24`. + - Writing `D^q K_L(0) - D^q phi(0) = S_L^-1 sum_(n != 0) D^q phi(L n) + D^q phi(0)(S_L^-1 - 1)`, every entry of + `C_L - C_0` is at most `E = (764 . 24^8 . 4374 + 105 . 54) 10^-125 <= 3.7e-108` in absolute value. The + unnormalized periodization has only the first term. +- **S.2 (sandwich).** `C_0 - I/4` is positive definite: an exact `LDL^T` has all pivots positive, and `lambda_min(C_0)` is + about `0.258` in `d = 3`. So `|x^T (C_L - C_0) x| <= N E |x|^2 <= eps x^T C_0 x`. +- **S.3 (transfer).** Write the (CU.2) integral at direction `u` as + `I_u(C) = integral db integral dW F(W) phi_C(b, 0, W)`, with `F = |Y|^(7/4) |Delta|^(1/4) 1{A < 0} >= 0` and `2d + 1` + coordinates fixed at 0. (`pi_0 E_0[F | b]` is exactly this slice integral of the joint density.) + - `F` is homogeneous of degree `2d - 1/4` in the free jets: `Y` has degree `d` and `Delta` degree `d - 1`. + Substituting `(b, W) -> sqrt(a) (b, W)` therefore gives `I_u(a C_0) = a^(-5/8) I_u(C_0)`. + - From `(1 - eps) C_0 <= C <= (1 + eps) C_0`: `det C >= (1 - eps)^N det C_0` and `C^-1 >= C_0^-1/(1 + eps)`, and + symmetrically `det C <= (1 + eps)^N det C_0` and `C^-1 <= C_0^-1/(1 - eps)`. This gives the pointwise sandwich + + [(1 - eps)/(1 + eps)]^(N/2) phi_((1-eps) C_0) <= phi_C <= [(1 + eps)/(1 - eps)]^(N/2) phi_((1+eps) C_0). + + - Integrating the nonnegative `F` gives `I_u(C_L) / I_u(C_0)` in `[rho^-M, rho^M]`, with `rho = (1 + eps)/(1 - eps)`. + - `log rho <= 4 eps` and `e^x <= 1 + 2x` on `[0, 1]`, so this lies within `8 M eps` of 1. + - `I_u(C_0)` does not depend on `u` (isotropy), so integrating over `S^(d-1)` preserves the bound. ∎ + +The constants, as exact rationals in `RESULTS.json` (`rules`), are as follows: + +| `d` | `N` | `mu` | `eps` | `eta` | +|---|---|---|---|---| +| `2` | `9` | `1/4` | `<= 1.4e-106` | `<= 5.9e-105` | +| `3` | `14` | `1/4` | `<= 2.1e-106` | `<= 1.4e-104` | + +In binary64 the factor `[1 - eta, 1 + eta]` is enclosed by the one-ulp interval around 1. + +**Cross-check of Lemma R.1.** The self-test builds `C_0` exactly and regresses the free jets on the pins in rational +arithmetic. The result must reproduce R.1 exactly: +- `E[f4 | b] = -3b`, `Var = 24`; +- `gamma ~ N(0, 2 I)`; +- `E[A_jj | b] = -b`, `Var A_jj = 2`, `Var A_jk = 1`; +- all cross-covariances 0; +- `det C_P = 12` and `(C_P^-1)_ff = 3/2`, which is `pi_0 = (2 pi)^(-(d+1)) 12^(-1/2) e^(-3b^2/4)`. + +## 4. Arithmetic + +- **`ia.py`** is byte-identical to Math-#217's (`frontiers/c8_elder_failure_coefficient_20261001/ia.py`). It provides + binary64 with one-ulp outward rounding, `pi` and `ln 2` from exact rational series, and `exp(-y)` with an a priori Taylor + bound. +- **`elem.py`** adds `log`/`log1p` (`2 atanh` series), `atan` (two halvings plus the alternating series), `sin`/`cos` + (reduction mod `pi/2`, alternating series; interval ranges include `+-1` at enclosed critical points), and the four + kernels with their monotonicity. +- **`consts.py`** computes `log Gamma` by Stirling's series at `q + 10`, with the remainder bounded by the first omitted + term, together with the exact rational product `prod (q + j)`. +- No libm transcendental is trusted. +- **Results are rounded outward.** + - Both representations of each bound are reduced to a single float. + - Every sum is an interval sum. + - `RESULTS.json` widens by `1e-12` relative and rounds outward to 16 digits. +- `--check` recomputes everything. It requires each recomputed interval to lie inside the published one, and the published + one to be at most `3e-12` relatively looser. It also requires identical rule data (panels, box counts, node counts). + +## 5. Parts + +| quantity | enclosure | +|---|---| +| `integral_0^oo g2(u) du` (`d = 2`) | `[0.70583309158145, 0.70583309158290]` | +| `P2` (so that `c1^(2) = P2 . integral g2`) | `[-0.3816749722932, -0.3816749722921]` | +| `T3` (`d = 3`) | `[45.039309847, 45.039309924]` | +| `T3`: boxes / `t > 15` tail | `[45.039309787, 45.039309864]` / `[5.948E-8, 5.995E-8]` | +| `P3` (so that `c1^(3) = P3 . T3`) | `[-0.004703632154886, -0.004703632154874]` | +| `K = pi/(2 Gamma(11/4) sin(7 pi/8))` | `[2.552096599674, 2.552096599680]` | +| `Gamma(5/8)`, `Gamma(11/4)` | `[1.434518848089, 1.434518848093]`, `[1.608359421983, 1.608359421988]` | + +`d = 2` panel rules, as `[a, b, n, rho]`: [0, 0.4, 13, 4], [0.4, 0.7, 13, 4], [0.7, 0.9, 14, 4], [0.9, 1, 13, 4], [1, 1.1, 13, 4], [1.1, 1.3, 13, 4], [1.3, 1.6, 16, 3], [1.6, 2, 15, 3], [2, 2.6, 15, 3], [2.6, 3.4, 14, 3], [3.4, 4.5, 13, 3], [4.5, 6, 15, 2.4], [6, 8, 14, 2.4], [8, 11, 13, 2.4], [11, 15, 11, 2.4]. + +## 6. Verification + +`python3 -B -S certificate.py --check --procs N` runs a self-test, recomputes everything (about 4 minutes on four cores) +and compares the result with `RESULTS.json`. + +The self-test checks: +1. Gauss–Legendre exactness on even monomials for `n = 4 ... 48`. +2. The real-interval safe forms against the direct complex-power formulas, in floating point, at ten points. +3. A deliberately coarse 6-point rule on `[0.85, 1.15]`, whose certified enclosure must contain a 40-point floating + reference. This checks the error bound itself. +4. The `t`-tail against its floating formula. +5. Lemma S: the exact regression of R.1 (§3b), the image-bound order against the largest covariance order, `T_8 = 764`, + the positive definiteness of `C_0 - I/4`, and the transcription of `c_(d,24)` against the repository copy of + `ENCLOSURE.json` (when reachable). + +The assembly also compares `integral g2` and `T3` with coarse floating quadratures of the safe form in libm floating +point. This guards normalizations and factors. + +Six seeded defects must each be rejected by the self-test (exit 1): +- `kernel-A`: a `1e-7` change in `A`; +- `gl-weight`: weights `x (1 + 1e-9)`; +- `bound-scale`: error bounds `x 1e-6`; +- `phase-sign`: the sign of the `a1^2` phase term in `f3`; +- `tail-drop`; +- `image-order6`: an image bound valid only to order 6, as in side24_v1 §2, which is too low for `Cov(f4, f4)`. + +The workflow `.github/workflows/cusp-coefficient-certified.yml` runs, for `-B -S` and `-B -O -S`: the manifest and pin +check, `--check`, the six mutants, and a clean-tree check. + +## 7. Floating controls (`controls.py`; not part of the certificate) + +These recompute the same quantities in libm floating point by routes that share only the Gauss–Legendre nodes +(`gauss.py`) with the certificate. They guard against a wrong reduction, not against rounding. + +1. **`d = 2` without the Mellin identity.** Given `(a = -A, gamma, b)`, `Y` is Gaussian with standard deviation `a/sqrt 6` + and mean `(ab - gamma^2)/4`, so `E|Y|^(7/4)` is a Kummer function. Averaging over `b` leaves a 2D integral over + `(a, gamma)`, done by a product Gauss rule. Result: `-0.269398825672854`, `3.5e-15` from the certified midpoint. +2. **`d = 3` by the same route.** The eigenvalues `a1 > a2` of `-A`, the Vandermonde factor and `gamma ~ N(0, 2 I)` give a + 4D integral. Result: `-0.211848346210999`, `3.8e-13` from the certified midpoint. + + Both controls lie inside the certified intervals. +3. **The sources' quoted values.** Each is the correct rounding of the certified midpoint: the difference is below half a + unit in its last digit. + +| quoted value | quoted | certified midpoint | difference | half-unit | +|---|---|---|---|---| +| [CU] §8: `c1`, `d = 2` | `-0.26939883` | `-0.2693988257` | `-4.3e-09` | `5e-09` | +| [CU] §8: `c1`, `d = 3` | `-0.21184835` | `-0.2118483462` | `-3.8e-09` | `5e-09` | +| [CU] §8: `c`, `d = 2` | `0.0734069193` | `0.0734069193` | `-6.0e-12` | `5e-11` | +| [CU] §8: `c`, `d = 3` | `0.0417759318` | `0.0417759318` | `-4.1e-11` | `5e-11` | +| [CU] §8: `c1/c`, `d = 2` | `-3.669938` | `-3.6699377691` | `-2.3e-07` | `5e-07` | +| [CU] §8: `c1/c`, `d = 3` | `-5.071062` | `-5.0710621374` | `1.4e-07` | `5e-07` | +| [CU] §8: `I^cand`, `d = 2` | `-0.1772744` | `-0.1772743968` | `-3.2e-09` | `5e-08` | +| [CU] §8: `I^cand - c1`, `d = 2` | `0.0921244` | `0.0921244289` | `-2.9e-08` | `5e-08` | +| [CU] §8: `I^cand`, `d = 3` | `-0.1394041` | `-0.1394040516` | `-4.8e-08` | `5e-08` | +| [CU] §8: `I^cand - c1`, `d = 3` | `0.0724443` | `0.0724442946` | `5.4e-09` | `5e-08` | +| Math-#216 table: `c1`, `d = 1` (Math-#214 `C_1`) | `-0.227606` | `-0.2276063588` | `3.6e-07` | `5e-07` | + +[CU]'s three `d = 3` Monte Carlo runs (§8, remark 3), against the certified value: +- shipped fixed-seed: `-0.211630 +- 0.000140` (`+1.56` standard errors); +- development (`2e8` samples): `-0.211880 +- 0.000100` (`-0.32` standard errors); +- referee (`5e8` samples): `-0.211848 +- 0.000024` (`+0.01` standard errors). + +4. **Lemma S, actual deviations.** A direct floating lattice sum over `|n_i| <= 2` computes the deviations + `D^g K_24(0) - D^g phi(0)` of all jet covariance entries (`d = 3`, standard frame, variance-one normalization). The + largest is `1.756e-114`, at `g = (8, 0, 0)`, the `d_u^8` entry `Cov(f4, f4)`. The certified entry bound is + `E = 3.678e-108`, which exceeds it by a factor of about `2e6`. + +## 8. What this does not do + +- **Sides `L >= 24` only.** + - For `L < 24` the premises of S.1 are not established, and no torus statement is made. + - `d = 1` is the Gaussian kernel only. + - The ratios `c1/c_(d,24)` consume side24_v1's enclosures of `c_(d,24)`. These are author-side (OpenAI), with + nonauthor review open; they are not reviewed here. +- **Not a review of Theorem CU.** The numbers are the integral (CU.2) defines. That `c1` is the coefficient of + `ell^(1/4)` in `nu_eld` is [CU]'s statement (author-side; it depends on Math-#191 and Math-#198), consumed, not + reviewed. +- No `c2` (Math-#216), no `d >= 4`, no rate. +- **The certificate is bound to its pinned trajectory, not a reusable interval library** (C54 reviews 5384790625 and + 5384863435). Two out-of-trajectory limits are recorded: + - `d3._sub` does not always keep its exact input endpoint (C54-A-BOUNDARY-01). For example, `(a, b, k) = + (0.024016010673783495, 0.056037358238828154, 2)` loses `2^-57`. The published run never reaches such a case: the + reviewers' exhaustive actual-call audit and reachability trace found no gap. + - `ia.mul((0, 0), (inf, inf))` returns NaN rather than raising (C54-A-FINITE-02). Every bound the published run + returns is finite and nonnegative. + + Any change to grids, `rho`, cutoffs or the search must pin `_sub`'s endpoints, add finite-result guards and renew the + coverage checks. +- Same GitHub account as every lane; zero organizational-independence credit. **I will not merge.** + +## 9. Provenance + +- **Sources** (`SOURCE_MAP.json`): + - [CU] Math-#207 `frontiers/cusp_second_order_20261001/PROOF.md` (blob `f6df5a73`). It is on `main` since `566b1a1` + (byte-identical to the consumed head `826b8be`) and is now checked on this tree (`current_required`); + - [SIDE24] `coefficients/side24_v1/PROOF.md` (1), and §§2, 4 as the model for Lemma S; + - [SIDE24-E] `coefficients/side24_v1/ENCLOSURE.json` (`c_(d,24)`); + - [LP] `imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md` §15. + + The last three are verified on this tree by the workflow. +- **Files:** + - the certificate: `certificate.py` (drivers, self-test, mutants, publication), `cusp.py` (integrands), `d3.py` + (`d = 3` boxes and tails), `side24.py` (Lemma S, exact rationals), `gauss.py`, `cx.py`, `elem.py`, `ia.py`, `consts.py`; + - `controls.py`; + - `RESULTS.json`, `SOURCE_MAP.json`, `SOURCE_FILES.json` (the manifest). diff --git a/frontiers/cusp_coefficient_certified_20261001/README.md b/frontiers/cusp_coefficient_certified_20261001/README.md new file mode 100644 index 00000000..8870c5bf --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/README.md @@ -0,0 +1,15 @@ +# Certified cusp coefficient `c1` (Math-#207, (CU.2)), Gaussian kernel + +`CL-CU-CUSP-COEFFICIENT-20261001-v1` (v1.1). Author-side certified numerics; scientific effect NONE. See `NOTE.md`. + +| `d` | `c1` | `c1 / c_(d,ref)` | +|---|---|---| +| `1` | `[-0.22760635877558, -0.22760635877504]` | | +| `2` | `[-0.269398825674, -0.269398825672]` | `[-3.66993776909, -3.66993776907]` | +| `3` | `[-0.211848347, -0.211848346]` | `[-5.07106215, -5.07106213]` | + +For the [P] torus field of every side `L >= 24` (Lemma S), including SIDE24: `c1/c_(3,24)` is in +`[-5.07106215, -5.07106213]` and `c1/c_(2,24)` in `[-3.66993776909, -3.66993776907]`. + +Replay: `python3 -B -S certificate.py --check --procs N` (self-test, then a full recomputation against `RESULTS.json`). +Standard library only. diff --git a/frontiers/cusp_coefficient_certified_20261001/RESULTS.json b/frontiers/cusp_coefficient_certified_20261001/RESULTS.json new file mode 100644 index 00000000..8bbd0b0c --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/RESULTS.json @@ -0,0 +1,285 @@ +{ + "object": "CL-CU-CUSP-COEFFICIENT-20261001-v1", + "scientific_effect": "NONE", + "statement": "interval enclosures (outward decimals) of the cusp coefficient c1 of Math-#207 (CU.2) for the Gaussian kernel, d = 1, 2, 3; see NOTE.md", + "values": { + "d=1": { + "c1": [ + "-0.2276063587755705", + "-0.2276063587750459" + ], + "I_cand = (3^(1/4)/2) c1": [ + "-0.1497734069836378", + "-0.1497734069832868" + ], + "I_cand - c1": [ + "0.07783295179170776", + "0.07783295179198402" + ] + }, + "d=2": { + "c1": [ + "-0.2693988256731888", + "-0.2693988256725268" + ], + "I_cand = (3^(1/4)/2) c1": [ + "-0.1772743967942046", + "-0.1772743967937623" + ], + "I_cand - c1": [ + "0.09212442887867683", + "0.09212442887907186" + ], + "c_(d,ref)": [ + "0.07340691930595442", + "0.07340691930611396" + ], + "c1 / c_(d,ref)": [ + "-3.669937769082580", + "-3.669937769072928" + ] + }, + "d=3": { + "c1": [ + "-0.2118483463910089", + "-0.2118483460317414" + ], + "I_cand = (3^(1/4)/2) c1": [ + "-0.1394040516860850", + "-0.1394040514496684" + ], + "I_cand - c1": [ + "0.07244429434593535", + "0.07244429494106150" + ], + "c_(d,ref)": [ + "0.04177593184055288", + "0.04177593184064373" + ], + "c1 / c_(d,ref)": [ + "-5.071062141698316", + "-5.071062133097566" + ] + }, + "[P] field, every L >= 24 (Lemma S)": { + "d=2: c1": [ + "-0.2693988256731890", + "-0.2693988256725267" + ], + "d=2: I_cand": [ + "-0.1772743967942047", + "-0.1772743967937622" + ], + "d=2: I_cand - c1": [ + "0.09212442887867663", + "0.09212442887907205" + ], + "d=2: c_(d,24) (side24_v1, consumed)": [ + "0.07340691930596084", + "0.07340691930610771" + ], + "d=2: c1 / c_(d,24), L = 24": [ + "-3.669937769082261", + "-3.669937769073238" + ], + "d=3: c1": [ + "-0.2118483463910090", + "-0.2118483460317414" + ], + "d=3: I_cand": [ + "-0.1394040516860851", + "-0.1394040514496684" + ], + "d=3: I_cand - c1": [ + "0.07244429434593524", + "0.07244429494106161" + ], + "d=3: c_(d,24) (side24_v1, consumed)": [ + "0.04177593184055655", + "0.04177593184064014" + ], + "d=3: c1 / c_(d,24), L = 24": [ + "-5.071062141697873", + "-5.071062133098000" + ] + }, + "parts": { + "d=2: int_0^oo g2(u) du": [ + "0.7058330915814585", + "0.7058330915829000" + ], + "d=2: tail u > 15": [ + "8.361094744137853E-10", + "8.361094817557882E-10" + ], + "d=2: P2": [ + "-0.3816749722931025", + "-0.3816749722921806" + ], + "d=3: T3": [ + "45.03930984738922", + "45.03930992375617" + ], + "d=3: boxes": [ + "45.03930978790707", + "45.03930986381351" + ], + "d=3: tail t > 15": [ + "5.948215822041717E-8", + "5.994264845299655E-8" + ], + "d=3: tail s > 3": [ + "-4.940656458412466E-324", + "6.400000000006401E-279" + ], + "d=3: P3": [ + "-0.004703632154885225", + "-0.004703632154874365" + ], + "K = pi/(2 Gamma(11/4) sin(7 pi/8))": [ + "2.552096599674470", + "2.552096599679964" + ], + "Gamma(5/8)": [ + "1.434518848089029", + "1.434518848092093" + ], + "Gamma(11/4)": [ + "1.608359421983816", + "1.608359421987275" + ], + "Gamma(11/8)": [ + "0.8889135691552759", + "0.8889135691571741" + ], + "Gamma(7/6)": [ + "0.9277193336290466", + "0.9277193336310296" + ] + } + }, + "rules": { + "d=2 panels [a, b, n, rho]": [ + [ + 0.0, + 0.4, + 13, + 4.0 + ], + [ + 0.4, + 0.7, + 13, + 4.0 + ], + [ + 0.7, + 0.9, + 14, + 4.0 + ], + [ + 0.9, + 1.0, + 13, + 4.0 + ], + [ + 1.0, + 1.1, + 13, + 4.0 + ], + [ + 1.1, + 1.3, + 13, + 4.0 + ], + [ + 1.3, + 1.6, + 16, + 3.0 + ], + [ + 1.6, + 2.0, + 15, + 3.0 + ], + [ + 2.0, + 2.6, + 15, + 3.0 + ], + [ + 2.6, + 3.4, + 14, + 3.0 + ], + [ + 3.4, + 4.5, + 13, + 3.0 + ], + [ + 4.5, + 6.0, + 15, + 2.4 + ], + [ + 6.0, + 8.0, + 14, + 2.4 + ], + [ + 8.0, + 11.0, + 13, + 2.4 + ], + [ + 11.0, + 15.0, + 11, + 2.4 + ] + ], + "d=3": { + "boxes": 247, + "crude": 76, + "gauss_legendre": 171, + "nodes": 1406081, + "T_MAX": 15.0, + "S_MAX": 3.0 + }, + "Lemma S (exact rationals)": { + "d=2": { + "N": "9", + "E": "183921320046627603/50000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", + "mu": "1/4", + "eps": "1655291880419648427/12500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", + "M": "11/2", + "eta": "18208210684616132697/3125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000" + }, + "d=3": { + "N": "14", + "E": "183921320046627603/50000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", + "mu": "1/4", + "eps": "1287449240326393221/6250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", + "M": "8", + "eta": "1287449240326393221/97656250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000" + } + } + }, + "publication": { + "margin_relative": 1e-12, + "digits": 16 + } +} diff --git a/frontiers/cusp_coefficient_certified_20261001/SOURCE_FILES.json b/frontiers/cusp_coefficient_certified_20261001/SOURCE_FILES.json new file mode 100644 index 00000000..8a329eb9 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/SOURCE_FILES.json @@ -0,0 +1,89 @@ +{ + "object": "CL-CU-CUSP-COEFFICIENT-20261001-v1", + "files": [ + { + "path": "NOTE.md", + "bytes": 27264, + "sha256": "9c74117f6b7ff2054fcfb5ddcfb4c8ad4a3a1278c201faafc5e98e7612abfc0c", + "git_blob": "b2866871c10cd0a0def41e10bdaf5ecbef18241c" + }, + { + "path": "README.md", + "bytes": 774, + "sha256": "58c23e1062ee0d8d8d3f17af4dd6bfca76fc897ab51fc1dcc08481e4bf5b3c39", + "git_blob": "8870c5bf794baa0e9d1d2eca4f1ac2a509cf30e3" + }, + { + "path": "RESULTS.json", + "bytes": 5350, + "sha256": "2de0a2381eb2f73b3f2ed7367c78d28856067b271f9e3484384c199a24cacdbe", + "git_blob": "8bbd0b0c69f6b29849fb2e3a3b2ad86974130987" + }, + { + "path": "SOURCE_MAP.json", + "bytes": 3336, + "sha256": "306b39a011472e6c86642c1226f8605f4b4ec9b01308e0c7b0233c24023f1390", + "git_blob": "7e9092f185c55efe22bc9587be5cdf2f5ef89bd1" + }, + { + "path": "certificate.py", + "bytes": 18392, + "sha256": "7472ac07d96210801b4070771713b0732bafb2cc479dd74d733eebc097b71838", + "git_blob": "7a4ba3e80acb3bb2423429e4808138574303b82d" + }, + { + "path": "consts.py", + "bytes": 4462, + "sha256": "5265668c0f489137d6362d230d7f2df0dfb25632abdadc0934a2c89c8594759e", + "git_blob": "5cac30c3bf086697b978f1161191764121298fd4" + }, + { + "path": "controls.py", + "bytes": 8690, + "sha256": "e2d59dfa162947e02d3ba696e5fb89aa926d2a658d4252f0ec1d0e6806ec3ddc", + "git_blob": "5e41d95cf30ab47581ad15e78577a80a9d7e5384" + }, + { + "path": "cusp.py", + "bytes": 7686, + "sha256": "695038fb8de5c7f13fd481384a29117df297c4d7125b7b761ee4114e7d0f19fc", + "git_blob": "fba24e5ee979ac23c6c00935af6ba1d1ecb09590" + }, + { + "path": "cx.py", + "bytes": 4716, + "sha256": "af2a32b3cf6b3826520ac3a14bc4238143447c5e3697ad0355d8a55e6b9aacd1", + "git_blob": "29e3a4a8a46f065346592320e0c6bb65c6ccb331" + }, + { + "path": "d3.py", + "bytes": 11136, + "sha256": "eafa1a6cf76cbdd9338f315377522aa3c23b361216c77bcd5e27ef610f57b0f7", + "git_blob": "bae664536aac0f673cd8cd211d8637b6ba34c50c" + }, + { + "path": "elem.py", + "bytes": 9825, + "sha256": "65490f314fd731ced8d893d137ed5d7a516cf3cf36ba68cf77d7845805c55fb0", + "git_blob": "99ad3332b40daa1eee50106dbaa8b4b726cd5ae6" + }, + { + "path": "gauss.py", + "bytes": 6228, + "sha256": "51cd9e56f1858f711bee0a4554727d847f22a0796e85c77be58b4e74f2278eb0", + "git_blob": "71f0b5709d5867497a5a1da1c08de48c30d68976" + }, + { + "path": "ia.py", + "bytes": 8930, + "sha256": "6d048d58c6ab18f21c0609fc888ca73497a60e234217bb4f7a309a1502fccc9e", + "git_blob": "cf4bbfa16f577fe65f4ab5fb3621943fa63cf9fa" + }, + { + "path": "side24.py", + "bytes": 10043, + "sha256": "85c330d82623dc89effed5a09b8d15a8986688f5eacc922f8b9b87754f9eef35", + "git_blob": "31da88a031f3df682d1411415bbc1a9a7f591850" + } + ] +} diff --git a/frontiers/cusp_coefficient_certified_20261001/SOURCE_MAP.json b/frontiers/cusp_coefficient_certified_20261001/SOURCE_MAP.json new file mode 100644 index 00000000..7e9092f1 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/SOURCE_MAP.json @@ -0,0 +1,47 @@ +{ + "object": "CL-CU-CUSP-COEFFICIENT-20261001-v1", + "delivered_under": "explicit owner instructions of Dylan Roy in Claude Code session session_015wNj8LPTKXsaT68G3DgPPh, given after OWNER_STOP (2026-09-27): 2026-09-28 \"you have my full permisions and approval to make the proper changes and uploads\"; 2026-09-30 \"Work on closing as many mathematical items as possible while working with and supporting the other models in the GitHub. Also advance the mathematics in every place that's justified based on what we currently know. Ultra tokens, ultra depth @GitHub @Superpowers\" (repeated 2026-09-30 as \"Continue to Work on closing ...\"). Pull requests only, never merged by their author; no register, GRAPH, STATUS, lemma_closed, prize or premise edits. OWNER_STOP.md addresses Cursor agents and automations and is untouched.", + "sources": [ + { + "label": "[CU]", + "path": "frontiers/cusp_second_order_20261001/PROOF.md", + "ref": "Math-#207, merged at 566b1a1 (byte-identical to the consumed head 826b8be)", + "consumed": "Theorem CU, (CU.2) and its parity form: the definition of c1 (jets, Y, Delta, E_0, v_0(b,0), pi_0, the sphere integral); section 8 remark 3 floating values (controls only); for Lemma S: section 0 (Theorem CU for the [P] field on R^d/(L Z^d), every L > 0) and section 8 remark 5 (only the law of (f4, gamma, A), hence c1, depends on the field)", + "current_required": true, + "bytes": 52605, + "sha256": "cffc1b83ac364da58bd650c3c301ea1babd38e19af54f74010bdc120edce38a8", + "git_blob": "f6df5a7356a54771b0052e42a864a140709903ad" + }, + { + "label": "[SIDE24]", + "path": "coefficients/side24_v1/PROOF.md", + "ref": "main 3e0a91b", + "consumed": "equation (1): c_(d,ref) = Gamma(7/6)(3/2)^(1/3) D_(d-1)/(2 sqrt3 pi^(d-1) sqrt(pi)), D_1 = 4/3, D_2 = 29/6 - sqrt6; section 2 (the periodic-image bound, extended here to order 8 in Lemma S) and section 4 (the density sandwich), as a model for Lemma S", + "current_required": true, + "bytes": 10272, + "sha256": "c06daccc4ba4b9168522b9888b76a7d599934fc3b91bd753ee5d492262917769", + "git_blob": "44b66f04f89fcd87383b3603fa69f1feb64cdddd" + }, + { + "label": "[SIDE24-E]", + "path": "coefficients/side24_v1/ENCLOSURE.json", + "ref": "main 3e0a91b (last change d83af85)", + "consumed": "dimensions: the outward decimal enclosures of c_(2,24) and c_(3,24), transcribed in side24.py C_24 (checked against this file by the self-test when the repository copy is reachable); used only for the ratios c1/c_(d,24)", + "current_required": true, + "bytes": 1090, + "sha256": "72b6cd92d31394cdaf5da8919a5d548e902228af1f095cc184158a71d8287811", + "git_blob": "57af39a05e14ed0ba8ebc00a9b4aca4dffb067c7" + }, + { + "label": "[LP]", + "path": "imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md", + "ref": "main 3e0a91b", + "consumed": "section 15: parity factorization pi_0 = p_(f,V_u)(b,0) p_G(0) phi_tau(12k), tau^2 = Var(t_u | G = 0)", + "current_required": true, + "bytes": 40261, + "sha256": "9350ad6eaba6626b93c3dedeef9e2ff816e5cdf1c8318e85fb27499141c84bc7", + "git_blob": "dfed3b8d318a3ab1950957f393307733a4bef3f2" + } + ], + "note": "current_required entries, including [CU] (Math-#207, merged at 566b1a1, byte-identical to the consumed head 826b8be), are verified byte-for-byte on the checked-out tree by the workflow." +} diff --git a/frontiers/cusp_coefficient_certified_20261001/certificate.py b/frontiers/cusp_coefficient_certified_20261001/certificate.py new file mode 100644 index 00000000..7a4ba3e8 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/certificate.py @@ -0,0 +1,479 @@ +"""Certificate for CL-CU-CUSP-COEFFICIENT-20261001: rigorous enclosures of the cusp coefficient c1 of Theorem CU +(Math-#207, (CU.2)) for the Gaussian kernel in d = 1, 2, 3, of I^cand = (3^(1/4)/2) c1, and of the ratios c1/c_(d,ref); +by Lemma S (side24.py) also for the [P] torus field of every side L >= 24 in d = 2, 3, with the ratios c1/c_(d,24). + + python3 -B -S certificate.py --write [--procs N] compute and write RESULTS.json + python3 -B -S certificate.py --check [--procs N] self-test, recompute, compare with RESULTS.json + python3 -B -S certificate.py --check --mutant NAME a seeded defect; the self-test must reject it (exit 1) + +Standard library only. All arithmetic is outward-rounded binary64 interval arithmetic (ia.py, elem.py, cx.py). +NOTE.md states the reduction and the error bounds.""" +import argparse +import json +import math +import os +import sys +import time +from fractions import Fraction as Fr +from multiprocessing import Pool + +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +sys.dont_write_bytecode = True + +import ia +from ia import dn, up, add, sub, mul, div, scal, from_frac +import elem as EL +import gauss as G +import cusp as CU +import cx as CX +import consts as K +import d3 +import side24 as S24 + +PACKET = 'CL-CU-CUSP-COEFFICIENT-20261001-v1' +U_EDGES = (0.0, 0.4, 0.7, 0.9, 1.0, 1.1, 1.3, 1.6, 2.0, 2.6, 3.4, 4.5, 6.0, 8.0, 11.0, 15.0) +U_MAX = U_EDGES[-1] +TOL2 = 5e-16 # per u-panel, d = 2 +TOL3 = 1e-11 # per (s, t) box, d = 3 +RHOS2 = (4.0, 3.0, 2.4, 2.0, 1.7, 1.5, 1.35, 1.25, 1.18, 1.12) +MARGIN = 1e-12 # relative widening before publication +DIGITS = 16 # significant digits of the published decimals +MUTANTS = ('kernel-A', 'gl-weight', 'bound-scale', 'phase-sign', 'tail-drop', 'image-order6') +MUTANT = None + + +# ---------------------------------------------------------------- mutants (seeded defects for the workflow) +def install_mutant(name): + global MUTANT + MUTANT = name + if name == 'kernel-A': + orig = EL.A_pt + + def A_pt(t): + return orig(t * 1.0000001) + EL.A_pt = A_pt + elif name == 'gl-weight': + orig_rule = G.rule + + def rule(n): + xs, ws = orig_rule(n) + return xs, [scal(1.0 + 1e-9, w) for w in ws] + G.rule = rule + elif name == 'bound-scale': + orig_eb = G.error_bound + + def error_bound(n, rho, M, h): + return orig_eb(n, rho, M, h) * 1e-6 + G.error_bound = error_bound + elif name == 'phase-sign': + CU.MUTANT = 'phase-sign' + elif name == 'tail-drop': + d3.tail_t = lambda: (0.0, 0.0) + elif name == 'image-order6': + S24.IMAGE_ORDER = 6 + elif name is not None: + raise SystemExit('unknown mutant: %s' % name) + + +# ---------------------------------------------------------------- independent floating references (self-test only) +def _g2_float(u): + w = u ** 4 + H = (1 + 1j * w) ** -0.5 * (1 - 1j * w / 3 + 19 * w * w / 36) ** -0.625 + return (1 - H.real) / (w * w) + + +def _f3_float(s, v, t): + a1 = s ** 4 + a2 = a1 * v ** 4 + w = t ** 4 / a1 + Phi = ((1 - 1j * w * a1) ** -0.5 * (1 - 1j * w * a2) ** -0.5 + * complex(math.e) ** (-1j * w * a1 * a2 * (a1 + a2) / 20 - 23 * w * w * (a1 * a2) ** 2 / 240)) + Gt = (1 - Phi.real) / (w * a1) ** 2 + return s ** 20 * v ** 4 * (1 - v ** 4) * math.exp(-s ** 8 * (2 * v ** 8 - v ** 4 + 2) / 10) * Gt + + +def _E1f(x): + return math.expm1(x) / x if x != 0.0 else 1.0 + + +def _Lf(x): + return math.log1p(x) / x if x != 0.0 else 1.0 + + +def _Af(x): + return math.atan(x) / x if x != 0.0 else 1.0 + + +def _Sf(x): + return math.sin(x) / x if x != 0.0 else 1.0 + + +def _Gf(Rt, It, R, I): + return -_E1f(R) * Rt + math.exp(R) * _Sf(I / 2) ** 2 * It * It / 2 + + +def _g2_safe_float(u): + """g2 by the safe formulation in plain floating point (libm), for the assembly check.""" + w = u ** 4 + w2 = w * w + q = 7 / 6 + 361 * w2 / 1296 + Rt = -0.25 * _Lf(w2) - 5 / 16 * q * _Lf(w2 * q) + v = 3 + 19 * w2 / 12 + It = -0.5 * _Af(w) + 5 / 8 * _Af(w / v) / v + return _Gf(Rt, It, w2 * Rt, w * It) + + +def _f3_safe_float(s, v, t): + a1sq = s ** 8 + v4, t4 = v ** 4, t ** 4 + v8, t8 = v4 * v4, t4 * t4 + Rt = -0.25 * (_Lf(t8) + v8 * _Lf(t8 * v8)) - 23 / 240 * a1sq * v8 + It = 0.5 * (_Af(t4) + v4 * _Af(t4 * v4)) - a1sq * v4 * (1 + v4) / 20 + G = _Gf(Rt, It, t8 * Rt, t4 * It) + return s ** 20 * v4 * (1 - v4) * math.exp(-a1sq * (2 * v8 - v4 + 2) / 10) * G + + +def _gl_float(f, a, b, n=40): + xs, ws = G.rule(n) + c, h = 0.5 * (a + b), 0.5 * (b - a) + return sum(0.5 * (w[0] + w[1]) * h * f(c + h * 0.5 * (x[0] + x[1])) for x, w in zip(xs, ws)) + + +def selftest(): + """Returns a list of failure strings (empty when all checks pass).""" + fails = [] + # (1) Gauss-Legendre exactness on even monomials, for a sample of orders: a defect detector, not the proof of the rules. + for n in (4, 6, 8, 12, 16, 20, 24, 32, 40, 48): + xs, ws = G.rule(n) + for k in range(0, n, max(1, n // 6)): + s = (0.0, 0.0) + for x, w in zip(xs, ws): + s = add(s, mul(w, ia.pw(x, 2 * k))) + if not (Fr(s[0]) <= Fr(2, 2 * k + 1) <= Fr(s[1])): + fails.append('GL n=%d exactness at degree %d' % (n, 2 * k)) + break + # (2) the safe-form integrands against the direct complex-power formulas (floating), at moderate arguments + for u in (0.3, 0.7, 1.0, 1.4, 2.2): + iv = CU.g2(CU.RealOps, (u, u)) + ref = _g2_float(u) + if not abs(0.5 * (iv[0] + iv[1]) - ref) <= 1e-10 * abs(ref): + fails.append('g2(%g): %r vs %r' % (u, iv, ref)) + for s, v, t in ((1.2, 0.6, 0.5), (1.4, 0.9, 1.2), (0.9, 0.3, 2.0), (1.6, 0.8, 0.8), (1.3, 0.5, 3.0)): + iv = CU.f3(CU.RealOps, (s, s), (v, v), (t, t)) + ref = _f3_float(s, v, t) + if not abs(0.5 * (iv[0] + iv[1]) - ref) <= 1e-9 * abs(ref): + fails.append('f3(%g, %g, %g): %r vs %r' % (s, v, t, iv, ref)) + # (3) a deliberately coarse rule: the certified enclosure must contain a 40-point floating reference. + iv, n, rho, M, e = panel2(0.85, 1.15, 1e-9, nmin=6) + ref = _gl_float(_g2_float, 0.85, 1.15) + if not (iv[0] <= ref <= iv[1]): + fails.append('coarse d = 2 panel: %r does not contain %r' % (iv, ref)) + # (4) the t-tail against its floating formula + T = d3.T_MAX + tt = d3.tail_t() + W = math.gamma(21 / 8) / 8 * _gl_float(lambda v: v ** 4 * (1 - v ** 4) * ((2 * v ** 8 - v ** 4 + 2) / 10) ** (-21 / 8), 0.0, 1.0) + ref = W * T ** -7 / 7 + if not (tt[0] <= ref <= tt[1]) or not tt[1] - tt[0] < 1e-2 * ref: + fails.append('t-tail: %r vs %r' % (tt, ref)) + # (5) Lemma S: the exact jet covariance C_0 reproduces Lemma R.1 (regression, det C_P = 12, (C_P^-1)_ff = 3/2), and + # the premises of the transfer hold (image bound order, matchings count, positive definiteness of C_0 - mu I) + for d in (2, 3): + fails.extend(S24.lemma_r_check(d)) + try: + S24.transfer(d) + except ValueError as e: + fails.append('Lemma S, d = %d: %s' % (d, e)) + if S24.matchings(8) != 764: + fails.append('T_8 != 764') + enc = S24.enclosure_check(HERE) + if enc: + fails.extend(enc) + return fails + + +# ---------------------------------------------------------------- assembly cross-checks (floating, coarse) +def assembly_check(name, iv, ref, rel): + """The certified enclosure must agree with an independent coarse floating quadrature of the same integral to the + relative tolerance rel. This guards the assembly (normalizations, factors, tails) and is not part of the bound.""" + mid = 0.5 * (iv[0] + iv[1]) + if not abs(mid - ref) <= rel * abs(ref): + raise RuntimeError('assembly check failed for %s: certified mid %r vs floating %r' % (name, mid, ref)) + + +def float_int_g2(): + """int_0^oo g2 by floating Gauss-Legendre (40 nodes per u-panel, safe formulation in libm floating point) plus the + main tail term.""" + return sum(_gl_float(_g2_safe_float, a, b) for a, b in zip(U_EDGES, U_EDGES[1:])) + U_MAX ** -7 / 7 + + +def float_T3(): + """T3 = 64 int int int f3 by coarse floating Gauss-Legendre with the safe formulation in libm floating point.""" + tot = 0.0 + xs, ws = G.rule(24) + def nodes(a, b): + c, h = 0.5 * (a + b), 0.5 * (b - a) + return [(c + h * 0.5 * (x[0] + x[1]), h * 0.5 * (w[0] + w[1])) for x, w in zip(xs, ws)] + s_nodes = [q for a, b in ((0.3, 1.0), (1.0, 1.4), (1.4, 1.8), (1.8, 2.4)) for q in nodes(a, b)] + v_nodes = [q for a, b in ((0.0, 0.5), (0.5, 0.8), (0.8, 1.0)) for q in nodes(a, b)] + t_nodes = [q for a, b in ((0.0, 0.7), (0.7, 1.2), (1.2, 2.5), (2.5, 6.0), (6.0, 15.0)) for q in nodes(a, b)] + for s, w1 in s_nodes: + for v, w2 in v_nodes: + inner = sum(w3 * _f3_safe_float(s, v, t) for t, w3 in t_nodes) + wgt = s ** 20 * v ** 4 * (1 - v ** 4) * math.exp(-s ** 8 * (2 * v ** 8 - v ** 4 + 2) / 10) + inner += wgt * 15.0 ** -7 / 7 + tot += w1 * w2 * inner + return 64.0 * tot + + +# ---------------------------------------------------------------- d = 2 +def _g2_safe_bound(box): + return CX.cabs_up(CU.g2(CU.CplxOps, box)) + + +def panel2(a, b, tol, nmin=4, m=24, nmax=60): + """Certified int_a^b g2 with the smallest admissible Gauss-Legendre rule (n >= nmin); returns (enclosure, n, rho, + M, err). For each rho the bound M is the smaller of the two representations' maxima over the slabs; a + representation counts only if it succeeds on every slab (NOTE.md, Lemma Q).""" + h = 0.5 * (b - a) + best = None + for rho in RHOS2: + slabs = G.ellipse_slabs(a, b, rho, m) + Ms = [] + for bf in (_g2_safe_bound, CU.g2_bound_direct): + try: + Ms.append(max(bf(bx) for bx in slabs)) + except d3.FAIL: + pass + if not Ms: + continue + M = min(Ms) + for n in range(nmin, nmax + 1): + e = G.error_bound(n, rho, M, h) + if e <= tol: + if best is None or (n, e) < (best[0], best[1]): + best = (n, e, rho, M) + break + if nmin > 4 and best is not None: + break + if best is None: + raise RuntimeError('d = 2 panel [%g, %g]: no admissible rule' % (a, b)) + n, e, rho, M = best + if nmin > 4: + n = nmin + e = G.error_bound(n, rho, M, h) + xs, ws = G.panel(n, a, b) + Q = (0.0, 0.0) + for x, w in zip(xs, ws): + Q = add(Q, mul(w, CU.g2(CU.RealOps, x))) + return (dn(Q[0] - e), up(Q[1] + e)), n, rho, M, e + + +def run_d2(log=print): + tot = (0.0, 0.0) + rules = [] + for a, b in zip(U_EDGES, U_EDGES[1:]): + iv, n, rho, M, e = panel2(a, b, TOL2) + tot = add(tot, iv) + rules.append([a, b, n, rho]) + # tail u > U_MAX: int g2 = U^-7/7 - int u^-8 Re H, |Re H(u^4)| <= (36/19)^(5/8) u^-7 < 1.5 u^-7 + U = U_MAX + main = div((1.0, 1.0), (7.0 * U ** 7, 7.0 * U ** 7)) + err = up(1.5 * U ** -14 / 14.0 * (1 + 1e-12)) + tail = (dn(main[0] - err), up(main[1] + err)) + tot = add(tot, tail) + log('d = 2: int g2 in [%.17g, %.17g]' % tot) + assembly_check('d = 2: int g2', tot, float_int_g2(), 1e-9) + return {'int_g2': tot, 'tail': tail, 'rules': rules} + + +# ---------------------------------------------------------------- d = 3 +def _box_job(b): + S, T = b + iv, info = d3.box(S, T, TOL3) + return iv, info + + +def run_d3(procs, log=print): + bx = d3.boxes() + t0 = time.time() + if procs > 1: + with Pool(procs) as pool: + res = pool.map(_box_job, bx, chunksize=1) + else: + res = [_box_job(b) for b in bx] + tot = (0.0, 0.0) + crude = gl = nodes = 0 + for iv, info in res: + tot = add(tot, iv) + if info.get('crude'): + crude += 1 + else: + gl += 1 + nodes += info['ns'] * info['nt'] * sum(info['nv']) + tt = d3.tail_t() + ts = d3.tail_s_bound() + T3 = scal(64.0, add(add(tot, tt), (0.0, ts))) + log('d = 3: T3 in [%.17g, %.17g] (%d boxes: %d crude, %d Gauss-Legendre, %d nodes; %.0fs)' % ( + T3[0], T3[1], len(bx), crude, gl, nodes, time.time() - t0)) + assembly_check('d = 3: T3', T3, float_T3(), 1e-6) + return {'T3': T3, 'boxes': scal(64.0, tot), 'tail_t': scal(64.0, tt), 'tail_s': (0.0, 64.0 * ts), + 'n_boxes': len(bx), 'n_crude': crude, + 'n_gl': gl, 'n_nodes': nodes} + + +# ---------------------------------------------------------------- assembly and publication +def certify(procs, log=print): + out = {} + out['d2'] = run_d2(log) + out['d3'] = run_d3(procs, log) + out['P2'] = K.prefactor_d2() + out['P3'] = K.prefactor_d3() + out['c1'] = {1: K.c1_d1(), 2: mul(out['P2'], out['d2']['int_g2']), 3: mul(out['P3'], out['d3']['T3'])} + out['cref'] = {2: K.c_ref(2), 3: K.c_ref(3)} + out['side24'] = {d: S24.transfer(d) for d in (2, 3)} + out['K_p'] = K.K_p() + out['Gamma'] = {'5/8': K.gamma(Fr(5, 8)), '11/4': K.gamma(Fr(11, 4)), '11/8': K.gamma(Fr(11, 8)), + '7/6': K.gamma(Fr(7, 6))} + return out + + +def _round_out(x, digits, upward): + if x == 0.0: + return '0' + from decimal import Decimal, ROUND_FLOOR, ROUND_CEILING + d = Decimal(x) + q = Decimal(1).scaleb(d.adjusted() - digits + 1) + return str(d.quantize(q, rounding=ROUND_CEILING if upward else ROUND_FLOOR)) + + +def pub(iv, digits=DIGITS): + lo = iv[0] - MARGIN * abs(iv[0]) + hi = iv[1] + MARGIN * abs(iv[1]) + return [_round_out(dn(lo), digits, False), _round_out(up(hi), digits, True)] + + +def derived(raw): + from decimal import Decimal + res = {} + for d in (1, 2, 3): + c1 = raw['c1'][d] + Ic = mul(K.I_CAND_FACTOR, c1) + res['d=%d' % d] = {'c1': c1, 'I_cand = (3^(1/4)/2) c1': Ic, 'I_cand - c1': sub(Ic, c1)} + if d in (2, 3): + res['d=%d' % d]['c_(d,ref)'] = raw['cref'][d] + res['d=%d' % d]['c1 / c_(d,ref)'] = div(c1, raw['cref'][d]) + sec = res['[P] field, every L >= 24 (Lemma S)'] = {} + for d in (2, 3): + eta = raw['side24'][d]['eta'] + fac = (from_frac(1 - eta)[0], from_frac(1 + eta)[1]) + c1L = mul(raw['c1'][d], fac) + IcL = mul(K.I_CAND_FACTOR, c1L) + lo, hi = S24.C_24[d] + c24 = (from_frac(Fr(lo))[0], from_frac(Fr(hi))[1]) + sec['d=%d: c1' % d] = c1L + sec['d=%d: I_cand' % d] = IcL + sec['d=%d: I_cand - c1' % d] = sub(IcL, c1L) + sec['d=%d: c_(d,24) (side24_v1, consumed)' % d] = c24 + sec['d=%d: c1 / c_(d,24), L = 24' % d] = div(c1L, c24) + res['parts'] = { + 'd=2: int_0^oo g2(u) du': raw['d2']['int_g2'], + 'd=2: tail u > %g' % U_MAX: raw['d2']['tail'], + 'd=2: P2': raw['P2'], + 'd=3: T3': raw['d3']['T3'], + 'd=3: boxes': raw['d3']['boxes'], + 'd=3: tail t > %g' % d3.T_MAX: raw['d3']['tail_t'], + 'd=3: tail s > %g' % d3.S_MAX: raw['d3']['tail_s'], + 'd=3: P3': raw['P3'], + 'K = pi/(2 Gamma(11/4) sin(7 pi/8))': raw['K_p'], + } + for k, v in raw['Gamma'].items(): + res['parts']['Gamma(%s)' % k] = v + return res + + +def document(raw): + der = derived(raw) + vals = {} + for sec, items in der.items(): + vals[sec] = {k: pub(v) for k, v in items.items()} + return { + 'object': PACKET, + 'scientific_effect': 'NONE', + 'statement': 'interval enclosures (outward decimals) of the cusp coefficient c1 of Math-#207 (CU.2) for the ' + 'Gaussian kernel, d = 1, 2, 3; see NOTE.md', + 'values': vals, + 'rules': { + 'd=2 panels [a, b, n, rho]': raw['d2']['rules'], + 'd=3': {'boxes': raw['d3']['n_boxes'], 'crude': raw['d3']['n_crude'], 'gauss_legendre': raw['d3']['n_gl'], + 'nodes': raw['d3']['n_nodes'], 'T_MAX': d3.T_MAX, 'S_MAX': d3.S_MAX}, + 'Lemma S (exact rationals)': { + 'd=%d' % d: {k: str(v) for k, v in raw['side24'][d].items()} for d in (2, 3)}, + }, + 'publication': {'margin_relative': MARGIN, 'digits': DIGITS}, + } + + +def compare(doc, raw): + bad = [] + new = document(raw) + if new['rules'] != doc['rules']: + bad.append('rule data differ: %r vs %r' % (new['rules']['d=3'], doc['rules']['d=3'])) + der = derived(raw) + from decimal import Decimal + for sec, items in der.items(): + for k, v in items.items(): + p = doc['values'].get(sec, {}).get(k) + if p is None: + bad.append('missing %s / %s' % (sec, k)) + continue + lo, hi = Decimal(p[0]), Decimal(p[1]) + if not (lo <= Decimal(v[0]) and Decimal(v[1]) <= hi): + bad.append('%s / %s: raw [%r, %r] not inside published %r' % (sec, k, v[0], v[1], p)) + slack = Decimal(3 * MARGIN) + if (Decimal(v[0]) - lo) > slack * abs(Decimal(v[0])) + Decimal('1e-300') or \ + (hi - Decimal(v[1])) > slack * abs(Decimal(v[1])) + Decimal('1e-300'): + bad.append('%s / %s: published %r looser than the replay [%r, %r]' % (sec, k, p, v[0], v[1])) + return bad + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument('--write', action='store_true') + ap.add_argument('--check', 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() + if a.mutant is not None and a.mutant not in MUTANTS: + print('unknown mutant', a.mutant) + return 2 + install_mutant(a.mutant) + fails = selftest() + print('self-test: %d failures' % len(fails)) + if fails: + for f in fails[:10]: + print('SELF-TEST FAILURE', f) + return 1 + if not (a.write or a.check): + return 0 + raw = certify(a.procs) + path = os.path.join(HERE, 'RESULTS.json') + if a.write: + with open(path, 'w') as fh: + json.dump(document(raw), fh, indent=1, sort_keys=False) + fh.write('\n') + print('wrote', path) + return 0 + 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 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/frontiers/cusp_coefficient_certified_20261001/consts.py b/frontiers/cusp_coefficient_certified_20261001/consts.py new file mode 100644 index 00000000..5cac30c3 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/consts.py @@ -0,0 +1,115 @@ +"""Rigorous constants for CL-CU-CUSP-COEFFICIENT-20261001: Gamma at rationals, the prefactors of c1 in d = 1, 2, 3 +(NOTE.md sections 2-4) and the reference leading coefficients c_(d,ref) of coefficients/side24_v1 (1). + +log Gamma(y) at y = q + N >= 10 (the shift used by lgamma_iv; the bound below holds for every real y > 0) uses Stirling's series + (y - 1/2) log y - y + log(2 pi)/2 + sum_(k=1..K) B_(2k)/(2k (2k-1) y^(2k-1)) + R_K, + |R_K| <= |B_(2K+2)|/((2K+2)(2K+1) y^(2K+1)) (real y > 0: the remainder is bounded by the first omitted term), +and the recurrence log Gamma(q) = log Gamma(q + N) - log prod_(j 0 (Stirling at y = q + N >= 10; remainder below 1e-20).""" + q = Fr(q) + y = q + N + ly = _log_rat(y) + Y = from_frac(y) + s = sub(mul(sub(Y, (0.5, 0.5)), ly), Y) + s = add(s, scal(0.5, LOG_2PI)) + corr = Fr(0) + for k in range(1, K + 1): + corr += BERN[k - 1] / (2 * k * (2 * k - 1) * y ** (2 * k - 1)) + rem = abs(BERN[K]) / ((2 * K + 2) * (2 * K + 1) * y ** (2 * K + 1)) + s = add(s, (from_frac(corr - rem)[0], from_frac(corr + rem)[1])) + prod = Fr(1) + for j in range(N): + prod *= q + j + return sub(s, _log_rat(prod)) + + +def gamma(q): + return exp_iv(lngamma(q)) + + +def rpow(base, expo): + """base^expo for a positive interval base and a rational exponent, via exp(expo log base).""" + return exp_iv(mul(from_frac(Fr(expo)), log_iv(base))) + + +def ipow(base, n): + out = (1.0, 1.0) + for _ in range(abs(n)): + out = mul(out, base) + return out if n >= 0 else div((1.0, 1.0), out) + + +def c(q): + return from_frac(Fr(q)) + + +SIN_PI_8 = scal(0.5, isqrt(sub((2.0, 2.0), isqrt((2.0, 2.0))))) # sin(7 pi/8) = sin(pi/8) = sqrt(2 - sqrt2)/2 + + +def K_p(): + """K = int_0^oo (1 - cos u) u^(-11/4) du = pi/(2 Gamma(11/4) sin(7 pi/8)).""" + return div(PI, mul((2.0, 2.0), mul(gamma(Fr(11, 4)), SIN_PI_8))) + + +def prefactor_d2(): + """c1^(2) = P2 * int_0^oo g2(u) du with P2 = -(96/7) 2^(1/4) (2 pi)^-2 12^(-1/2) Gamma(5/8) Gamma(11/4) sin(7 pi/8)/pi + * 4 C0, C0 = (3/16)^(-5/8).""" + P = mul(c(Fr(-96, 7)), rpow((2.0, 2.0), Fr(1, 4))) + P = mul(P, div((1.0, 1.0), sqr(TWO_PI))) + P = mul(P, rpow((12.0, 12.0), Fr(-1, 2))) + P = mul(P, mul(gamma(Fr(5, 8)), gamma(Fr(11, 4)))) + P = mul(P, div(SIN_PI_8, PI)) + C0 = rpow(c(Fr(3, 16)), Fr(-5, 8)) + return mul(P, scal(4.0, C0)) + + +def prefactor_d3(): + """c1^(3) = P3 * T3 with P3 = -(1536/7) 2^(1/4) pi C_A Gamma(11/4) sin(7 pi/8), + C_A = (2 pi)^(-11/2) 12^(-1/2) sqrt(pi/5).""" + CA = mul(rpow(TWO_PI, Fr(-11, 2)), rpow((12.0, 12.0), Fr(-1, 2))) + CA = mul(CA, isqrt(div(PI, (5.0, 5.0)))) + P = mul(c(Fr(-1536, 7)), rpow((2.0, 2.0), Fr(1, 4))) + P = mul(P, mul(PI, CA)) + return mul(P, mul(gamma(Fr(11, 4)), SIN_PI_8)) + + +def c1_d1(): + """c1^(1) = -(192/7) 2^(1/4) * 2 * 12^(-7/4) (2 pi)^-2 12^(-1/2) sqrt(4 pi/3) 30^(7/8) 2^(7/8) Gamma(11/8)/sqrt(pi).""" + v = mul(c(Fr(-384, 7)), rpow((2.0, 2.0), Fr(1, 4))) + v = mul(v, rpow((12.0, 12.0), Fr(-9, 4))) + v = mul(v, div((1.0, 1.0), sqr(TWO_PI))) + v = mul(v, isqrt(div(scal(4.0, PI), (3.0, 3.0)))) + v = mul(v, rpow((30.0, 30.0), Fr(7, 8))) + v = mul(v, rpow((2.0, 2.0), Fr(7, 8))) + v = mul(v, div(gamma(Fr(11, 8)), isqrt(PI))) + return v + + +def c_ref(d): + """c_(d,ref) = Gamma(7/6) (3/2)^(1/3) D_(d-1)/(2 sqrt3 pi^(d-1) sqrt(pi)), D_1 = 4/3, D_2 = 29/6 - sqrt6 + (coefficients/side24_v1/PROOF.md (1)).""" + D = c(Fr(4, 3)) if d == 2 else sub(c(Fr(29, 6)), isqrt((6.0, 6.0))) + v = mul(gamma(Fr(7, 6)), rpow(c(Fr(3, 2)), Fr(1, 3))) + v = mul(v, D) + den = mul(scal(2.0, isqrt((3.0, 3.0))), mul(ipow(PI, d - 1), isqrt(PI))) + return div(v, den) + + +I_CAND_FACTOR = scal(0.5, rpow((3.0, 3.0), Fr(1, 4))) # I^cand = (3^(1/4)/2) c1 diff --git a/frontiers/cusp_coefficient_certified_20261001/controls.py b/frontiers/cusp_coefficient_certified_20261001/controls.py new file mode 100644 index 00000000..5e41d95c --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/controls.py @@ -0,0 +1,211 @@ +"""Floating controls for CL-CU-CUSP-COEFFICIENT-20261001. These are not part of the certificate. + +1. d = 2 by a route that does not use the Mellin identity. Given (a = -A, gamma, b), Y is Gaussian with + sigma = a/sqrt6 and mean (a b - gamma^2)/4, so E|Y|^(7/4) = sigma^(7/4) M(mu/sigma), M(t) = E|Z + t|^(7/4) (Kummer + series). b ~ N(a/4, 1/2) under the weight, so the b-average is again of M type with variance factor 19/16. This + leaves a 2D integral over (a, gamma). +2. d = 3 by the same route: eigenvalues a1 > a2 of -A, the Vandermonde, gamma ~ N(0, 2 I), and the b-average with + factor 23/20. This leaves a 4D integral, done here by a floating product Gauss rule. +3. Math-#207's floating values and Monte Carlo runs (head 826b8be, PROOF.md section 8, remark 3) and the d = 1 value of + Math-#214/#216. +4. Lemma S: the actual deviations D^g K_24(0) - D^g phi(0) of the jet covariance entries (standard frame, d = 3), by a + direct floating lattice sum over |n_i| <= 2, against the certified entry bound E. + +Run: python3 -B controls.py (pure Python, a few minutes).""" +import json +import math +import os +import sys + +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +sys.dont_write_bytecode = True +import gauss as G + +P = 1.75 +CP = 2 ** (P / 2) * math.gamma((P + 1) / 2) / math.sqrt(math.pi) + + +def M(t): + """E|Z + t|^(7/4): Kummer's 1F1((1+p)/2; 1/2; t^2/2) e^(-t^2/2) series for |t| <= 26, else the asymptotic series + t^p sum_j binom(p, 2j) (2j - 1)!! t^(-2j) through j = 4 (relative error below 1e-14 there).""" + t = abs(t) + if t > 26.0: + c2 = P * (P - 1) / 2 + c4 = P * (P - 1) * (P - 2) * (P - 3) / 8 + c6 = P * (P - 1) * (P - 2) * (P - 3) * (P - 4) * (P - 5) / 48 + c8 = P * (P - 1) * (P - 2) * (P - 3) * (P - 4) * (P - 5) * (P - 6) * (P - 7) / 384 + return t ** P * (1 + c2 / t ** 2 + c4 / t ** 4 + c6 / t ** 6 + c8 / t ** 8) + z = 0.5 * t * t + term = 1.0 + s = 1.0 + a, b = (1 + P) / 2, 0.5 + n = 0 + while True: + term *= (a + n) / (b + n) * z / (n + 1) + s += term + n += 1 + if term < 1e-17 * s and n > z: + break + return CP * math.exp(-z) * s + + +def gl(n, a, b): + xs, ws = G.rule(n) + c, h = 0.5 * (a + b), 0.5 * (b - a) + return [(c + h * 0.5 * (x[0] + x[1]), h * 0.5 * (w[0] + w[1])) for x, w in zip(xs, ws)] + + +def gamma_nodes(scale, gmax=12.0, n=16): + """Composite nodes for gamma in [0, gmax], refined near 0 at the scale sqrt(a) of the integrand.""" + edges = [0.0] + e = 0.25 * scale + while e < gmax: + edges.append(e) + e *= 2.0 + edges.append(gmax) + out = [] + for a, b in zip(edges, edges[1:]): + out.extend(gl(n, a, b)) + return out + + +def phi2(g): + """density of N(0, 2)""" + return math.exp(-g * g / 4) / math.sqrt(4 * math.pi) + + +def control_d2(): + c = math.sqrt(6 / 19) + tot = 0.0 + for s, ws in gl(80, 0.0, 2.7): + a = s ** 4 + jac = 4 * s ** 3 + Eg = 0.0 + for g, wg in gamma_nodes(math.sqrt(a)): + Eg += 2 * wg * phi2(g) * M(c * (a / 4 - g * g / a)) + tot += ws * jac * a * a * math.exp(-3 * a * a / 16) * Eg + pref = -(192 / 7) * 2 ** 0.25 * 2 * math.pi * (2 * math.pi) ** -3 * 12 ** -0.5 * 0.5 * (19 / 16) ** (7 / 8) * 6 ** (-7 / 8) + return pref * tot + + +def control_d3(): + c = math.sqrt(15 / 46) + tot = 0.0 + s_nodes = gl(24, 0.0, 1.1) + gl(24, 1.1, 1.6) + gl(24, 1.6, 2.6) + for s, ws in s_nodes: + a1 = s ** 4 + for w, ww in gl(28, 0.0, 1.0): + a2 = a1 * w ** 4 + if a2 <= 0.0: + continue + jac = 4 * s ** 3 * a1 * 4 * w ** 3 + g1n = gamma_nodes(math.sqrt(a1), n=10) + g2n = gamma_nodes(math.sqrt(a2), n=10) + Eg = 0.0 + for g1, w1 in g1n: + p1 = w1 * 2 * phi2(g1) + x1 = g1 * g1 / a1 + for g2, w2 in g2n: + Eg += p1 * w2 * 2 * phi2(g2) * M(c * (x1 + g2 * g2 / a2 - (a1 + a2) / 5)) + tot += 2 * ws * ww * jac * (a1 - a2) * (a1 * a2) ** 2 * math.exp((a1 + a2) ** 2 / 20 - (a1 * a1 + a2 * a2) / 4) * Eg + pref = (-(192 / 7) * 2 ** 0.25 * 4 * math.pi * (2 * math.pi) ** -4 * 12 ** -0.5 * 6 ** (-7 / 8) * (2 * math.pi) ** -1.5 + * 0.5 * math.sqrt(4 * math.pi / 5) * (23 / 20) ** (7 / 8) * math.pi * 0.5) + return pref * tot + + +def _he(n, x): + """Probabilists' Hermite polynomial He_n(x), so that d^n/dx^n e^(-x^2/2) = (-1)^n He_n(x) e^(-x^2/2).""" + a, b = 1.0, x + if n == 0: + return a + for k in range(1, n): + a, b = b, x * b - k * a + return b + + +def control_side24(L=24.0, d=3): + """max over the covariance entries of |D^g K_L(0) - D^g phi(0)|, K_L normalized to variance one.""" + import side24 as S24 + J = S24.jets(d) + gs = sorted({tuple(x + y for x, y in zip(a, b)) for _, a, _ in J for _, b, _ in J}) + rng = range(-2, 3) + pts = [(i, j, k) for i in rng for j in rng for k in rng] if d == 3 else [(i, j) for i in rng for j in rng] + pts = [n for n in pts if any(n)] + S = 1.0 + sum(math.exp(-0.5 * L * L * sum(c * c for c in n)) for n in pts) + worst = (0.0, None) + for g in gs: + far = 0.0 + for n in pts: + v = 1.0 + for gi, ni in zip(g, n): + x = L * ni + v *= (-1) ** gi * _he(gi, x) * math.exp(-0.5 * x * x) + far += v + dev = abs(far / S + S24.dK0(g) * (1.0 / S - 1.0)) + if dev > worst[0]: + worst = (dev, g) + return worst + + +REFERENCE = { + 'Math-#207 (826b8be) PROOF.md section 8: c1, d = 2': -0.26939883, + 'Math-#207 (826b8be) PROOF.md section 8: c1, d = 3': -0.21184835, + 'Math-#207: c1/c, d = 2': -3.669938, + 'Math-#207: c1/c, d = 3': -5.071062, + 'Math-#207: I_cand, d = 2': -0.1772744, + 'Math-#207: I_cand - c1, d = 2': 0.0921244, + 'Math-#207: I_cand, d = 3': -0.1394041, + 'Math-#207: I_cand - c1, d = 3': 0.0724443, + 'Math-#207: c, d = 2': 0.0734069193, + 'Math-#207: c, d = 3': 0.0417759318, + 'Math-#207: shipped fixed-seed Monte Carlo c1, d = 3': (-0.21163, 0.00014), + 'Math-#207: development Monte Carlo c1, d = 3 (2e8 samples)': (-0.21188, 0.00010), + 'Math-#207: referee Monte Carlo c1, d = 3 (5e8 samples)': (-0.211848, 0.000024), + 'Math-#216 table: c1, d = 1 (= Math-#214 C_1)': -0.227606, +} + + +def main(): + with open(os.path.join(HERE, 'RESULTS.json')) as fh: + res = json.load(fh)['values'] + def mid(p): + return 0.5 * (float(p[0]) + float(p[1])) + print('1. d = 2, independent 2D route (no Mellin identity):') + c2 = control_d2() + cert = res['d=2']['c1'] + print(' control %.15f certified [%s, %s] difference %.1e' % (c2, cert[0], cert[1], c2 - mid(cert))) + print('2. d = 3, independent 4D route (no Mellin identity):') + c3 = control_d3() + cert = res['d=3']['c1'] + print(' control %.15f certified [%s, %s] difference %.1e' % (c3, cert[0], cert[1], c3 - mid(cert))) + print('3. Floating values of the sources against the certified enclosures:') + rows = [ + ('Math-#207 (826b8be) PROOF.md section 8: c1, d = 2', res['d=2']['c1']), + ('Math-#207 (826b8be) PROOF.md section 8: c1, d = 3', res['d=3']['c1']), + ('Math-#207: c, d = 2', res['d=2']['c_(d,ref)']), + ('Math-#207: c, d = 3', res['d=3']['c_(d,ref)']), + ('Math-#207: c1/c, d = 2', res['d=2']['c1 / c_(d,ref)']), + ('Math-#207: c1/c, d = 3', res['d=3']['c1 / c_(d,ref)']), + ('Math-#207: I_cand, d = 2', res['d=2']['I_cand = (3^(1/4)/2) c1']), + ('Math-#207: I_cand - c1, d = 2', res['d=2']['I_cand - c1']), + ('Math-#207: I_cand, d = 3', res['d=3']['I_cand = (3^(1/4)/2) c1']), + ('Math-#207: I_cand - c1, d = 3', res['d=3']['I_cand - c1']), + ('Math-#216 table: c1, d = 1 (= Math-#214 C_1)', res['d=1']['c1']), + ] + for name, cert in rows: + v = REFERENCE[name] + print(' %-52s %12s certified mid %+.10f rounded value differs by %.1e' % (name, v, mid(cert), v - mid(cert))) + for name in ('Math-#207: shipped fixed-seed Monte Carlo c1, d = 3', + 'Math-#207: development Monte Carlo c1, d = 3 (2e8 samples)', + 'Math-#207: referee Monte Carlo c1, d = 3 (5e8 samples)'): + mc, se = REFERENCE[name] + print(' %-58s %.6f +- %.6f -> (MC - certified)/s.e. = %+.2f' % (name, mc, se, (mc - mid(res['d=3']['c1'])) / se)) + import side24 as S24 + dev, g = control_side24() + print('4. Lemma S, d = 3, L = 24: largest |D^g K_24(0) - D^g phi(0)| = %.3e (at g = %s); certified entry bound E = %.3e' + % (dev, g, float(S24.image_bound()))) + + +if __name__ == '__main__': + main() diff --git a/frontiers/cusp_coefficient_certified_20261001/cusp.py b/frontiers/cusp_coefficient_certified_20261001/cusp.py new file mode 100644 index 00000000..fba24e5e --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/cusp.py @@ -0,0 +1,189 @@ +"""The integrands of CL-CU-CUSP-COEFFICIENT-20261001, written once against an operations object: +RealOps (sharp real intervals, for the quadrature sums) and CplxOps (complex boxes, for the Bernstein-ellipse bounds). + +d = 2 (NOTE.md section 3). With omega = u^4 and + H(omega) = (1 + i omega)^(-1/2) (1 - i omega/3 + 19 omega^2/36)^(-5/8), +let R = Re log H and I = Im log H. Then + R/omega^2 = -(1/4) L(omega^2) - (5/16) q L(omega^2 q), q = 7/6 + 361 omega^2/1296, + I/omega = -(1/2) A(omega) + (5/8) A(omega/v)/v, v = 3 + 19 omega^2/12, +and + g2(u) = (1 - Re H)/omega^2 = -E1(R) (R/omega^2) + e^R Sinc(I/2)^2 (I/omega)^2/2. +Moreover J = 4 C0 int_0^oo g2(u) du with C0 = (3/16)^(-5/8). + +d = 3 (NOTE.md section 4). With a1 = s^4 >= a2 = a1 v^4 the transverse eigenvalue moduli and omega = t^4/a1: + Rt = -(1/4)(L(t^8) + v^8 L(t^8 v^8)) - 23 a1^2 v^8/240, R = t^8 Rt, + It = (1/2)(A(t^4) + v^4 A(t^4 v^4)) - a1^2 v^4 (1 + v^4)/20, I = t^4 It, + Gt = -E1(R) Rt + e^R Sinc(I/2)^2 It^2/2, + f3(s, v, t) = s^20 v^4 (1 - v^4) exp(-s^8 (2 v^8 - v^4 + 2)/10) Gt. +Moreover T3 = 64 int_0^oo ds int_0^1 dv int_0^oo dt f3.""" +from fractions import Fraction as Fr +from ia import dn, up, add, sub, neg, mul, div, sqr, pw, from_frac +import elem as EL +import cx as CX + +_CONST = {} +MUTANT = None # 'phase-sign' flips the sign of the a1^2 phase term of f3 (certificate.py mutants) + + +def _c(q): + q = Fr(q) + r = _CONST.get(q) + if r is None: + r = from_frac(q) + _CONST[q] = r + return r + + +class RealOps: + @staticmethod + def c(q): + return _c(q) + add = staticmethod(add) + sub = staticmethod(sub) + mul = staticmethod(mul) + div = staticmethod(div) + neg = staticmethod(neg) + sq = staticmethod(sqr) + pw = staticmethod(pw) + exp = staticmethod(EL.exp_iv) + E1 = staticmethod(EL.E1_iv) + L = staticmethod(EL.L_iv) + A = staticmethod(EL.A_iv) + Sinc = staticmethod(EL.Sinc_iv) + + +class CplxOps: + @staticmethod + def c(q): + return (_c(q), (0.0, 0.0)) + add = staticmethod(CX.cadd) + sub = staticmethod(CX.csub) + mul = staticmethod(CX.cmul) + div = staticmethod(CX.cdiv) + + @staticmethod + def neg(z): + return (neg(z[0]), neg(z[1])) + + @staticmethod + def sq(z): + return CX.cmul(z, z) + pw = staticmethod(CX.cpow) + exp = staticmethod(CX.cexp) + E1 = staticmethod(CX.E1c) + L = staticmethod(CX.Lc) + A = staticmethod(CX.Ac) + Sinc = staticmethod(CX.Sincc) + + +def _G(o, R, I, Rt, It): + t1 = o.neg(o.mul(o.E1(R), Rt)) + t2 = o.mul(o.mul(o.exp(R), o.sq(o.Sinc(o.mul(o.c(Fr(1, 2)), I)))), o.mul(o.c(Fr(1, 2)), o.sq(It))) + return o.add(t1, t2) + + +def g2(o, u): + w = o.pw(u, 4) + w2 = o.pw(u, 8) + q = o.add(o.c(Fr(7, 6)), o.mul(o.c(Fr(361, 1296)), w2)) + Rt = o.neg(o.add(o.mul(o.c(Fr(1, 4)), o.L(w2)), o.mul(o.mul(o.c(Fr(5, 16)), q), o.L(o.mul(w2, q))))) + v = o.add(o.c(3), o.mul(o.c(Fr(19, 12)), w2)) + It = o.sub(o.mul(o.c(Fr(5, 8)), o.div(o.A(o.div(w, v)), v)), o.mul(o.c(Fr(1, 2)), o.A(w))) + return _G(o, o.mul(w2, Rt), o.mul(w, It), Rt, It) + + +def f3(o, s, v, t): + s4 = o.pw(s, 4) + s8 = o.pw(s, 8) + v4 = o.pw(v, 4) + v8 = o.pw(v, 8) + t4 = o.pw(t, 4) + t8 = o.pw(t, 8) + Rt = o.neg(o.add(o.mul(o.c(Fr(1, 4)), o.add(o.L(t8), o.mul(v8, o.L(o.mul(t8, v8))))), + o.mul(o.c(Fr(23, 240)), o.mul(s8, v8)))) + ph = o.mul(o.c(Fr(1, 20)), o.mul(o.mul(s8, v4), o.add(o.c(1), v4))) + half = o.mul(o.c(Fr(1, 2)), o.add(o.A(t4), o.mul(v4, o.A(o.mul(t4, v4))))) + It = o.add(half, ph) if MUTANT == 'phase-sign' else o.sub(half, ph) + G = _G(o, o.mul(t8, Rt), o.mul(t4, It), Rt, It) + E = o.mul(o.c(Fr(1, 10)), o.add(o.sub(o.mul(o.c(2), v8), v4), o.c(2))) + wgt = o.mul(o.mul(o.pw(s, 20), o.mul(v4, o.sub(o.c(1), v4))), o.exp(o.neg(o.mul(s8, E)))) + return o.mul(wgt, G) + + +# ---------------------------------------------------------------- direct bounds for the Bernstein ellipses +# On the real axis Re H = (H + H*)/2 with H*(w) = conj(H(conj w)), so the analytic continuation of g2 is +# D(u) = (1 - (H(w) + H*(w))/2)/w^2, w = u^4, with principal powers. It is analytic wherever no base of a principal +# power meets the closed half-line (-oo, 0] and w != 0. For a real exponent, |z^a| = |z|^a, hence +# |D(u)| <= (1 + (|H(w)| + |H*(w)|)/2)/|w|^2. +# A box passes _off_cut when it cannot meet (-oo, 0]. When every box covering an ellipse passes, D is analytic on +# the ellipse and equals the integrand there (both are analytic and agree on the real panel). The same holds for the +# d = 3 direct form. + +def _off_cut(z): + return z[0][0] > 0.0 or z[1][0] > 0.0 or z[1][1] < 0.0 + + +def _pow_abs_up(z, a): + """Upper bound of |z|^a over a box, a real; requires |z| > 0 when a < 0.""" + lo, hi = CX.cabs_dn(z), CX.cabs_up(z) + if a < 0: + if not lo > 0.0: + raise ZeroDivisionError('modulus may vanish') + return EL.exp_iv(mul(_c(Fr(a)), EL.log_iv((lo, lo))))[1] + return EL.exp_iv(mul(_c(Fr(a)), EL.log_iv((hi, hi))))[1] + + +def g2_bound_direct(u): + """Upper bound of |D(u)| over the complex box u (raises if the branch conditions cannot be certified).""" + w = CX.cpow(u, 4) + iw = CX.cmul(((0.0, 0.0), (1.0, 1.0)), w) + one = ((1.0, 1.0), (0.0, 0.0)) + w2 = CX.cpow(u, 8) + quad = CX.cscal(_c(Fr(19, 36)), w2) + third = CX.cscal(_c(Fr(1, 3)), iw) + bases_h = (CX.cadd(one, iw), CX.cadd(CX.csub(one, third), quad)) + bases_hs = (CX.csub(one, iw), CX.cadd(CX.cadd(one, third), quad)) + for b in bases_h + bases_hs: + if not _off_cut(b): + raise ValueError('branch condition not certified') + habs = up(_pow_abs_up(bases_h[0], -0.5) * _pow_abs_up(bases_h[1], -0.625)) + hsabs = up(_pow_abs_up(bases_hs[0], -0.5) * _pow_abs_up(bases_hs[1], -0.625)) + num = up(1.0 + up(0.5 * up(habs + hsabs))) + wmin = CX.cabs_dn(w) + if not wmin > 0.0: + raise ZeroDivisionError('w may vanish') + return up(num / dn(wmin * wmin)) + + +def f3_bound_direct(s, v, t): + """Upper bound of the modulus of the d = 3 integrand's direct continuation over complex boxes (s, v, t): + f3 = wgt (1 - (Phi + Phi*)/2)/t^8, + Phi = (1 - i t^4)^(-1/2) (1 - i t^4 v^4)^(-1/2) exp(-i t^4 c1 - t^8 c2), + c1 = a1^2 v^4 (1 + v^4)/20, c2 = 23 a1^2 v^8/240, a1 = s^4, and Phi* takes i -> -i.""" + one = ((1.0, 1.0), (0.0, 0.0)) + ii = ((0.0, 0.0), (1.0, 1.0)) + s4, s8 = CX.cpow(s, 4), CX.cpow(s, 8) + v4, v8 = CX.cpow(v, 4), CX.cpow(v, 8) + t4, t8 = CX.cpow(t, 4), CX.cpow(t, 8) + c1 = CX.cscal(_c(Fr(1, 20)), CX.cmul(CX.cmul(s8, v4), CX.cadd(one, v4))) + c2 = CX.cscal(_c(Fr(23, 240)), CX.cmul(s8, v8)) + it4 = CX.cmul(ii, t4) + it4v4 = CX.cmul(it4, v4) + total = 0.0 + for sign in (1, -1): + b1 = CX.csub(one, it4) if sign == 1 else CX.cadd(one, it4) + b2 = CX.csub(one, it4v4) if sign == 1 else CX.cadd(one, it4v4) + for b in (b1, b2): + if not _off_cut(b): + raise ValueError('branch condition not certified') + ex = CX.csub(CX.cscal((-float(sign), -float(sign)), CX.cmul(it4, c1)), CX.cmul(t8, c2)) + mod = up(up(_pow_abs_up(b1, -0.5) * _pow_abs_up(b2, -0.5)) * EL.exp_iv(ex[0])[1]) + total = up(total + mod) + num = up(1.0 + up(0.5 * total)) + tmin = CX.cabs_dn(t4) + if not tmin > 0.0: + raise ZeroDivisionError('t may vanish') + G = up(num / dn(tmin * tmin)) + E = CX.cscal(_c(Fr(1, 10)), CX.cadd(CX.csub(CX.cscal((2.0, 2.0), v8), v4), CX.cscal((2.0, 2.0), one))) + wgt = CX.cmul(CX.cmul(CX.cpow(s, 20), CX.cmul(v4, CX.csub(one, v4))), CX.cexp(CX.cscal((-1.0, -1.0), CX.cmul(s8, E)))) + return up(CX.cabs_up(wgt) * G) diff --git a/frontiers/cusp_coefficient_certified_20261001/cx.py b/frontiers/cusp_coefficient_certified_20261001/cx.py new file mode 100644 index 00000000..29e3a4a8 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/cx.py @@ -0,0 +1,159 @@ +"""Complex boxes for the analytic-continuation bounds of CL-CU-CUSP-COEFFICIENT-20261001. + +A complex box is a pair (re, im) of real intervals (ia.py conventions) and stands for {x + iy : x in re, y in im}. +Every operation returns a box containing all results. These boxes are used only to bound |f| on Bernstein ellipses; +the quadrature sums themselves are real-interval computations. + +The four kernels have integral representations over tau in [0, 1]: + + E1(z) = int e^(tau z), L(z) = int 1/(1 + tau z), A(z) = int 1/(1 + tau^2 z^2), Sinc(z) = int cos(tau z). + +So each value lies in the convex hull of the integrand's range. The box of the integrand with tau replaced by the +interval [0, 1] contains that range, and boxes are convex. L and A are analytic where 1 + tau z or 1 + tau^2 z^2 +stays away from 0 for all tau; cdiv raises if a denominator box meets 0. + +For larger |z| these hulls are needlessly wide, so the closed forms are used instead: +- L(z) = Log(1 + z)/z when Re(1 + z) > 0 (principal Log, Arg(1 + z) = atan(Im/Re)); +- A(z) = (pi/2 - A(1/z)/z)/z when Re z > 0, since atan z = pi/2 - atan(1/z) there; +- E1(z) = (e^z - 1)/z; +- Sinc(z) = sin z/z. +Each is analytic on the boxes where it is used. Both forms enclose the same analytic function, so either is a valid +bound; the closed form is chosen when min |z| >= 1.""" +from ia import dn, up, add, sub, neg, mul, div, scal, sqr, isqrt, PI +from elem import exp_iv, sin_iv, cos_iv, sinh_iv, cosh_iv, log_iv, atan_iv + +ZERO = (0.0, 0.0) +ONE = (1.0, 1.0) +TAU = (0.0, 1.0) + + +def creal(x): + return (x, ZERO) + + +def cadd(a, b): + return (add(a[0], b[0]), add(a[1], b[1])) + + +def csub(a, b): + return (sub(a[0], b[0]), sub(a[1], b[1])) + + +def cmul(a, b): + return (sub(mul(a[0], b[0]), mul(a[1], b[1])), add(mul(a[0], b[1]), mul(a[1], b[0]))) + + +def cscal(r, a): + """real interval r times complex box a.""" + return (mul(r, a[0]), mul(r, a[1])) + + +def cabs2(a): + return add(sqr(a[0]), sqr(a[1])) + + +def cdiv(a, b): + d = cabs2(b) + if d[0] <= 0.0: + raise ZeroDivisionError('complex box division by a box that may contain 0') + num = cmul(a, (b[0], neg(b[1]))) + return (div(num[0], d), div(num[1], d)) + + +def cexp(a): + m = exp_iv(a[0]) + return (mul(m, cos_iv(a[1])), mul(m, sin_iv(a[1]))) + + +def ccos(a): + """cos(x + iy) = cos x cosh y - i sin x sinh y.""" + return (mul(cos_iv(a[0]), cosh_iv(a[1])), neg(mul(sin_iv(a[0]), sinh_iv(a[1])))) + + +def _abs2_clamped(a): + d = cabs2(a) + return (max(0.0, d[0]), d[1]) + + +def cabs_up(a): + """An upper bound of |z| over the box.""" + return isqrt(_abs2_clamped(a))[1] + + +def cabs_dn(a): + """A lower bound of |z| over the box.""" + return isqrt(_abs2_clamped(a))[0] + + +def cpow(a, n): + """z^n for an integer n >= 1. On a box in the right half-plane Re > 0, the polar form r^n (cos n theta, + sin n theta) is used: r over [min |z|, max |z|], and theta = atan(y/x) over the box, which is the interval atan of + the interval quotient. This avoids the wrapping of repeated box products. Elsewhere it uses repeated products.""" + if a[0][0] > 0.0: + r = (cabs_dn(a), cabs_up(a)) + th = atan_iv(div(a[1], a[0])) + rn = _pow_iv(r, n) + nth = (dn(n * th[0]), up(n * th[1])) + return (mul(rn, cos_iv(nth)), mul(rn, sin_iv(nth))) + out = a + for _ in range(n - 1): + out = cmul(out, a) + return out + + +def _pow_iv(r, n): + """[r0^n, r1^n] for 0 <= r0 <= r1, outward.""" + lo = 1.0 + hi = 1.0 + for _ in range(n): + lo = dn(lo * r[0]) + hi = up(hi * r[1]) + return (max(0.0, lo), hi) + + +def csin(a): + """sin(x + iy) = sin x cosh y + i cos x sinh y.""" + return (mul(sin_iv(a[0]), cosh_iv(a[1])), mul(cos_iv(a[0]), sinh_iv(a[1]))) + + +def clog_right(a): + """Principal Log over a box in the right half-plane Re > 0.""" + if not a[0][0] > 0.0: + raise ValueError('clog_right needs Re > 0') + mod = scal(0.5, log_iv(cabs2(a))) + ratio = div(a[1], a[0]) + return (mod, atan_iv(ratio)) + + +def _big(z): + return cabs_dn(z) >= 1.0 + + +def E1c(z): + if _big(z): + return cdiv(csub(cexp(z), creal(ONE)), z) + return cexp(cscal(TAU, z)) + + +def Lc(z): + if _big(z) and z[0][0] > -0.5: + return cdiv(clog_right(cadd(creal(ONE), z)), z) + return cdiv(creal(ONE), cadd(creal(ONE), cscal(TAU, z))) + + +def _Ac_hull(z): + return cdiv(creal(ONE), cadd(creal(ONE), cscal(TAU, cmul(z, z)))) + + +def Ac(z): + if _big(z) and z[0][0] > 0.0: + w = cdiv(creal(ONE), z) + atan_w = cmul(w, _Ac_hull(w)) + return cdiv(csub(creal(scal(0.5, PI)), atan_w), z) + return _Ac_hull(z) + + +def Sincc(z): + if _big(z): + return cdiv(csin(z), z) + return ccos(cscal(TAU, z)) diff --git a/frontiers/cusp_coefficient_certified_20261001/d3.py b/frontiers/cusp_coefficient_certified_20261001/d3.py new file mode 100644 index 00000000..bae66453 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/d3.py @@ -0,0 +1,276 @@ +"""The d = 3 integral T3 of CL-CU-CUSP-COEFFICIENT-20261001 (NOTE.md section 4): + + T3 = 64 int_0^oo ds int_0^oo dt int_0^1 dv f3(s, v, t) (cusp.f3). + +Order of integration. v is innermost: for fixed real (s, t), f3 is analytic in v on the ellipses used. The outer +variables (s, t) run over boxes S x T with t <= T_MAX. The tail t > T_MAX is done separately (tail_integral). + +Box enclosure (Lemma B of NOTE.md). Let Q = Q_S Q_T Q_V be the tensor Gauss-Legendre rule (Q_V composite over the +v-panels). Then + |I - Q| <= |(I_S - Q_S) I_T I_V f| + |Q_S (I_T - Q_T) I_V f| + |Q_S Q_T (I_V - Q_V) f|. +Each term takes the one-dimensional bound of gauss.error_bound in its own direction: +- for s and t, the bound uses sup |I_V f| on the ellipse in that direction, with the other outer variable real; +- for v, it uses sup |f| on each v-ellipse, with s and t real. +The sums of positive weights are |S| and |T|. Each sup is taken over complex boxes: ellipse slabs times real +sub-intervals of the other variables. For each (sub-interval, panel) one representation of f is used on every slab, +the safe form (cusp.f3 with CplxOps) or the direct form (cusp.f3_bound_direct). Its success on every slab certifies +that it is analytic there, so it equals f on the ellipse.""" +import math +from fractions import Fraction as Fr +from ia import dn, up, add, sub, mul, div, scal, from_frac +import gauss as G +import cusp as CU +import cx as CX +import elem as EL + +R, C = CU.RealOps, CU.CplxOps +ZERO = (0.0, 0.0) +FAIL = (ValueError, ZeroDivisionError, OverflowError) + +S_EDGES = (0.0, 0.6, 0.9, 1.1, 1.25, 1.4, 1.55, 1.7, 1.9, 2.1, 2.35, 2.6, 3.0) +T_EDGES = (0.0, 0.4, 0.7, 0.9, 1.0, 1.1, 1.3, 1.6, 2.0, 2.6, 3.4, 4.5, 6.0, 8.0, 11.0, 15.0) +V_EDGES = (0.0, 0.35, 0.55, 0.65, 0.75, 0.85, 0.93, 1.0) +T_MAX = T_EDGES[-1] +S_MAX = S_EDGES[-1] +RHOS = (3.0, 2.2, 1.7, 1.4, 1.2) +SLABS = 8 +NMAX = 48 + + +def _sub(a, b, k): + return [(a + (b - a) * j / k, a + (b - a) * (j + 1) / k) for j in range(k)] + + +def _merge(pts, top): + out = [pts[0]] + for x in pts[1:]: + if x - out[-1] >= 0.02: + out.append(x) + elif x == top: + out[-1] = top + return tuple(out) + + +def s_edges(tb): + """s-panel edges for an outer box with t <= tb. For large t the integrand contains e^(-(23/240) s^8 v^8 t^8), and + for complex s with arg s > pi/16 that factor grows once |s| > 1.3/t. So extra edges c/tb are added below 2, + and the s-ellipses near the origin shrink with 1/t.""" + pts = sorted(set(S_EDGES) | {c / tb for c in (0.5, 0.8, 1.1, 1.5, 2.0, 2.8) if c / tb < 2.0}) + return _merge(pts, S_EDGES[-1]) + + +def v_edges(tb): + """v-panel edges for an outer box with t <= tb. For large t the integrand has branch points at |v| = 1/t, about + 0.38/t from the real axis, so extra edges c/tb are added below 0.9. Edges closer than 0.02 are merged, keeping the + first.""" + pts = sorted(set(V_EDGES) | {c / tb for c in (0.3, 0.5, 0.7, 1.0, 1.4, 2.0, 2.8, 4.0) if c / tb < 0.9}) + out = [pts[0]] + for x in pts[1:]: + if x - out[-1] >= 0.02 or x == 1.0: + if x == 1.0 and 1.0 - out[-1] < 0.02: + out[-1] = 1.0 + else: + out.append(x) + return tuple(out) + + +def _bnd(form, s, v, t): + if form == 'safe': + return CX.cabs_up(CU.f3(C, s, v, t)) + return CU.f3_bound_direct(s, v, t) + + +def _forms(t_lo): + """Representation order: the direct form first for t >= 1 (tighter there), the safe form first below.""" + return ('direct', 'safe') if t_lo >= 1.0 else ('safe', 'direct') + + +def _sup_over_slabs(slabs, mk, t_lo): + """max over the slabs of one representation (the first that succeeds on every slab); None if neither does.""" + for form in _forms(t_lo): + try: + m = 0.0 + for sl in slabs: + m = max(m, _bnd(form, *mk(sl))) + return m + except FAIL: + continue + return None + + +GOOD_N = 16 + + +def _search(h, tol, M_of): + """Try rho in RHOS (descending); M_of(rho) returns sup |f| on that ellipse or None. Stop at the first rho giving an + admissible n <= GOOD_N, else return the smallest admissible (n, err, rho) found. Deterministic.""" + best = None + for rho in RHOS: + M = M_of(rho) + if M is None: + continue + for n in range(4, NMAX + 1): + e = G.error_bound(n, rho, M, h) + if e <= tol: + if best is None or (n, e) < (best[0], best[1]): + best = (n, e, rho) + break + if best is not None and best[0] <= GOOD_N: + break + if best is None: + raise RuntimeError('no admissible Gauss-Legendre rule') + return best + + +def _v_rules(S, T, tol_v): + """Per v-panel: (n, err, rho) for the inner rule, with sup over s in S, t in T (real sub-intervals) and v on the + ellipse.""" + rules = [] + ve = v_edges(T[1]) + for (va, vb) in zip(ve, ve[1:]): + def M_of(rho, va=va, vb=vb): + slabs = G.ellipse_slabs(va, vb, rho, SLABS) + M = 0.0 + for si in _sub(S[0], S[1], 2): + for ti in _sub(T[0], T[1], 2): + m = _sup_over_slabs(slabs, lambda sl: ((si, ZERO), sl, (ti, ZERO)), ti[0]) + if m is None: + return None + M = max(M, m) + return M + rules.append(_search(0.5 * (vb - va), tol_v, M_of)) + return rules + + +def _int_v_bound(s_box, t_box): + """Upper bound of int_0^1 |f3(s, v, t)| dv over complex s_box, t_box (one of them real), by v-panels. For each + v sub-interval one representation is used (the first that succeeds).""" + tot = 0.0 + ve = v_edges(max(t_box[0][1], 1e-300)) + for (va, vb) in zip(ve, ve[1:]): + for vi in _sub(va, vb, 2): + m = None + for form in _forms(t_box[0][0]): + try: + m = _bnd(form, s_box, (vi, ZERO), t_box) + break + except FAIL: + continue + if m is None: + raise ValueError('no bound') + tot = up(tot + up(m * up(vi[1] - vi[0]))) + return tot + + +def _outer_rule(a, b, tol, other, which): + """Rule in the outer direction `which` ('s' or 't') on [a, b], with the other outer variable real over `other`.""" + def M_of(rho): + slabs = G.ellipse_slabs(a, b, rho, SLABS) + try: + M = 0.0 + for oi in _sub(other[0], other[1], 2): + for sl in slabs: + if which == 's': + M = max(M, _int_v_bound(sl, (oi, ZERO))) + else: + M = max(M, _int_v_bound((oi, ZERO), sl)) + return M + except FAIL: + return None + return _search(0.5 * (b - a), tol, M_of) + + +def box(S, T, tol): + """Certified enclosure of int_S int_T int_0^1 f3, with the rule data. Boxes whose crude enclosure is below + tol/10 use it directly.""" + cr = crude(S, T) + if cr is not None and cr[1] <= tol / 10.0: + return cr, {'crude': True, 'err': cr[1]} + lenS, lenT = S[1] - S[0], T[1] - T[0] + ve = v_edges(T[1]) + tol_v = tol / (3.0 * lenS * lenT * (len(ve) - 1)) + vr = _v_rules(S, T, tol_v) + ns, es, rs = _outer_rule(S[0], S[1], tol / (3.0 * lenT), T, 's') + nt, et, rt = _outer_rule(T[0], T[1], tol / (3.0 * lenS), S, 't') + xs, ws = G.panel(ns, S[0], S[1]) + xt, wt = G.panel(nt, T[0], T[1]) + vpanels = [] + for (va, vb), (nv, ev, rv) in zip(zip(ve, ve[1:]), vr): + vpanels.append(G.panel(nv, va, vb)) + Q = (0.0, 0.0) + for x1, w1 in zip(xs, ws): + for x2, w2 in zip(xt, wt): + inner = (0.0, 0.0) + for xv, wv in vpanels: + for x3, w3 in zip(xv, wv): + inner = add(inner, mul(w3, CU.f3(R, x1, x3, x2))) + Q = add(Q, mul(mul(w1, w2), inner)) + err_v = sum(e for (_, e, _) in vr) + err = up(up(es * lenT * (1 + 1e-12)) + up(et * lenS * (1 + 1e-12)) + up(err_v * lenS * lenT * (1 + 1e-12))) + return (dn(Q[0] - err), up(Q[1] + err)), {'ns': ns, 'nt': nt, 'nv': [r[0] for r in vr], + 'rho_s': rs, 'rho_t': rt, 'err': err} + + +def boxes(): + out = [] + for (ta, tb) in zip(T_EDGES, T_EDGES[1:]): + se = s_edges(tb) + out.extend(((sa, sb), (ta, tb)) for (sa, sb) in zip(se, se[1:])) + return out + + +def crude(S, T): + """0 <= int_box f3 <= sum over 2 x 2 x 4 sub-boxes of volume * sup f3, since f3 >= 0: the weight is nonnegative + and 1 - Re Phi >= 0 because |Phi| <= 1. Each sup is the upper end of one real interval evaluation of f3 over the + sub-box.""" + hi = 0.0 + try: + for si in _sub(S[0], S[1], 2): + for ti in _sub(T[0], T[1], 2): + for vi in _sub(0.0, 1.0, 4): + val = CU.f3(R, si, vi, ti) + vol = up(up(up(si[1] - si[0]) * up(ti[1] - ti[0])) * up(vi[1] - vi[0])) + hi = up(hi + up(max(0.0, val[1]) * vol)) + except FAIL: + return None + return (0.0, hi) + + +# ---------------------------------------------------------------- tails t > T_MAX and s > S_MAX +def tail_t(): + """int_0^oo ds int_0^1 dv int_(T_MAX)^oo dt f3, for real (s, v) (without the factor 64 of T3). + + Here f3 = wgt (1 - Re Phi) t^-8, and on the real axis |Phi| <= |1 - i t^4|^(-1/2) <= t^-2. So the t-integral is + wgt (T^-7/7 + theta T^-9/9) with |theta| <= 1. The s-integral is exact: + int_0^oo s^20 e^(-s^8 E) ds = Gamma(21/8)/(8 E^(21/8)), E(v) = (2 v^8 - v^4 + 2)/10 >= 3/16. + This leaves Gamma(21/8)/8 int_0^1 v^4 (1 - v^4) E(v)^(-21/8) dv. That integral is enclosed by interval evaluation + on 20000 cells; the result is multiplied by T^-7/7 < 4e-10.""" + import consts as K + n = 20000 + tot = (0.0, 0.0) + for j in range(n): + V = (j / n, (j + 1) / n) + v4 = mul(mul(V, V), mul(V, V)) + v4 = (max(0.0, v4[0]), v4[1]) + poly = mul(v4, sub((1.0, 1.0), v4)) + E = scal(0.1, add(sub(scal(2.0, mul(v4, v4)), v4), (2.0, 2.0))) + E = (max(E[0], 0.1874999), E[1]) + val = mul(poly, EL.exp_iv(mul(from_frac(Fr(-21, 8)), EL.log_iv(E)))) + tot = add(tot, mul(val, (dn(V[1] - V[0]), up(V[1] - V[0])))) + W = mul(div(K.gamma(Fr(21, 8)), (8.0, 8.0)), tot) + T = T_MAX + main = div((1.0, 1.0), (7.0 * T ** 7, 7.0 * T ** 7)) + e = up(1.0 / (9.0 * T ** 9) * (1 + 1e-12)) + return mul(W, (dn(main[0] - e), up(main[1] + e))) + + +def tail_s_bound(): + """Bound of int_(S_MAX)^oo ds int_0^1 dv int_0^oo dt f3 >= 0 (without the factor 64; NOTE.md section 4). + + Here 1 - Re Phi <= omega^2 E[Y^2]/2 and <= 2, so Gt <= min(E[Y^2]/(2 a1^2), 2 t^-8). With + E[Y^2] <= a1^6/50 + 3 a1^4 + 3 a1^2 and a1 = s^4 >= 81, int Gt dt <= a1^4 + 2/7 <= 2 s^16. Hence the s-integrand + is at most s^20 e^(-3 s^8/16) (1/4) 2 s^16 <= e^(-s^8/10)/2 for s >= 3 (s^36 <= e^(0.0875 s^8)), and + int_3^oo e^(-s^8/10)/2 ds <= e^(-656.1)/(1.6 * 3^7) < 1e-280.""" + if S_MAX != 3.0: + raise ValueError('the s-tail bound is proved for S_MAX = 3 only') + return 1e-280 diff --git a/frontiers/cusp_coefficient_certified_20261001/elem.py b/frontiers/cusp_coefficient_certified_20261001/elem.py new file mode 100644 index 00000000..99ad3332 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/elem.py @@ -0,0 +1,308 @@ +"""Real-interval elementary functions for CL-CU-CUSP-COEFFICIENT-20261001, built on ia.py (its conventions apply: an +interval is a pair (lo, hi) of floats, every operation is rounded outward). + +Point enclosures come from series with explicit remainder bounds; the only constants are pi and ln 2 of ia.py. +Interval extensions use the monotonicity stated with each function; sin and cos take their endpoint values and +-1 +at every critical point that may lie inside the interval. + + exp, log (t > 0), log1p (t > -1), atan, sin, cos, sinh, cosh, + E1(x) = (e^x - 1)/x, L(x) = log(1 + x)/x, A(x) = atan(x)/x, Sinc(x) = sin(x)/x (each = 1 at x = 0). + +E1 is increasing on R, L is decreasing on (-1, oo), A is even and decreasing in |x|, Sinc is even and decreasing on +[0, pi]. The four kernels are analytic, so the integrands of certificate.py never divide by a small quantity.""" +import math +from ia import (dn, up, add, sub, neg, mul, div, scal, divc, sqr, isqrt, hull, absmax, PI, LN2, exp_neg) + +ONE = (1.0, 1.0) +HALF_PI = scal(0.5, PI) + + +def _with_err(s, e): + """s widened by the nonnegative error bound e.""" + return (dn(s[0] - e), up(s[1] + e)) + + +# ---------------------------------------------------------------- exp +def exp_pt(t): + """Enclosure of e^t for a float t with t <= 700.""" + if t <= 0.0: + return exp_neg((-t, -t)) + if t > 700.0: + raise OverflowError('exp argument too large: %r' % t) + return div(ONE, exp_neg((t, t))) + + +def exp_iv(x): + return (exp_pt(x[0])[0], exp_pt(x[1])[1]) + + +# ---------------------------------------------------------------- log, log1p via atanh +def _atanh_series(u, K=26): + """atanh(u) = sum_k u^(2k+1)/(2k+1) for an interval u with |u| <= 1/3. The tail after K terms is at most + m^(2K+3)/((2K+3)(1 - m^2)) with m = max|u|; it is added with a safety factor 2.""" + m = absmax(u) + if m > 0.34: + raise ValueError('atanh series outside its range') + u2 = sqr(u) + term = u + s = u + for k in range(1, K + 1): + term = mul(term, u2) + s = add(s, divc(term, float(2 * k + 1))) + tail = 2.0 * m ** (2 * K + 3) / ((2 * K + 3) * (1.0 - m * m)) + return _with_err(s, tail) + + +def log_pt(t): + """Enclosure of log t for a float t > 0: t = m 2^e with m in [1/sqrt2, sqrt2), log m = 2 atanh((m-1)/(m+1)).""" + if not t > 0.0: + raise ValueError('log of a nonpositive number') + m, e = math.frexp(t) + if m < 0.7071067811865476: + m *= 2.0 + e -= 1 + M = (m, m) + u = div(sub(M, ONE), add(M, ONE)) + r = scal(2.0, _atanh_series(u)) + return add(r, mul((float(e), float(e)), LN2)) + + +def log_iv(x): + return (log_pt(x[0])[0], log_pt(x[1])[1]) + + +def log1p_pt(t): + """Enclosure of log(1 + t) for a float t > -1, with full relative accuracy near 0 (2 atanh(t/(2 + t))).""" + if not t > -1.0: + raise ValueError('log1p argument <= -1') + if abs(t) <= 0.5: + T = (t, t) + u = div(T, add((2.0, 2.0), T)) + return scal(2.0, _atanh_series(u)) + return log_iv((dn(1.0 + t), up(1.0 + t))) + + +def log1p_iv(x): + return (log1p_pt(x[0])[0], log1p_pt(x[1])[1]) + + +# ---------------------------------------------------------------- atan +def _atan_small(y, K=14): + """Alternating series for atan on an interval y within [0, 0.21]: the error after K terms is at most the next + term.""" + y2 = sqr(y) + term = y + s = y + for k in range(1, K + 1): + term = mul(term, y2) + t = divc(term, float(2 * k + 1)) + s = sub(s, t) if k % 2 else add(s, t) + m = absmax(y) + nxt = m ** (2 * K + 3) / (2 * K + 3) + return _with_err(s, 2.0 * nxt) + + +def atan_pt(t): + """Enclosure of atan t for a float t: odd symmetry, atan t = pi/2 - atan(1/t) for t > 1, and two halvings + atan x = 2 atan(x/(1 + sqrt(1 + x^2))) bringing the argument below tan(pi/16) < 0.199.""" + if t < 0.0: + return neg(atan_pt(-t)) + if t == 0.0: + return (0.0, 0.0) + if t > 1.0: + x = div(ONE, (t, t)) + return sub(HALF_PI, _atan_reduced(x)) + return _atan_reduced((t, t)) + + +def _atan_reduced(x): + y = x + for _ in range(2): + y = div(y, add(ONE, isqrt(add(ONE, sqr(y))))) + if y[1] > 0.21: + raise ValueError('atan reduction failed') + return scal(4.0, _atan_small(y)) + + +def atan_iv(x): + return (atan_pt(x[0])[0], atan_pt(x[1])[1]) + + +# ---------------------------------------------------------------- sin, cos +def _sincos_small(r, K=12): + """sin and cos on an interval r within [-0.8, 0.8] by their alternating Taylor series; the error after K terms is + at most the next term.""" + r2 = sqr(r) + s = r + c = ONE + ts = r + tc = ONE + for k in range(1, K + 1): + ts = divc(mul(ts, r2), float((2 * k) * (2 * k + 1))) + tc = divc(mul(tc, r2), float((2 * k - 1) * (2 * k))) + if k % 2: + s = sub(s, ts) + c = sub(c, tc) + else: + s = add(s, ts) + c = add(c, tc) + m = absmax(r) + es = m ** (2 * K + 3) / math.factorial(2 * K + 3) + ec = m ** (2 * K + 2) / math.factorial(2 * K + 2) + return _with_err(s, 2.0 * es), _with_err(c, 2.0 * ec) + + +def sincos_pt(t): + """Enclosures of (sin t, cos t) for a float t with |t| <= 1e6: t = k pi/2 + r with |r| <= pi/4 + 1e-9.""" + if abs(t) > 1e6: + raise ValueError('sin/cos argument too large') + k = round(t / 1.5707963267948966) + r = sub((t, t), mul((float(k), float(k)), HALF_PI)) + if absmax(r) > 0.8: + raise ValueError('sin/cos reduction failed') + s, c = _sincos_small(r) + q = k % 4 + if q == 0: + return s, c + if q == 1: + return c, neg(s) + if q == 2: + return neg(s), neg(c) + return neg(c), s + + +def _clip1(x): + return (max(-1.0, x[0]), min(1.0, x[1])) + + +def _crit_inside(x, offset): + """Integers k such that the enclosure of offset + k pi meets [x0, x1] (offset 0 or pi/2), with their parity.""" + out = [] + k0 = math.floor(x[0] / 3.141592653589793) - 2 + k1 = math.ceil(x[1] / 3.141592653589793) + 2 + for k in range(k0, k1 + 1): + c = add(offset, mul((float(k), float(k)), PI)) + if c[0] <= x[1] and c[1] >= x[0]: + out.append(k) + return out + + +def sin_iv(x): + """Range of sin over an interval: the endpoint values, together with +-1 at every pi/2 + k pi that may lie + inside (decided against the interval enclosure of pi).""" + if x[1] - x[0] >= 6.3: + return (-1.0, 1.0) + a, b = sincos_pt(x[0])[0], sincos_pt(x[1])[0] + lo, hi = min(a[0], b[0]), max(a[1], b[1]) + for k in _crit_inside(x, HALF_PI): + if k % 2 == 0: + hi = 1.0 + else: + lo = -1.0 + return _clip1((lo, hi)) + + +def cos_iv(x): + """Range of cos over an interval: the endpoint values, together with +-1 at every k pi that may lie inside.""" + if x[1] - x[0] >= 6.3: + return (-1.0, 1.0) + a, b = sincos_pt(x[0])[1], sincos_pt(x[1])[1] + lo, hi = min(a[0], b[0]), max(a[1], b[1]) + for k in _crit_inside(x, (0.0, 0.0)): + if k % 2 == 0: + hi = 1.0 + else: + lo = -1.0 + return _clip1((lo, hi)) + + +# ---------------------------------------------------------------- sinh, cosh (used only in complex bounds) +def cosh_iv(x): + """cosh is even and increasing in |x|.""" + lo = 0.0 if x[0] <= 0.0 <= x[1] else min(abs(x[0]), abs(x[1])) + hi = max(abs(x[0]), abs(x[1])) + def ch(t): + return scal(0.5, add(exp_pt(t), exp_pt(-t))) + return (ch(lo)[0], ch(hi)[1]) + + +def sinh_iv(x): + """sinh is odd and increasing.""" + def sh(t): + return scal(0.5, sub(exp_pt(t), exp_pt(-t))) + return (sh(x[0])[0], sh(x[1])[1]) + + +# ---------------------------------------------------------------- the four analytic kernels +def E1_pt(t): + """(e^t - 1)/t: series sum t^k/(k+1)! for |t| <= 1/2 (tail <= 2|t|^(K+1)/(K+2)!), else (e^t - 1)/t.""" + if abs(t) <= 0.5: + T = (t, t) + s = ONE + term = ONE + K = 18 + for k in range(1, K + 1): + term = divc(mul(term, T), float(k + 1)) + s = add(s, term) + err = 2.0 * abs(t) ** (K + 1) / math.factorial(K + 2) + return _with_err(s, err) + return div(sub(exp_pt(t), ONE), (t, t)) + + +def E1_iv(x): + """E1 is increasing.""" + return (E1_pt(x[0])[0], E1_pt(x[1])[1]) + + +def L_pt(t): + """log(1 + t)/t for t > -1 (1 at t = 0).""" + if t == 0.0: + return ONE + return div(log1p_pt(t), (t, t)) + + +def L_iv(x): + """L is decreasing on (-1, oo).""" + return (L_pt(x[1])[0], L_pt(x[0])[1]) + + +def A_pt(t): + """atan(t)/t (1 at t = 0).""" + if t == 0.0: + return ONE + return div(atan_pt(t), (t, t)) + + +def A_iv(x): + """A is even and decreasing in |x|.""" + lo = 0.0 if x[0] <= 0.0 <= x[1] else min(abs(x[0]), abs(x[1])) + hi = max(abs(x[0]), abs(x[1])) + return (A_pt(hi)[0], A_pt(lo)[1]) + + +def Sinc_pt(t): + """sin(t)/t: alternating series for |t| <= 1/2, else sin(t)/t.""" + if abs(t) <= 0.5: + T = (t, t) + t2 = sqr(T) + s = ONE + term = ONE + K = 10 + for k in range(1, K + 1): + term = divc(mul(term, t2), float((2 * k) * (2 * k + 1))) + s = sub(s, term) if k % 2 else add(s, term) + err = 2.0 * abs(t) ** (2 * K + 2) / math.factorial(2 * K + 3) + return _with_err(s, err) + return div(sincos_pt(t)[0], (t, t)) + + +def Sinc_iv(x): + """Sinc is even and decreasing on [0, pi]; outside [-3, 3] it is sin(x)/x, and |Sinc| <= 1 everywhere with + Sinc >= -0.2173 (its global minimum is -0.21723... at |x| = 4.4934...).""" + hi = max(abs(x[0]), abs(x[1])) + if hi <= 3.0: + lo = 0.0 if x[0] <= 0.0 <= x[1] else min(abs(x[0]), abs(x[1])) + return (Sinc_pt(hi)[0], Sinc_pt(lo)[1]) + if x[0] > 0.0 or x[1] < 0.0: + return div(sin_iv(x), x) + return (-0.2173, 1.0) diff --git a/frontiers/cusp_coefficient_certified_20261001/gauss.py b/frontiers/cusp_coefficient_certified_20261001/gauss.py new file mode 100644 index 00000000..71f0b570 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/gauss.py @@ -0,0 +1,147 @@ +"""Gauss-Legendre rules with rigorous node and weight enclosures, and the Bernstein-ellipse error bound used by +certificate.py (CL-CU-CUSP-COEFFICIENT-20261001). + +Nodes. P_n is evaluated exactly in rational arithmetic by (k+1) P_(k+1) = (2k+1) x P_k - k P_(k-1). A floating +Newton root is refined by exact Newton steps on the grid 2^-200 to a rational r, and accepted only when P_n has +strictly opposite signs at r - 2^-150 and r + 2^-150. The n enclosing intervals are +disjoint, so they contain all n roots, one each. The rule is symmetric, so only the roots in (0, 1) are verified; the +others are their negatives (and 0 itself when n is odd). + +Weights. At a root, (1 - x^2) P_n'(x) = n P_(n-1)(x), so w = 2 (1 - x^2)/(n P_(n-1)(x))^2. It is enclosed over the +node interval by its exact midpoint value plus a Lipschitz bound (_weight). + +Error bound (Lemma Q of NOTE.md). Let f be analytic in the open Bernstein ellipse E_rho of [-1, 1] and bounded by M +there. Then its Chebyshev coefficients satisfy |a_k| <= 2 M rho^(-k). The n-point rule is exact up to degree 2n - 1 +and is exact on odd T_k by symmetry. For even k, |int T_k| <= 2/(k^2 - 1) and |Q_n(T_k)| <= sum w = 2. Hence + |int f - Q_n f| <= 4 M (1 + 1/(4n^2 - 1)) rho^(-2n) / (1 - rho^(-2)), +and on [a, b] the right side is multiplied by h = (b - a)/2, for the ellipse c + h E_rho.""" +import math +from fractions import Fraction as Fr +from ia import dn, up, add, sub, mul, div, scal, divc, sqr, from_frac + +_CACHE = {} + + +def _legendre_exact(x, n): + """(P_(n-1)(x), P_n(x)) in exact rational arithmetic at a rational x, n >= 1.""" + p0 = Fr(1) + p1 = x + for k in range(1, n): + p0, p1 = p1, ((2 * k + 1) * x * p1 - k * p0) / (k + 1) + return p0, p1 + + +def _legendre_float(x, n): + p0, p1 = 1.0, x + for k in range(1, n): + p0, p1 = p1, ((2 * k + 1) * x * p1 - k * p0) / (k + 1) + return p0, p1 + + +def _weight(lo, hi, n): + """Enclosure of w(x) = 2 (1 - x^2)/(n P_(n-1)(x))^2 over the rational interval [lo, hi]: exact value at the + midpoint plus a Lipschitz bound. On [-1, 1], |P_(n-1)'| <= P_(n-1)'(1) = n (n - 1)/2, so |P_(n-1)| >= m on + [lo, hi] with m = min(|P_(n-1)(lo)|, |P_(n-1)(hi)|) - n (n - 1)/2 (hi - lo), and |w'| <= (4/n^2)(1/m^2 + (n (n - 1)/2)/m^3).""" + D1 = Fr(n * (n - 1), 2) + M0 = (lo + hi) / 2 + plo = abs(_legendre_exact(lo, n)[0]) + phi = abs(_legendre_exact(hi, n)[0]) + m = min(plo, phi) - D1 * (hi - lo) + if m <= 0: + raise RuntimeError('Gauss-Legendre weight: P_(n-1) not bounded away from 0') + pm = _legendre_exact(M0, n)[0] + wm = 2 * (1 - M0 * M0) / (n * n * pm * pm) + lip = Fr(4, n * n) * (1 / (m * m) + D1 / (m * m * m)) + err = lip * (hi - lo) / 2 + return (from_frac(wm - err)[0], from_frac(wm + err)[1]) + + +_GRID = Fr(1, 2 ** 200) +_EPS = Fr(1, 2 ** 150) + + +def _snap(r): + """Round the rational r to the grid 2^-200 (keeps the exact Newton iterates small).""" + return Fr(round(r / _GRID)) * _GRID + + +def rule(n): + """Nodes and weights on [-1, 1] as lists of intervals, ascending. Each root is refined by exact rational Newton + steps on the grid 2^-200, then certified by a strict sign change of P_n at r -+ 2^-150 (exact arithmetic); the + node interval is the outward float hull of [r - 2^-150, r + 2^-150].""" + if n in _CACHE: + return _CACHE[n] + pos = [] + for i in range(1, n // 2 + 1): + x = math.cos(math.pi * (i - 0.25) / (n + 0.5)) + for _ in range(100): + pm1, pn = _legendre_float(x, n) + dp = n * (x * pn - pm1) / (x * x - 1.0) + dx = pn / dp + x -= dx + if abs(dx) < 1e-17: + break + r = _snap(Fr(x)) + for _ in range(4): + pm1, pn = _legendre_exact(r, n) + dp = n * (r * pn - pm1) / (r * r - 1) + r = _snap(r - pn / dp) + lo, hi = r - _EPS, r + _EPS + if not _legendre_exact(lo, n)[1] * _legendre_exact(hi, n)[1] < 0: + raise RuntimeError('Gauss-Legendre node %d of %d not verified' % (i, n)) + pos.append(((from_frac(lo)[0], from_frac(hi)[1]), _weight(lo, hi, n))) + pos.sort(key=lambda e: e[0][0]) + for j in range(1, len(pos)): + if not pos[j - 1][0][1] < pos[j][0][0]: + raise RuntimeError('Gauss-Legendre node enclosures overlap') + if pos and not pos[0][0][0] > 0.0: + raise RuntimeError('Gauss-Legendre positive node enclosure meets 0') + nodes, weights = [], [] + for X, w in reversed(pos): + nodes.append((-X[1], -X[0])) + weights.append(w) + if n % 2: + pm1 = _legendre_exact(Fr(0), n)[0] + w0 = Fr(2, n * n) / (pm1 * pm1) + nodes.append((0.0, 0.0)) + weights.append(from_frac(w0)) + for X, w in pos: + nodes.append(X) + weights.append(w) + _CACHE[n] = (nodes, weights) + return nodes, weights + + +def panel(n, a, b): + """Rule on [a, b] (floats): nodes c + h x, weights h w, as intervals.""" + nodes, weights = rule(n) + A, B = (a, a), (b, b) + c = scal(0.5, add(A, B)) + h = scal(0.5, sub(B, A)) + return [add(c, mul(h, x)) for x in nodes], [mul(h, w) for w in weights] + + +def error_bound(n, rho, M, h): + """Upper bound of |int_[a,b] f - Q_n f| for |f| <= M on the ellipse c + h E_rho (h = (b - a)/2, rho > 1).""" + num = 4.0 * M * (1.0 + 1.0 / (4 * n * n - 1)) * rho ** (-2 * n) + return up(up(num / (1.0 - rho ** -2)) * h * (1.0 + 1e-12)) + + +def ellipse_slabs(a, b, rho, m): + """m complex boxes (re, im) covering the closed ellipse c + h E_rho of [a, b]: the real extent + c +- h (rho + 1/rho)/2 is cut into m slabs, and over each slab the imaginary half-height is the ellipse's at the slab + point nearest to c.""" + c = 0.5 * (a + b) + h = 0.5 * (b - a) + ax = h * (rho + 1.0 / rho) / 2.0 * (1.0 + 1e-12) + by = h * (rho - 1.0 / rho) / 2.0 * (1.0 + 1e-12) + out = [] + for j in range(m): + x0 = c - ax + 2.0 * ax * j / m + x1 = c - ax + 2.0 * ax * (j + 1) / m + x0, x1 = dn(x0), up(x1) + near = 0.0 if x0 <= c <= x1 else min(abs(x0 - c), abs(x1 - c)) + frac = max(0.0, 1.0 - (near / ax) ** 2) + yh = up(by * math.sqrt(frac) * (1.0 + 1e-12)) + out.append(((x0, x1), (-yh, yh))) + return out diff --git a/frontiers/cusp_coefficient_certified_20261001/ia.py b/frontiers/cusp_coefficient_certified_20261001/ia.py new file mode 100644 index 00000000..cf4bbfa1 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_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/cusp_coefficient_certified_20261001/side24.py b/frontiers/cusp_coefficient_certified_20261001/side24.py new file mode 100644 index 00000000..31da88a0 --- /dev/null +++ b/frontiers/cusp_coefficient_certified_20261001/side24.py @@ -0,0 +1,225 @@ +"""Lemma S: (CU.2) for the [P] torus field of side L >= 24, in particular SIDE24 (NOTE.md section 3b). Exact rational +arithmetic only. + +The [P] field of side L has covariance K_L(z) = sum_n phi(z + L n) / S_L, phi(x) = exp(-|x|^2/2), S_L = sum_n phi(L n) +(the normalization of coefficients/side24_v1; the unnormalized periodization is covered too). In an orthonormal frame +(u, Theta) the jets entering (CU.2) are the 2d + 2 pinned coordinates (f, d_u f, d_u^2 f, d_u^3 f, grad_Theta f, +grad_Theta d_u f) and the free ones (d_u^4 f, gamma = grad_Theta d_u^2 f, A = D_Theta^2 f). Every covariance entry is +(-1)^|b| D^(|a|+|b|) K(0) contracted with frame vectors, of total order at most 8. + +S.1 (image bound). For unit vectors v_1..v_q, D^q phi(x)[v] = phi(x) sum over partial matchings of products of + - and -, so |D^q phi(x)[v]| <= phi(x) sum_j q! |x|^(q-2j)/(2^j j! (q-2j)!). For q <= 8 this is at most + T_8 |x|^8 phi(x) when |x| >= 1 (T_8 = 764 matchings of 8 points), and at most 7!! = 105 at x = 0. For d <= 3 the shell + max|n_i| = j has at most 27 j^3 points with |n|^2 <= 3 j^2, and consecutive shell terms shrink by at least 1/2, so + sum_(n != 0) |n|^8 e^(-288 |n|^2) <= 2 . 27 . 81 e^(-288), S_24 - 1 <= 2 . 27 e^(-288), e^(-288) < 10^-125. + L^8 e^(-L^2 |n|^2/2) decreases in L for L |n| >= 3, so the bound at L = 24 holds for every L >= 24. Hence every entry of + C_L - C_0 has absolute value at most E = (764 . 24^8 . 4374 + 105 . 54) 10^-125. +S.2 (sandwich). C_0 - mu I is positive definite (exact LDL^T), so |x^T (C_L - C_0) x| <= N E |x|^2 <= eps x^T C_0 x with + eps = N E / mu, i.e. (1 - eps) C_0 <= C_L <= (1 + eps) C_0 in every frame. +S.3 (transfer). On the slice where 2d + 1 coordinates vanish, the (CU.2) integrand is a nonnegative function of the free + jets, homogeneous of degree 2d - 1/4, so the slice integral at a C_0 is a^(-5/8) times the one at C_0. The density + sandwich [(1 - eps)/(1 + eps)]^(N/2) phi_((1-eps) C_0) <= phi_(C_L) <= [(1 + eps)/(1 - eps)]^(N/2) phi_((1+eps) C_0) + then gives c1^(L)/c1^(ref) in [rho^-M, rho^M], rho = (1 + eps)/(1 - eps), M = N/2 + 1. Since log rho <= 4 eps and + e^x <= 1 + 2x on [0, 1], this is within eta = 8 M eps of 1.""" +import json +import os +from fractions import Fraction as Fr +from itertools import combinations + +L_MIN = 24 +IMAGE_ORDER = 8 # S.1 covers contractions of total order <= IMAGE_ORDER +MU = Fr(1, 4) # certified lower bound of the smallest eigenvalue of C_0 (S.2) + +# c_(d,24) of coefficients/side24_v1/ENCLOSURE.json ("dimensions"), consumed (SOURCE_MAP.json [SIDE24-E]). +C_24 = {2: ('0.07340691930603427103', '0.07340691930603427104'), + 3: ('0.04177593184059834334', '0.04177593184059834335')} + + +def _dfact(n): + """(n - 1)!! for even n >= 0: the number of perfect matchings of n points.""" + r = 1 + for k in range(n - 1, 0, -2): + r *= k + return r + + +def dK0(g): + """D^g phi(0) for the multi-index g: prod (-1)^(g_i/2) (g_i - 1)!! if every g_i is even, else 0.""" + if any(x % 2 for x in g): + return 0 + r = 1 + for x in g: + r *= (-1) ** (x // 2) * _dfact(x) + return r + + +def jets(d): + """The jets of (CU.2) in the frame (u = e_1, Theta = e_2..e_d): name, multi-index, role.""" + out = [] + for k in range(5): + out.append(('d_u^%d f' % k, (k,) + (0,) * (d - 1), 'b' if k == 0 else ('pin' if k < 4 else 'free'))) + for j in range(1, d): + for k in range(3): + g = [0] * d + g[0], g[j] = k, 1 + out.append(('d_y%d d_u^%d f' % (j, k), tuple(g), 'pin' if k < 2 else 'free')) + for j in range(1, d): + for k in range(j, d): + g = [0] * d + g[j] += 1 + g[k] += 1 + out.append(('d_y%d d_y%d f' % (j, k), tuple(g), 'free')) + return out + + +def cov0(d): + """Exact covariance of the jets for phi: Cov(D^a f, D^b f) = (-1)^|b| D^(a+b) phi(0).""" + J = jets(d) + return [[Fr((-1) ** sum(b) * dK0(tuple(x + y for x, y in zip(a, b)))) for _, b, _ in J] for _, a, _ in J] + + +def max_order(d): + J = jets(d) + return max(sum(a) + sum(b) for _, a, _ in J for _, b, _ in J) + + +def ldl_positive(M): + """True iff the symmetric rational matrix M is positive definite (exact LDL^T: every pivot > 0).""" + n = len(M) + A = [row[:] for row in M] + for k in range(n): + if A[k][k] <= 0: + return False + for i in range(k + 1, n): + f = A[i][k] / A[k][k] + if f: + for j in range(k, n): + A[i][j] -= f * A[k][j] + return True + + +def matchings(q): + """Number of partial matchings (involutions) of q points, by enumeration: T_q.""" + def count(pts): + if not pts: + return 1 + first, rest = pts[0], pts[1:] + tot = count(rest) + for i in range(len(rest)): + tot += count(rest[:i] + rest[i + 1:]) + return tot + return count(tuple(range(q))) + + +def exp_lower(x, K=20): + """A lower bound of e^x for rational x > 0: the positive Taylor sum through x^K/K!.""" + s, t = Fr(0), Fr(1) + for k in range(K + 1): + s += t + t = t * x / (k + 1) + return s + + +def image_bound(): + """E of S.1 (exact rational) after checking its numeric premises.""" + if not exp_lower(Fr(288, 125)) > 10: # e^(-288) < 10^-125 + raise ValueError('e^(288/125) > 10 not certified') + e288 = Fr(1, 10 ** 125) + if not 2 ** 11 * e288 ** 3 < Fr(1, 2) or not 2 ** 3 * e288 ** 3 < Fr(1, 2): + raise ValueError('shell ratio bound fails') # e^(-864) < 10^-375 + T = matchings(IMAGE_ORDER) + dfmax = max(_dfact(q) for q in range(0, IMAGE_ORDER + 1, 2)) + lattice = 2 * 27 * 3 ** (IMAGE_ORDER // 2) # sum |n|^q e^(-288|n|^2) <= 2 . 27 . 3^(q/2) e^(-288) + norm = 2 * 27 # S_24 - 1 <= 2 . 27 e^(-288) + return (T * L_MIN ** IMAGE_ORDER * lattice + dfmax * norm) * e288 + + +def transfer(d): + """S.2 and S.3 for dimension d: the exact constants and eta.""" + C0 = cov0(d) + N = len(C0) + if max_order(d) > IMAGE_ORDER: + raise ValueError('covariance entries of order %d exceed the image bound order %d' % (max_order(d), IMAGE_ORDER)) + if not ldl_positive([[C0[i][j] - (MU if i == j else 0) for j in range(N)] for i in range(N)]): + raise ValueError('C_0 - mu I is not positive definite') + E = image_bound() + eps = N * E / MU + M = Fr(N, 2) + 1 + if not 4 * M * eps <= 1: + raise ValueError('eps too large for the bound e^x <= 1 + 2x') + eta = 8 * M * eps + return {'N': N, 'E': E, 'mu': MU, 'eps': eps, 'M': M, 'eta': eta} + + +def lemma_r_check(d): + """Exact regression from C_0: the conditional law of the free jets given the pins, against Lemma R.1, and the pin + density constants det C_P = 12, (C_P^-1)_(f,f) = 3/2. Returns a list of failures.""" + J = jets(d) + C0 = cov0(d) + P = [i for i, j in enumerate(J) if j[2] in ('b', 'pin')] + F = [i for i, j in enumerate(J) if j[2] == 'free'] + CP = [[C0[i][j] for j in P] for i in P] + n = len(P) + # exact inverse of C_P by Gauss-Jordan + A = [CP[i][:] + [Fr(int(i == j)) for j in range(n)] for i in range(n)] + det = Fr(1) + for k in range(n): + piv = next(r for r in range(k, n) if A[r][k] != 0) + if piv != k: + A[k], A[piv] = A[piv], A[k] + det = -det + det *= A[k][k] + p = A[k][k] + A[k] = [x / p for x in A[k]] + for r in range(n): + if r != k and A[r][k] != 0: + f = A[r][k] + A[r] = [x - f * y for x, y in zip(A[r], A[k])] + inv = [row[n:] for row in A] + fails = [] + if det != 12: + fails.append('d = %d: det C_P = %s, expected 12' % (d, det)) + if inv[0][0] != Fr(3, 2): + fails.append('d = %d: (C_P^-1)_ff = %s, expected 3/2' % (d, inv[0][0])) + # regression coefficients on the b coordinate (pins at 0): mean = Cov(W, P) C_P^-1 e_b b; covariance = Schur complement + CFP = [[C0[i][j] for j in P] for i in F] + mean_b = [sum(CFP[r][k] * inv[k][0] for k in range(n)) for r in range(len(F))] + schur = [[C0[F[r]][F[s]] - sum(CFP[r][k] * inv[k][l] * CFP[s][l] for k in range(n) for l in range(n)) + for s in range(len(F))] for r in range(len(F))] + names = [J[i][0] for i in F] + m = d - 1 + exp_mean, exp_cov = {}, {} + exp_mean['d_u^4 f'] = -3 + exp_cov[('d_u^4 f', 'd_u^4 f')] = 24 + for j in range(1, d): + exp_mean['d_y%d d_u^2 f' % j] = 0 + exp_cov[('d_y%d d_u^2 f' % j, 'd_y%d d_u^2 f' % j)] = 2 + exp_mean['d_y%d d_y%d f' % (j, j)] = -1 + exp_cov[('d_y%d d_y%d f' % (j, j), 'd_y%d d_y%d f' % (j, j))] = 2 + for k in range(j + 1, d): + exp_mean['d_y%d d_y%d f' % (j, k)] = 0 + exp_cov[('d_y%d d_y%d f' % (j, k), 'd_y%d d_y%d f' % (j, k))] = 1 + for r, a in enumerate(names): + if mean_b[r] != exp_mean[a]: + fails.append('d = %d: E[%s | b] = %s b, expected %s b' % (d, a, mean_b[r], exp_mean[a])) + for s, b in enumerate(names): + want = exp_cov.get((a, b), exp_cov.get((b, a), 0)) + if schur[r][s] != want: + fails.append('d = %d: Cov(%s, %s | pins) = %s, expected %s' % (d, a, b, schur[r][s], want)) + if m != len([1 for a in names if a.endswith('d_u^2 f')]): + fails.append('d = %d: gamma has the wrong dimension' % d) + return fails + + +def enclosure_check(here): + """Compare C_24 with the repository copy of coefficients/side24_v1/ENCLOSURE.json when it is reachable.""" + path = os.path.join(here, '..', '..', 'coefficients', 'side24_v1', 'ENCLOSURE.json') + if not os.path.exists(path): + return None + with open(path) as fh: + dims = json.load(fh)['dimensions'] + bad = [] + for d, (lo, hi) in C_24.items(): + if (dims[str(d)]['lower'], dims[str(d)]['upper']) != (lo, hi): + bad.append('c_(%d,24) transcription differs from ENCLOSURE.json' % d) + return bad