diff --git a/.github/workflows/side24-d4-coefficient.yml b/.github/workflows/side24-d4-coefficient.yml new file mode 100644 index 000000000..e5d25415c --- /dev/null +++ b/.github/workflows/side24-d4-coefficient.yml @@ -0,0 +1,90 @@ +name: SIDE24 d4 coefficient + +'on': + pull_request: + paths: + - 'coefficients/side24_d4_20260930/**' + - '.github/workflows/side24-d4-coefficient.yml' + - 'imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md' + - 'coefficients/side24_v1/PROOF.md' + - 'coefficients/side24_v1/ENCLOSURE.json' + workflow_dispatch: + +permissions: + contents: read + +jobs: + certify-and-reject-mutants: + runs-on: ubuntu-latest + timeout-minutes: 30 + steps: + - uses: actions/checkout@11d5960a326750d5838078e36cf38b85af677262 + with: + persist-credentials: false + - uses: actions/setup-python@a26af69be951a213d495a4c3e4e4022e16d87065 + with: + python-version: '3.11.16' + - name: Verify sources and pins, replay both modes, reject mutants + env: + PYTHONDONTWRITEBYTECODE: '1' + run: | + set -euo pipefail + python -B -S - <<'PY' + import concurrent.futures, hashlib, json, pathlib, subprocess, sys + root = pathlib.Path('coefficients/side24_d4_20260930').resolve() + manifest_doc = json.loads((root / 'SOURCE_FILES.json').read_text()) + manifest = manifest_doc['files'] + present = sorted(p.name for p in root.iterdir()) + expected_names = sorted([e['path'] for e in manifest] + ['SOURCE_FILES.json']) + if present != expected_names: + 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() + if len(data) != entry['bytes'] or hashlib.sha256(data).hexdigest() != entry['sha256']: + raise ValueError('source identity mismatch: ' + entry['path']) + for pin in manifest_doc['pins_on_main']: + if not pin.get('workflow_checked'): + continue + p = pathlib.Path(pin['path']) + if p.is_symlink() or not p.is_file(): + raise ValueError('pinned source missing: ' + pin['path']) + data = p.read_bytes() + blob = hashlib.sha1(b'blob %d\0' % len(data) + data).hexdigest() + if blob != pin['git_blob']: + raise ValueError('pinned source drifted: ' + pin['path'] + ' ' + blob) + # the SIDE24 d = 2, 3 enclosures quoted inside the script must be the ones on main + enc = json.loads(pathlib.Path('coefficients/side24_v1/ENCLOSURE.json').read_text())['dimensions'] + src = (root / 'cone_moment_d4.py').read_text() + for d in ('2', '3'): + if ('"%s", "%s"' % (enc[d]['lower'], enc[d]['upper'])) not in src: + raise ValueError('SIDE24 d = ' + d + ' enclosure quoted in the script differs from ENCLOSURE.json') + script = str(root / 'cone_moment_d4.py') + expected = (root / 'RESULTS.json').read_bytes() + mutants = ['shift-variance', 'vandermonde', 'erfc-branch', 'remainder-dropped', 'mehta', 'sphere-area', 'image-shells', 'sandwich-exponent', 'recurrence'] + def mutant_rc(args): + flags, m = args + r = subprocess.run([sys.executable, *flags, script, '--mutant', m], + capture_output=True, timeout=1500) + try: + ok = r.returncode == 1 and json.loads(r.stdout)['passed'] is False + except (ValueError, KeyError, TypeError): + ok = False + return m, flags, ok + for flags in (['-B', '-S'], ['-B', '-O', '-S']): + out = subprocess.run([sys.executable, *flags, script], + check=True, capture_output=True, timeout=1500).stdout + if out != expected: + raise RuntimeError('output differs from RESULTS.json under ' + ' '.join(flags)) + doc = json.loads(out) + if not (doc['passed'] and len(doc['checks']) == 13 and all(doc['checks'].values())): + raise RuntimeError('certification did not pass under ' + ' '.join(flags)) + with concurrent.futures.ThreadPoolExecutor(max_workers=4) as pool: + for m, fl, ok in pool.map(mutant_rc, [(flags, m) for m in mutants]): + if not ok: + raise RuntimeError('mutant ' + m + ' not rejected by its rules under ' + ' '.join(fl)) + print(json.dumps({'scientific_effect': 'NONE', 'mathematical_acceptance': False, 'certified_enclosures': True, 'side24_v1': 'unchanged'})) + PY + git diff --exit-code diff --git a/coefficients/side24_d4_20260930/NOTE.md b/coefficients/side24_d4_20260930/NOTE.md new file mode 100644 index 000000000..e00de4665 --- /dev/null +++ b/coefficients/side24_d4_20260930/NOTE.md @@ -0,0 +1,247 @@ +# SIDE24 in dimension four: the birth-integrated cone moment `D_3`, certified, with `c_{4,ref}` and `c_{4,24}` + +**Object:** `CL-SIDE24-D4-COEFFICIENT-20260930-v1`. **Author lane:** Anthropic / Claude (Claude Code +`session_017Mi3hxjaxV45x6zo6o1ee3`). **Kind:** certified numerical constant with its exact reduction: interval arithmetic +over `decimal` (60 digits, directed rounding), exact integer roots, series with explicit remainders, and a one-dimensional +Taylor quadrature whose remainder is bounded by Cauchy's estimate. **Scientific effect:** NONE. No register, catalog, +GRAPH or STATUS change; `coefficients/side24_v1` is unchanged; catalog entry C8 stays OPEN. Scope claimed on Math-#197 +(comment 5916664092). Same GitHub account as every lane: zero organizational-independence credit. + +## 0. Statement + +`coefficients/side24_v1` evaluates the coefficient (15.2) of [LP] for the periodized Gaussian kernel of side `L = 24` in +`d = 2, 3`. Its reference law (section 1 there) is, in dimension `d = m + 1`: the transverse Hessian given `V = 0` is +`A = Q + sqrt(2/3) Z I_m`, `Q` the `m x m` GOE with density proportional to `exp(-tr Q^2/4)` (diagonal variance `2`, +off-diagonal `1`) and `Z` an independent standard normal, and the cone moment of [LP] (15.1) is + + D_m = E[ det(A)^2 1{A negative definite} ], D_1 = 4/3, D_2 = 29/6 - sqrt6. + +**Certified.** `D_3 = 5.3231802688898968949231987396872606357843885294...` (enclosure width `6.5e-47`), whence, by [LP] +(15.2) with `|S^3| = 2 pi^2`, `p_G(0) p_V(0) = (2 pi)^-4/sqrt3`, `tau^2 = 6`, + + c_{4,ref} = Gamma(7/6) (3/2)^(1/3) D_3 / (8 sqrt3 pi^(5/2)) = 0.023321666002952835094521194952856885153101372363... + +and for the torus kernel `K_24` in `d = 4`, by the covariance sandwich of SIDE24 sections 3-4 with the `d = 4` constants, + + |c_{4,24} / c_{4,ref} - 1| <= delta = 1.5702e-108, c_{4,24} = 0.023321666002952835094521194952856885153101372363... + +The same code path at `d = 2, 3` returns enclosures inside SIDE24's intervals for `c_{2,24}` and `c_{3,24}` (rule +SIDE24_CONSISTENT), and its `m = 2` cone moment encloses `29/6 - sqrt6` (rule D2_EXACT), which tests every step of the +reduction that `m = 3` shares. + +## 1. The cone moment as a shifted GOE integral + +Let `lambda_1, ..., lambda_m` be the eigenvalues of `Q`, with joint density `|Delta(lambda)| exp(-|lambda|^2/4) / Z_m` on +`R^m` (unordered), `Z_m = int_{R^m} prod_{i 0`), so +`|Delta| = p q (p+q)` and `|lambda|^2/4 = (3/4)(t - s/3)^2 + (p^2 + pq + q^2)/6` with `s = 2p + q`; the `t`-integral is +`2 sqrt(pi/3)`, and in `(s, q)` the rest is `s^2/24 + q^2/8` with `p q (p+q) = q (s^2 - q^2)/4`. Hence +`Z_3 = 6 . 2 sqrt(pi/3) . (1/8) int_0^inf e^{-s^2/24} [s^2 J_1(s) - J_3(s)] ds = (3/2) sqrt(pi/3) int_0^inf e^{-s^2/24} (4 s^2 - 32 + 32 e^{-s^2/8}) ds`, +and `int_0^inf s^2 e^{-s^2/24} ds = 12 sqrt6 sqrt pi`, `int_0^inf e^{-s^2/24} ds = sqrt6 sqrt pi`, `int_0^inf e^{-s^2/6} ds = (1/2) sqrt6 sqrt pi` +give `(3/2) sqrt(pi/3) . 32 sqrt6 sqrt pi = 48 sqrt2 pi`. Rule MEHTA_Z3 certifies this integral by the section 3 quadrature and +requires it to enclose both the Mehta value and `48 sqrt2 pi`. + +## 2. `m = 3`: exact reduction to one dimension + +On the ordered sector `mu_1 < mu_2 < mu_3 < 0` put `mu_3 = -a`, `mu_2 = -a - p`, `mu_1 = -a - p - q` (`a, p, q > 0`, +Jacobian `1`). Then `prod mu_i^2 = a^2 (a+p)^2 (a+p+q)^2`, `|Delta| = p q (p+q)`, and + + |mu|^2/4 - S^2/24 = Q(a, p, q) = (3/8) a^2 + (1/3) p^2 + (5/24) q^2 + (1/2) a p + (1/4) a q + (1/3) p q, + +so that `D_3 = (1/sqrt2) (6/Z_3) I_3 = I_3 / (16 pi)`, `I_3 = int_{(0,inf)^3} a^2 (a+p)^2 (a+p+q)^2 p q (p+q) e^{-Q} da dp dq`. + +**The `a`-integral.** With `alpha = 3/8`, `beta = p/2 + q/4`, `G_n(beta) = int_0^inf a^n e^{-alpha a^2 - beta a} da`: + + G_0 = (1/2) sqrt(pi/alpha) e^{beta^2/(4 alpha)} erfc(beta/(2 sqrt alpha)), G_1 = (1 - beta G_0)/(2 alpha), + G_{n+1} = (n G_{n-1} - beta G_n)/(2 alpha) (integrate (a^n e^{...})' over (0, inf)), + +so `G_n = P_n(beta) G_0 + R_n(beta)` with polynomials `P_n` (degree `n`) and `R_n` (degree `n-1`) of rational +coefficients. Expanding `a^2 (a+p)^2 (a+p+q)^2 = sum_n w_n(p,q) a^n` (`n = 2..6`) and multiplying by `p q (p+q)`, + + I_3 = int_{(0,inf)^2} [ (1/2) sqrt(8 pi/3) P_E(p,q) e^{-(p^2 + pq + q^2)/6} erfc((2p+q)/(2 sqrt6)) + P_R(p,q) e^{-(p^2/3 + 5q^2/24 + pq/3)} ] dp dq, + +because `beta^2/(4 alpha) = p^2/6 + pq/6 + q^2/24` cancels against the `(p,q)` part of `Q` to leave `-(p^2 + pq + q^2)/6`, +and `beta/(2 sqrt alpha) = (2p + q)/(2 sqrt6)`; `P_E = sum_n w_n P_n(beta) . pq(p+q)` (degree `9`), `P_R = sum_n w_n R_n(beta) . pq(p+q)` +(degree `8`), all with rational coefficients computed exactly by the script. + +**Decoupling.** With `s = 2p + q` and `q` (so `p = (s - q)/2`, `dp dq = ds dq/2`, region `0 < q < s`): + + p^2 + pq + q^2 = s^2/4 + 3 q^2/4, p^2/3 + 5q^2/24 + pq/3 = s^2/12 + q^2/8, + +so both exponents split, `-(s^2/24 + q^2/8)` and `-(s^2/12 + q^2/8)`, and the `erfc` depends on `s` alone (the script +verifies this structure symbolically before evaluating anything). Writing `P_E((s-q)/2, q) = sum_k e_k(s) q^k`, +`P_R((s-q)/2, q) = sum_k r_k(s) q^k` (`k <= 9`) and `J_k(s) = int_0^s q^k e^{-q^2/8} dq`, + + I_3 = (1/2) int_0^inf [ (1/2) sqrt(8 pi/3) e^{-s^2/24} erfc(s/(2 sqrt6)) sum_k e_k(s) J_k(s) + e^{-s^2/12} sum_k r_k(s) J_k(s) ] ds, (2.1) + J_0 = sqrt(2 pi) erf(s/(2 sqrt2)), J_1 = 4 (1 - e^{-s^2/8}), J_{k+2} = -4 s^{k+1} e^{-s^2/8} + 4(k+1) J_k. + +Every factor in (2.1) is an entire function of `s`: polynomials, `exp(-gamma s^2)`, `erf`/`erfc` of a multiple of `s`, and +the `J_k`. + +## 3. Certified quadrature (Taylor expansion with Cauchy remainders) + +On each cell `[s_0 - h/2, s_0 + h/2]` (`h = 1/2`, `s_0` rational, `0 < s < T = 60`) the integrand `F` is expanded to order +`K = 40` at `s_0`: the Taylor coefficients of `exp(-gamma s^2)` come from `(n+1) y_{n+1} = -2 gamma (s_0 y_n + y_{n-1})` +with `y_0 = exp(-gamma s_0^2)` (an interval), those of `erf(kappa s)`/`erfc(kappa s)` from the Gaussian's +(`(2 kappa/sqrt pi) exp(-kappa^2 s^2)`) by division by `n + 1`, with `erf(kappa s_0)` from the positive series, those of the +`J_k` from their recurrence in Taylor arithmetic, and polynomials exactly; products are Cauchy products of interval +coefficient lists. The cell integral of the degree-`K` polynomial is exact (odd powers vanish). For the remainder: `F` is entire, so its +Taylor series at `s_0`, `F(s_0 + x) = sum_n c_n x^n`, converges for every `x`, and Cauchy's estimate at the centre gives +`|c_n| <= M_F(rho)/rho^n` for every `n`, where `M_F(rho)` bounds `|F(z)|` on `|z - s_0| <= rho` (`rho = 6`); by the +triangle inequality `M_F` is the corresponding sum of products of factor bounds, each elementary on the disc: + + |exp(-gamma z^2)| <= exp(gamma rho^2) exp(-gamma (|s_0| - rho)_+^2), + |erf(kappa z)| <= (2 kappa/sqrt pi) |z| exp(kappa^2 rho^2), |erfc(kappa z)| <= 1 + that, + |J_k(z)| = |int_0^1 (tz)^k e^{-t^2 z^2/8} z dt| <= |z|^(k+1) exp(rho^2/8)/(k+1), |poly(z)| <= sum |c_j| |z|^j, + +with `|z| <= |s_0| + rho` (for `exp`: `|e^{-gamma z^2}| = e^{-gamma(x^2 - y^2)}`, `|y| <= rho`, `|x| >= (|s_0| - rho)_+`; for +`erf` and `J_k`: integrate along the segment from `0` to `z`, on which `|e^{-t^2}| <= e^{(Im t)^2}`). The omitted tail +of the series integrates to at most + + int_{-h/2}^{h/2} sum_{n > K} |c_n| |x|^n dx <= M_F(rho) sum_{n > K} rho^{-n} 2 (h/2)^{n+1}/(n+1) + <= M_F(rho) rho^{-(K+1)} . 2 (h/2)^{K+2}/(K+2) . 1/(1 - h/(2 rho)) + +per cell (`1/(n+1) <= 1/(K+2)` and the geometric series in `h/(2 rho) = 1/24`; at most `2.2e-48` per cell here, `3.3e-47` summed over the 120 cells, about half the `D_3` width). The v1.0 text +applied Cauchy's estimate at an off-centre point `xi` with the full radius `rho`, which is not justified (the disc around +`xi` inside the bounded disc has radius `rho - |xi - s_0|`); the centre-based tail bound above replaces it (Codex +4148001628) and is what the script computes. +Beyond `T` each branch is bounded by `|const| . tail_poly(s) e^{-gamma s^2}` with `erfc <= 2`, `|J_k| <= J_k(inf) = (1/2) 8^{(k+1)/2} Gamma((k+1)/2)`, +and `int_T^inf s^j e^{-gamma s^2} ds = (1/2) gamma^{-(j+1)/2} Gamma((j+1)/2, gamma T^2)`, the upper incomplete gamma +bounded by `x^{a-1} e^{-x}/(1 - (a-1)/x)` (`a > 1`) or `x^{a-1} e^{-x}` (`a <= 1`); at `T = 60` this is below `1e-45`. +The `m = 1` and `m = 2` integrals go through the same routine (their integrands are `poly . exp` and +`poly . exp . erfc`); rule TRUNCATION_NESTING re-runs `m = 3` with `K = 20`, `h = 1` and requires the coarser enclosure +to intersect and to be wider. + +## 4. The coefficient and the torus transfer in `d = 4` + +[LP] (15.2) for the reference kernel: `c_{d,ref} = Gamma(7/6) (24)^{-1/3} pi^{-1/2} |S^{d-1}| p_G(0) p_V(0) tau^{4/3} D_{d-1}` +with `p_G(0) p_V(0) = (2 pi)^{-d}/sqrt3`, `tau^{4/3} = 6^{2/3}`, `|S^{d-1}| = 2 pi^{d/2}/Gamma(d/2)`, i.e. + + c_{d,ref} = Gamma(7/6) (3/2)^(1/3) |S^{d-1}| D_{d-1} / (sqrt3 sqrt pi (2 pi)^d): + d = 2: /(2 sqrt3 pi^(3/2)), d = 3: /(2 sqrt3 pi^(5/2)) (SIDE24 (1)), d = 4: c_{4,ref} = Gamma(7/6) (3/2)^(1/3) D_3 / (8 sqrt3 pi^(5/2)). + +The script evaluates the general formula (rule CLOSED_FORM_D4 compares it with the displayed `d = 4` form) and at +`d = 2, 3` with the certified `D_1`, `D_2` lands inside SIDE24's `c_{2,24}`, `c_{3,24}` intervals widened by their +`1e-106` periodization allowance (rule SIDE24_CONSISTENT). + +**Transfer to `K_24`, `d = 4`.** SIDE24 sections 2-4 are dimension-generic except for three constants, restated here. +(i) Image bound: for `|x| >= 1` every unit-direction contraction of order `q <= 6` of `phi = exp(-|x|^2/2)` is bounded by +`76 |x|^6 phi(x)`, and by `15` at `0` (SIDE24 section 2). In `d = 4` the shell `|n|_inf = j` holds `(2j+1)^4 - (2j-1)^4 = 64 j^3 + 16 j <= 80 j^3` +points with `|n|^2 <= 4 j^2`, so `sum_{n != 0} |n|^6 e^{-288 |n|^2} <= 5120 sum_j j^9 e^{-288 j^2} <= 10240 e^{-288}` (successive +terms have ratio below `1/2`), and every covariance entry of the 3-jet of `K_24` is within +`E_4 = 10240 (76 . 24^6 + 15) e^{-288} = 1.2462e-111` of its reference value (normalization `S >= 1`). +(ii) Sandwich: the joint covariance of `(G, t, svec H)` has dimension `n = d + 1 + d(d+1)/2 = 15`, so the spectral norm of +the difference is at most `2 n E_4 = 30 E_4`; the reference matrix satisfies `C_ref >= I/3` in every `d` (Hessian block +eigenvalues `d + 2` and `2`; odd block eigenvalues `1` and `8 +- sqrt58 > 1/3`), hence `(1 - eps) C_ref <= C_24 <= (1 + eps) C_ref` +with `eps = 90 E_4 = 1.1215e-109`, in every frame, and the same for marginals and Schur complements. +(iii) Ratios: as in SIDE24 section 4, the integrand of (15.2) changes by a factor between `(1-eps)^a/(1+eps)^b` and +`(1+eps)^a/(1-eps)^b` with `a = m + 2/3 + n_A/2 = 20/3` (`h` of degree `2m` under covariance scaling, `tau^{4/3}`, the +density-comparison constant with `n_A = m(m+1)/2 = 6`) and `b = d + n_A/2 = 7` (the two `d`-dimensional densities at `0` +and the determinant factor). Since `(1+eps)^{20/3} <= (1+eps)^7`, the ratio lies in `[1/(1+delta), 1 + delta]` with +`delta = (1+eps)^7/(1-eps)^7 - 1 = [sum_k (C(7,k) - (-1)^k C(7,k)) eps^k]/(1-eps)^7 = 14 eps + O(eps^2)`, computed without +cancellation (rule TRANSFER_BOUND certifies `13 eps <= delta <= 15 eps` for every represented value, `delta.lo >= 13 eps.hi` +and `delta.hi <= 15 eps.lo`, and records that the enclosures of `delta` and `14 eps` intersect: the offset +`delta - 14 eps = 98 eps^2 + O(eps^3) ~ 1.2e-216` lies far below the arithmetic resolution of the two enclosures, so `14 eps <= delta` +is not separable and is not claimed as certified); positivity of every reference integrand lets the +bound pass through the angular integral. So `|c_{4,24}/c_{4,ref} - 1| <= 1.5702e-108`, far below the arithmetic width of +`c_{4,ref}`, and the printed digits of `c_{4,24}` are those of `c_{4,ref}`. The exact periodic constant is not claimed equal to +the reference constant; it is enclosed. + +## 5. Values (`RESULTS.json`) + +| quantity | enclosure (leading digits) | width | +|---|---|---| +| `D_1` | `4/3` (rule D1_EXACT) | `< 1e-40` | +| `D_2` | `2.38384359055015523513604925862744194136738585267666...` = `29/6 - sqrt6` (rule D2_EXACT) | `1e-50` | +| `D_3` | `5.3231802688898968949231987396872606357843885294...` | `6.5e-47` | +| `c_{2,ref}` | `0.07340691930603427103013596295777...` (inside SIDE24's `c_{2,24}` interval) | `1.3e-54` | +| `c_{3,ref}` | `0.04177593184059834334293666542857...` (inside SIDE24's `c_{3,24}` interval) | `1.9e-52` | +| `c_{4,ref}` | `0.023321666002952835094521194952856885153101372363...` | `2.9e-49` | +| `c_{4,24}` | the same digits; `|c_{4,24}/c_{4,ref} - 1| <= 1.5702e-108` | | +| `E_4`, `eps`, `delta` | `1.2462e-111`, `1.1215e-109`, `1.5702e-108` | | +| `Z_3` | `48 sqrt2 pi = 213.25838103...` (floating cube quadrature `212.2`) | | + +Floating controls (not part of the certificate; printed to ten digits): orthant midpoint quadrature of (1.1) at +`m = 1, 2, 3` (`1.333333333`, `2.383737`, `5.3098` on a `90^3` grid) and the `(p, q)` midpoint quadrature of the +`a`-integrated integrand (`5.32335`, `1200^2` grid, error `O(h^2)`), which does not use the `s`-substitution, the +`q`-integration or the `J_k`. Offline, the same `(p, q)` quadrature at `2400^2` and `4800^2` points (`5.3232355`, +`5.3231941`) Richardson-extrapolates to `5.3231802687`, within `2e-10` of the certified `D_3`. + +## 6. Rules, mutants, verification + +`python3 -B -S cone_moment_d4.py` (also `-B -O -S`; about one minute) prints `RESULTS.json` and exits `0` only if all +thirteen rules hold: FLOAT_INSIDE (the four floating controls and the `Z_3` cube within their tolerances), D1_EXACT, +D2_EXACT, WIDTHS (`D_1`, `D_2` below `1e-40`, `D_3` and `c_{4,ref}` below `1e-30`), SIDE24_CONSISTENT, CLOSED_FORM_D4, +TRUNCATION_NESTING, MONOTONE (`D_1 < D_2 < D_3`, `c_{2,ref} > c_{3,ref} > c_{4,ref} > 0`), IMAGE_BOUND (`E_4` in +`(1.2e-111, 1.3e-111)`), TRANSFER_BOUND, MEHTA_Z3 (section 1), LIBRARY_EXACT (the decimal module's `exp` and `sqrt`, and the interval negation, +against exact rational brackets at six arguments), PINNED (forty-six digits of `D_3`, forty-eight of `c_{4,ref}`, fifty of +`pi`). Mutants (`--mutant`, each must exit `1`): `shift-variance` (`S^2/20` for `S^2/24` at `m = 3`), `vandermonde` +(the Vandermonde factor dropped), `erfc-branch` (the `erfc` branch halved), `remainder-dropped` (order `6`, cells of +width `1`, no Cauchy remainder), `mehta` (`Z_4` for `Z_3`), `sphere-area` (`4 pi` for `2 pi^2`), `image-shells` (half the +shell count), `sandwich-exponent` (`b = 5`), `recurrence` (`P_3` scaled by `101/100`). The workflow +`side24-d4-coefficient.yml` replays the manifest, the three main-resident pins, checks that the SIDE24 intervals quoted in +the script are those of `ENCLOSURE.json`, runs both interpreter modes byte for byte against `RESULTS.json` and the nine +mutants in each, on a clean tree. + +## 7. What this does not do + +Not a proof or review of [LP] or of SIDE24's sections 1-4, which are consumed at their stated scope (the reference law, +the identification of (15.2) with the persistence coefficient, the covariance sandwich). Not `d >= 5` (the `m = 4` sector +integral has three remaining variables after the `a`-integration; whether a linear substitution decouples them as +`s = 2p + q` does here is not examined). Not a window coefficient, a finite-radius quantity, `C`, `r_*` or `z_*`; not a +closed form for `D_3`: none is claimed here. Math-#201 proposes `(50 pi + 200 arctan 2 - 228)/(9 pi)`, which lies inside +this enclosure; that identity is Math-#201's claim and is neither proved nor certified by this record. No register, catalog, GRAPH or STATUS change; C8 stays OPEN. +Claude reads of this packet count for nothing; the author will not merge. + +## 8. Provenance + +Pins on `main 3e0a91b` ([LP] `dfed3b8d`, SIDE24 `PROOF.md` `44b66f04`, `ENCLOSURE.json` `57af39a0`) verified by the +workflow; read, not pinned: the SIDE24 remainder record (`0bcf0a6f`) and the C8 catalog entry (`02455f05`). +`SOURCE_FILES.json` fingerprints the packet. The general `Z_m` is Mehta's integral (the Gaussian specialization of +Selberg's integral, Mehta, *Random Matrices*, chapter 17); at `m <= 3`, the only cases used, it is verified inside the +record (`m = 1` trivially, `m = 2` by the elementary two-variable integral, `m = 3` by rule MEHTA_Z3). Cauchy's estimate +and the Lagrange remainder are elementary. The interval class is the repaired class of Math-#190 v1.3 / Math-#197 (exact +negation, explicit contexts, rule LIBRARY_EXACT). No external numerical library is used. + +## 9. Revisions + +- **v1.1 (Codex review thread 4148001628 on `54dd458`, taken).** The v1.0 text applied Cauchy's estimate off-centre without + shrinking the radius; the remainder is now the centre-based Taylor tail with the geometric factor `1/(1 - h/(2 rho))` + (section 3, `enclose_integral`). Values unchanged. +- **v1.2 (Codex thread 4148778791 on Math-#200, the same construction there; applied here for consistency).** Rule + TRANSFER_BOUND compared endpoints in the non-proving direction and section 4 described it as checking + `14 eps <= delta <= 15 eps`; it now certifies `13 eps <= delta <= 15 eps` for every represented value (`delta.lo >= 13 eps.hi`, + `delta.hi <= 15 eps.lo`) and records the intersection of the `delta` and `14 eps` enclosures (section 4). No certified value + changed. Math-#201 gives `D_3` in closed form, `(50 pi + 200 arctan 2 - 228)/(9 pi)`, at 47 common digits with this record's + enclosure; this record is otherwise unchanged and section 7 now cites it without claiming it. +- **v1.3 (nonauthor read S1, Grok Bot agent 1, Math-#199 review 5403487705, on `617236d`; all findings P3).** No certified + value or `RESULTS.json` byte changed. F1: the displayed `D_3` digits end at the certified prefix `…3885294…` (the old trailing + `09` lay past it; the true digits are `…38852944092…`), and the width is `6.5e-47`. F2: the widths of `c_{2,ref}` (`1.3e-54`) + and `c_{4,ref}` (`2.9e-49`) and the runtime (about one minute) are updated. F3: the per-cell remainder is at most `2.2e-48`, + not of order `1e-33`. F4: `delta - 14 eps = 98 eps^2 + O(eps^3)`, in the NOTE and in the script comment. F5: section 7 cites + Math-#201's closed form without claiming it. F6: the workflow accepts a mutant only if it exits 1 with JSON `passed: false` + (rejection by rules, not a crash). F7: the script checks `gamma T^2 > a - 1` before the `Gamma(a, x)` tail bound for `a > 1` + (it holds here, `x >= 150`). diff --git a/coefficients/side24_d4_20260930/RESULTS.json b/coefficients/side24_d4_20260930/RESULTS.json new file mode 100644 index 000000000..187552f59 --- /dev/null +++ b/coefficients/side24_d4_20260930/RESULTS.json @@ -0,0 +1,143 @@ +{ + "object": "CL-SIDE24-D4-COEFFICIENT-20260930-v1", + "scientific_effect": "NONE", + "certified": true, + "mutant": null, + "passed": true, + "checks": { + "FLOAT_INSIDE": true, + "D1_EXACT": true, + "D2_EXACT": true, + "WIDTHS": true, + "SIDE24_CONSISTENT": true, + "CLOSED_FORM_D4": true, + "TRUNCATION_NESTING": true, + "MONOTONE": true, + "IMAGE_BOUND": true, + "TRANSFER_BOUND": true, + "MEHTA_Z3": true, + "LIBRARY_EXACT": true, + "PINNED": true + }, + "cone_moments": { + "D_1": [ + "1.33333333333333333333333333333333333333333333333333332172755", + "1.33333333333333333333333333333333333333333333333333334493939" + ], + "D_2": [ + "2.38384359055015523513604925862744194136738585267665775580671", + "2.38384359055015523513604925862744194136738585267666865399451" + ], + "D_2_exact_29/6-sqrt6": [ + "2.38384359055015523513604925862744194136738585267666320490062", + "2.38384359055015523513604925862744194136738585267666320490067" + ], + "D_3": [ + "5.32318026888989689492319873968726063578438852940845576405060", + "5.32318026888989689492319873968726063578438852947338433618854" + ], + "D_3_coarse": [ + "5.32318026888935064250203575239838607829572344863496106068363", + "5.32318026889044314734430434611621544000761990335365792727914" + ] + }, + "coefficients": { + "c_2_ref": [ + "0.0734069193060342710301359629577740500176642446843326711937355", + "0.0734069193060342710301359629577740500176642446843326724906200" + ], + "c_3_ref": [ + "0.0417759318405983433429366654285755564666815196610446423810796", + "0.0417759318405983433429366654285755564666815196610448333783896" + ], + "c_4_ref": [ + "0.0233216660029528350945211949528568851531013723630981034059901", + "0.0233216660029528350945211949528568851531013723633825654059788" + ], + "c_4_ref_closed": [ + "0.0233216660029528350945211949528568851531013723630981034060010", + "0.0233216660029528350945211949528568851531013723633825654059682" + ], + "c_4_24": [ + "0.0233216660029528350945211949528568851531013723630981034059898", + "0.0233216660029528350945211949528568851531013723633825654059791" + ] + }, + "transfer": { + "E_4": [ + "1.24614968346183326033273338980201191527161500226085180392811E-111", + "1.24614968346183326033273338980201191527161500226085180392812E-111" + ], + "epsilon": [ + "1.12153471511564993429946005082181072374445350203476662353529E-109", + "1.12153471511564993429946005082181072374445350203476662353531E-109" + ], + "delta": [ + "1.5701486011619099080192440711505350132422349028486732729494E-108", + "1.57014860116190990801924407115053501324223490284867327294949E-108" + ], + "n_jet": 15, + "exponents": { + "a": "20/3 (<= 7 used)", + "b": 7 + } + }, + "constants": { + "pi": [ + "3.14159265358979323846264338327950288419716939937510582097457", + "3.14159265358979323846264338327950288419716939937510582097528" + ], + "Gamma(7/6)": [ + "0.927719333630039200708349482534621018566466519145475576700351", + "0.927719333630039200708349482534621018566466519145475576939062" + ], + "(3/2)^(1/3)": [ + "1.14471424255333186780804221193967700891590692078793107209905", + "1.14471424255333186780804221193967700891590692078793107209906" + ], + "Z_1": [ + "3.54490770181103205459633496668229036559509891224477425642683", + "3.54490770181103205459633496668229036559509891224477425642838" + ], + "Z_2": [ + "20.0530261970480040193261222784883620240558939248795065330326", + "20.0530261970480040193261222784883620240558939248795065330462" + ], + "Z_3": [ + "213.258381031601579856762287522913297533501841090033130707983", + "213.258381031601579856762287522913297533501841090033130708211" + ], + "Z_3_quadrature": [ + "213.258381031601579856762287522913297533501841090033130694877", + "213.258381031601579856762287522913297533501841090033130721327" + ], + "|S^3|": [ + "19.7392088021787172376689819997523022706273988144815812528219", + "19.7392088021787172376689819997523022706273988144815812528310" + ] + }, + "quadrature": { + "K": 40, + "cell": "1/2", + "T": "60", + "rho": "6" + }, + "float_reference_10sig": { + "D_1_orthant": "1.333333333", + "D_2_orthant": "2.383737252", + "D_3_orthant": "5.309848377", + "D_3_pq_quadrature": "5.32335234", + "Z_3_cube": "212.2117026" + }, + "exact_relations": { + "D_m": "sqrt(3/(m+3)) Z_m^-1 int_{mu<0} prod mu_i^2 |Delta(mu)| exp(-|mu|^2/4 + (sum mu)^2/(4(m+3))) dmu", + "Z_m": "2^(m/2 + m(m-1)/4) (2 pi)^(m/2) prod_{j<=m} Gamma(1 + j/2) / Gamma(3/2)^m; Z_3 = (3/2) sqrt(pi/3) int_0^inf e^{-s^2/24} (4 s^2 - 32 + 32 e^{-s^2/8}) ds = 48 sqrt2 pi", + "m=3 sector": "mu_3 = -a, mu_2 = -a-p, mu_1 = -a-p-q; Q = 3a^2/8 + p^2/3 + 5q^2/24 + ap/2 + aq/4 + pq/3", + "a-integral": "G_0 = (1/2) sqrt(8 pi/3) exp(p^2/6 + pq/6 + q^2/24) erfc((2p+q)/(2 sqrt6)); G_1 = (1 - beta G_0)/(2 alpha); G_{n+1} = (n G_{n-1} - beta G_n)/(2 alpha)", + "decoupling": "s = 2p + q: erfc branch exp(-s^2/24 - q^2/8), polynomial branch exp(-s^2/12 - q^2/8); dp dq = ds dq/2, 0 < q < s", + "q-integral": "J_0 = sqrt(2 pi) erf(s/(2 sqrt2)), J_1 = 4(1 - e^{-s^2/8}), J_{k+2} = -4 s^{k+1} e^{-s^2/8} + 4(k+1) J_k", + "D_3": "I_3 / (16 pi)", + "c_d_ref": "Gamma(7/6) (3/2)^(1/3) |S^{d-1}| D_{d-1} / (sqrt3 sqrt pi (2 pi)^d); d = 4: Gamma(7/6) (3/2)^(1/3) D_3 / (8 sqrt3 pi^(5/2))", + "transfer": "E_4 = 10240 (76 24^6 + 15) e^-288; eps = 90 E_4; c_{4,24}/c_{4,ref} in [(1-eps)^7/(1+eps)^7, (1+eps)^7/(1-eps)^7]" + } +} diff --git a/coefficients/side24_d4_20260930/SOURCE_FILES.json b/coefficients/side24_d4_20260930/SOURCE_FILES.json new file mode 100644 index 000000000..06e5cb11e --- /dev/null +++ b/coefficients/side24_d4_20260930/SOURCE_FILES.json @@ -0,0 +1,67 @@ +{ + "schema": 1, + "object": "CL-SIDE24-D4-COEFFICIENT-20260930-v1", + "scientific_effect": "NONE", + "certified": true, + "review_status": "AUTHOR_SIDE_CERTIFIED_CONSTANT_NONAUTHOR_READS_OPEN", + "disposition": "CERTIFIED_ENCLOSURE_OF_THE_STATED_EXPRESSION_IN_D4_NO_REGISTER_CHANGE_SIDE24_V1_UNCHANGED", + "declarative": true, + "executed": false, + "base": "3e0a91b69892b9717564c592e7a966492a81566a", + "claim": "https://github.com/d6g8k5htny-coder/Math-/pull/197#issuecomment-5916664092", + "pins_on_main": [ + { + "label": "LP", + "path": "imports/lifetime_parent_20260925/UNIFORM_MATRIX_CAP_AND_LIFETIME.md", + "git_blob": "dfed3b8d318a3ab1950957f393307733a4bef3f2", + "role": "(15.1) the cone moment D_u; (15.2) the coefficient c_{d,L} with the angular integral, p_G(0) p_Vu(0), tau_u^(4/3)", + "workflow_checked": true + }, + { + "label": "S24", + "path": "coefficients/side24_v1/PROOF.md", + "git_blob": "44b66f04f89fcd87383b3603fa69f1feb64cdddd", + "role": "section 1 reference law A = Q + sqrt(2/3) Z I_m, p_G p_V = (2 pi)^-d/sqrt3, tau^2 = 6, D_1 = 4/3, D_2 = 29/6 - sqrt6; sections 2-4 image bound, covariance sandwich and cone-moment comparison generalized to d = 4 in NOTE section 4", + "workflow_checked": true + }, + { + "label": "ENC", + "path": "coefficients/side24_v1/ENCLOSURE.json", + "git_blob": "57af39a05e14ed0ba8ebc00a9b4aca4dffb067c7", + "role": "the d = 2, 3 enclosures quoted in the script for rule SIDE24_CONSISTENT", + "workflow_checked": true + } + ], + "read_not_pinned": [ + { + "path": "coefficients/side24_remainder_20260929/PROOF.md", + "git_blob_at_base": "0bcf0a6f5a608ce369e769103cd07fe71e0da844", + "role": "Proposition V (first nearest-image variation for every d): not used; the transfer here is the section 3-4 sandwich" + }, + { + "path": "reviews/candidates_pending_20260928/CANDIDATES.md", + "git_blob_at_base": "02455f0561d51990e1486e5d4fbb4229c8a34f5a", + "role": "catalog entry C8 listing c_{d,L} among the unsupplied constants; C8 stays OPEN" + } + ], + "files": [ + { + "path": "NOTE.md", + "bytes": 20408, + "sha256": "338ddfe869542d08e3595a34d11ca1c8cabc9248e167663e09a178ccb3066abc", + "git_blob": "e00de4665594bad39545046433989bfa5178cd65" + }, + { + "path": "RESULTS.json", + "bytes": 5325, + "sha256": "804200e99d45e4b019a0e46709e912a14f0964730a461b448e4beb1d70ce2710", + "git_blob": "187552f596576c9dcdecd1deabc0304412d8a1f7" + }, + { + "path": "cone_moment_d4.py", + "bytes": 37638, + "sha256": "6b13502fe00f331af679c02baa2816e1b668fd82c0c37c7e6102b2f3a156e7c3", + "git_blob": "3d9f42a833871e45c4a34152062238feb241154c" + } + ] +} diff --git a/coefficients/side24_d4_20260930/cone_moment_d4.py b/coefficients/side24_d4_20260930/cone_moment_d4.py new file mode 100644 index 000000000..3d9f42a83 --- /dev/null +++ b/coefficients/side24_d4_20260930/cone_moment_d4.py @@ -0,0 +1,874 @@ +"""SIDE24's coefficient in d = 4: the birth-integrated cone moment D_3, certified, with c_{4,ref} and c_{4,24}. + +Standard library only. Run from anywhere: python -B -S cone_moment_d4.py [--mutant NAME] +Every printed interval [lo, hi] is a rigorous enclosure: interval arithmetic over `decimal` (60 digits, directed rounding; +exp and sqrt of the module widened by two units in the last place and checked against exact rational brackets by rule +LIBRARY_EXACT), pi by Machin, fractional powers by exact integer roots, erf by its positive series with a geometric tail, +and the one-dimensional integrals by an order-K Taylor expansion at each cell centre whose remainder is bounded by +Cauchy's estimate on a disc of radius RHO around the centre (every factor of the integrand is entire; its modulus on the +disc is bounded elementarily). Nothing is fitted, sampled or extrapolated. + +Object. For the Euclidean reference kernel of coefficients/side24_v1 section 1 in dimension d = m + 1, the transverse +Hessian given V = 0 is A = Q + sqrt(2/3) Z I_m with Q the m x m GOE of density proportional to exp(-tr Q^2/4) and Z an +independent standard normal, and the cone moment of [LP] (15.1) is D_m = E[det(A)^2 1{A < 0}]. Integrating Z against the +GOE eigenvalue density (NOTE section 1) gives, with S = sum mu_i, + + D_m = sqrt(3/(m+3)) / Z_m . int_{mu in (-inf,0)^m} prod mu_i^2 prod_{i 0) the exponent is +Q = (3/8) a^2 + (1/3) p^2 + (5/24) q^2 + (1/2) a p + (1/4) a q + (1/3) p q, the a-integral of a^n exp(-alpha a^2 - beta a), +alpha = 3/8, beta = p/2 + q/4, is exact (G_0 = (1/2) sqrt(pi/alpha) e^{beta^2/(4 alpha)} erfc(beta/(2 sqrt alpha)), +G_1 = (1 - beta G_0)/(2 alpha), G_{n+1} = (n G_{n-1} - beta G_n)/(2 alpha)), and in s = 2p + q the remaining quadratic forms +decouple: s^2/24 + q^2/8 on the erfc branch, s^2/12 + q^2/8 on the polynomial branch, so the q-integral over (0, s) is exact +through J_k(s) = int_0^s q^k e^{-q^2/8} dq. Then D_3 = I_3/(16 pi) with I_3 a one-dimensional integral (NOTE section 2). + +Coefficient. [LP] (15.2) with |S^{d-1}| = 2 pi^{d/2}/Gamma(d/2), p_G(0) p_V(0) = (2 pi)^{-d}/sqrt3 and tau^2 = 6: +c_{d,ref} = Gamma(7/6) (3/2)^{1/3} |S^{d-1}| D_{d-1} / (sqrt3 sqrt pi (2 pi)^d); d = 2, 3 reproduce SIDE24 (1) and its intervals +(a rule); d = 4: c_{4,ref} = Gamma(7/6) (3/2)^{1/3} D_3 / (8 sqrt3 pi^{5/2}). The torus value c_{4,24} follows from the SIDE24 +section 3-4 covariance sandwich with the d = 4 constants (NOTE section 4). Scientific effect NONE. +""" +import argparse +import itertools +import json +import math +import sys +import time +from decimal import Decimal, Context, ROUND_FLOOR, ROUND_CEILING, ROUND_HALF_EVEN, setcontext +from fractions import Fraction as Fr + +MUTANTS = ("shift-variance", "vandermonde", "erfc-branch", "remainder-dropped", "mehta", "sphere-area", "image-shells", + "sandwich-exponent", "recurrence") +MUT = None + +PREC = 60 +CF = Context(prec=PREC, rounding=ROUND_FLOOR, Emin=-9999999, Emax=9999999) +CC = Context(prec=PREC, rounding=ROUND_CEILING, Emin=-9999999, Emax=9999999) +CE = Context(prec=PREC, rounding=ROUND_HALF_EVEN, Emin=-9999999, Emax=9999999) +# Every arithmetic step names one of the three contexts; operator syntax on Decimals would round to the current context +# (default 28 digits), so it is widened as a safety net and rule LIBRARY_EXACT checks exp, sqrt and the negation. +setcontext(Context(prec=2 * PREC + 20, rounding=ROUND_HALF_EVEN, Emin=-9999999, Emax=9999999)) +ZERO, ONE, HALF = Decimal(0), Decimal(1), Decimal("0.5") +ROOT_DIGITS = 70 +K_ORDER, RHO, CELL, T_CUT = 40, Fr(6), Fr(1, 2), Fr(60) +SIDE24 = {2: ("0.07340691930603427103", "0.07340691930603427104"), 3: ("0.04177593184059834334", "0.04177593184059834335")} +PINNED = {"D_3": "5.323180268889896894923198739687260635784388529", + "c_4_ref": "0.023321666002952835094521194952856885153101372363", + "pi": "3.14159265358979323846264338327950288419716939937510"} + + +def ulp(x): + if x == 0: + return Decimal(1).scaleb(-PREC - 30) + return Decimal(1).scaleb(x.adjusted() - PREC + 1) + + +class Iv: + """Closed interval [lo, hi] of Decimals; every operation encloses the exact result.""" + __slots__ = ("lo", "hi") + + def __init__(self, lo, hi=None): + if hi is None: + hi = lo + if lo > hi: + raise ValueError("empty interval") + self.lo, self.hi = lo, hi + + @staticmethod + def frac(q): + q = Fr(q) + n, d = Decimal(q.numerator), Decimal(q.denominator) + return Iv(CF.divide(n, d), CC.divide(n, d)) + + @staticmethod + def of(x): + return x if isinstance(x, Iv) else Iv.frac(x) + + def __add__(a, b): + b = Iv.of(b) + return Iv(CF.add(a.lo, b.lo), CC.add(a.hi, b.hi)) + __radd__ = __add__ + + def __neg__(a): + return Iv(a.hi.copy_negate(), a.lo.copy_negate()) + + def __sub__(a, b): + b = Iv.of(b) + return Iv(CF.subtract(a.lo, b.hi), CC.subtract(a.hi, b.lo)) + + def __rsub__(a, b): + return Iv.of(b) - a + + def __mul__(a, b): + b = Iv.of(b) + c = [(x, y) for x in (a.lo, a.hi) for y in (b.lo, b.hi)] + return Iv(min(CF.multiply(x, y) for x, y in c), max(CC.multiply(x, y) for x, y in c)) + __rmul__ = __mul__ + + def __truediv__(a, b): + b = Iv.of(b) + if b.lo <= 0 <= b.hi: + raise ZeroDivisionError("interval division by an interval containing 0") + c = [(x, y) for x in (a.lo, a.hi) for y in (b.lo, b.hi)] + return Iv(min(CF.divide(x, y) for x, y in c), max(CC.divide(x, y) for x, y in c)) + + def __rtruediv__(a, b): + return Iv.of(b) / a + + def sq(a): + if a.lo <= 0 <= a.hi: + return Iv(ZERO, max(CC.multiply(a.lo, a.lo), CC.multiply(a.hi, a.hi))) + return a * a + + def exp(a): + lo, hi = CE.exp(a.lo), CE.exp(a.hi) + return Iv(max(ZERO, CF.subtract(lo, CC.multiply(Decimal(2), ulp(lo)))), CC.add(hi, CC.multiply(Decimal(2), ulp(hi)))) + + def sqrt(a): + if a.lo < 0: + raise ValueError("sqrt of an interval with negative part") + lo, hi = CE.sqrt(a.lo), CE.sqrt(a.hi) + return Iv(max(ZERO, CF.subtract(lo, CC.multiply(Decimal(2), ulp(lo)))), CC.add(hi, CC.multiply(Decimal(2), ulp(hi)))) + + def width(a): + return CC.subtract(a.hi, a.lo) + + def mag(a): + return max(abs(a.lo), abs(a.hi)) + + def contains(a, x): + x = Iv.of(x) + return a.lo <= x.lo and x.hi <= a.hi + + def intersects(a, b): + return a.lo <= b.hi and b.lo <= a.hi + + def pair(a): + return [str(a.lo), str(a.hi)] + + +def hull(a, b): + return Iv(min(a.lo, b.lo), max(a.hi, b.hi)) + + +def ipow(x, n): + r = Iv(ONE) + for _ in range(n): + r = r * x + return r + + +# ---------------- constants and special functions +def inroot(N, n): + if N < 2: + return N + x = 1 << ((N.bit_length() + n - 1) // n) + while True: + y = ((n - 1) * x + N // pow(x, n - 1)) // n + if y >= x: + break + x = y + if not (pow(x, n) <= N < pow(x + 1, n)): + raise ValueError("integer root failed") + return x + + +def root_frac(q, n): + q = Fr(q) + if q <= 0: + raise ValueError("root of a nonpositive number") + p, d = q.numerator, q.denominator + scale = 10 ** ROOT_DIGITS + r = inroot(p * d ** (n - 1) * scale ** n, n) + den = d * scale + return Iv(Iv.frac(Fr(r, den)).lo, Iv.frac(Fr(r + 1, den)).hi) + + +def pi_interval(): + def atan_inv(n, terms): + x = Fr(1, n) + partial, sums = Iv(ZERO), [] + for m in range(terms + 2): + t = Iv.frac(x ** (2 * m + 1) / (2 * m + 1)) + partial = partial + t if m % 2 == 0 else partial - t + sums.append(partial) + return hull(sums[-1], sums[-2]) + return atan_inv(5, 44) * 16 - atan_inv(239, 22) * 4 + + +PI = pi_interval() +SQRTPI = PI.sqrt() +SQRT2 = Iv.frac(2).sqrt() +SQRT3 = Iv.frac(3).sqrt() +SQRT6 = Iv.frac(6).sqrt() +SQRT2PI = (2 * PI).sqrt() +TINY = Decimal(10) ** (-PREC - 8) + + +def erf_point(x): + """erf(x) at a Decimal point by the positive series (2/sqrt pi) e^{-x^2} sum 2^n x^{2n+1}/(2n+1)!!; once the term ratio + 2x^2/(2n+3) is at most 1/2 the omitted tail is at most twice the first omitted term.""" + if x < 0: + return -erf_point(-x) + if x == 0: + return Iv(ZERO) + X = Iv(x) + x2 = X * X + t, S, n = X, Iv(ZERO), 0 + while True: + S = S + t + ratio = x2 * 2 / (2 * n + 3) + t = t * ratio + n += 1 + if ratio.hi <= HALF and t.hi < S.lo * TINY: + break + if n > 200000: + raise RuntimeError("erf series did not converge") + return (S + Iv(ZERO, CC.multiply(Decimal(2), t.hi))) * (-x2).exp() * 2 / SQRTPI + + +def erf_iv(x): + return Iv(erf_point(x.lo).lo, erf_point(x.hi).hi) + + +def gamma_lower_point(x, a): + """gamma(a, x) by the positive series x^a e^-x sum_n x^n/(a (a+1) ... (a+n)); a = 7/6, x a Decimal > 0.""" + X = Iv(x) + t, S, n = Iv.frac(1 / a), Iv(ZERO), 0 + while True: + S = S + t + n += 1 + ratio = X / Iv.frac(a + n) + t = t * ratio + if ratio.hi <= HALF and t.hi < S.lo * TINY: + break + xa = X * root_frac(Fr(x), 6) + return xa * (-X).exp() * (S + Iv(ZERO, CC.multiply(Decimal(2), t.hi))) + + +def gamma_complete(a, X=Decimal(120)): + """Gamma(a) = gamma(a, X) + Gamma(a, X), 1 < a < 2, X^(a-1) e^-X <= Gamma(a, X) <= X^(a-1) e^-X / (1 - (a-1)/X).""" + low = gamma_lower_point(X, a) + up = root_frac(Fr(X), 6) * (-Iv(X)).exp() + return low + Iv(up.lo, (up / (1 - Iv.frac(a - 1) / Iv(X))).hi) + + +GAMMA76 = gamma_complete(Fr(7, 6)) +CBRT32 = root_frac(Fr(12), 3) / 2 # (3/2)^(1/3) + + +def gamma_half_integer(twice): + """Gamma(twice/2) for a positive integer `twice`, exactly in terms of sqrt pi.""" + if twice % 2 == 0: + return Iv.frac(math.factorial(twice // 2 - 1)) + k = (twice - 1) // 2 # Gamma(k + 1/2) = (2k)! sqrt pi / (4^k k!) + return Iv.frac(Fr(math.factorial(2 * k), 4 ** k * math.factorial(k))) * SQRTPI + + +def mehta_Z(m): + """Z_m = int_{R^m} prod_{i Fraction) +def padd(u, v): + w = dict(u) + for k, c in v.items(): + w[k] = w.get(k, Fr(0)) + c + return {k: c for k, c in w.items() if c != 0} + + +def pmul(u, v): + w = {} + for k1, c1 in u.items(): + for k2, c2 in v.items(): + k = tuple(a + b for a, b in zip(k1, k2)) + w[k] = w.get(k, Fr(0)) + c1 * c2 + return {k: c for k, c in w.items() if c != 0} + + +def pscale(u, c): + return {k: v * c for k, v in u.items() if v * c != 0} + + +def pconst(c, nvars): + return {(0,) * nvars: Fr(c)} if c != 0 else {} + + +def pvar(i, nvars): + e = [0] * nvars + e[i] = 1 + return {tuple(e): Fr(1)} + + +def ppow(u, n, nvars): + r = pconst(1, nvars) + for _ in range(n): + r = pmul(r, u) + return r + + +def psubst(u, images, nvars_out): + """Substitute variable i -> images[i] (polynomials in nvars_out variables).""" + out = {} + for k, c in u.items(): + term = pconst(c, nvars_out) + for i, e in enumerate(k): + if e: + term = pmul(term, ppow(images[i], e, nvars_out)) + out = padd(out, term) + return out + + +def univariate(u): + """dict {(n,): c} -> list of Fractions by power.""" + if not u: + return [Fr(0)] + deg = max(k[0] for k in u) + return [u.get((n,), Fr(0)) for n in range(deg + 1)] + + +def reduce_sector(m): + """The m-dimensional sector integrand in (a, p_1, ..., p_{m-1}) with mu_m = -a, mu_{m-1} = -a - p_1, ...: + returns (polynomial, quadratic form) with polynomial = prod mu_i^2 prod_{i 0: + acc = padd(acc, pvar(i, nv)) + mus.append(pscale(acc, Fr(-1))) + poly = pconst(1, nv) + for mu in mus: + poly = pmul(poly, pmul(mu, mu)) + if MUT != "vandermonde": + for i in range(m): + for j in range(i + 1, m): + poly = pmul(poly, padd(mus[i], pscale(mus[j], Fr(-1)))) # mu_(later index) is larger: mus[i] - mus[j] > 0 + S = {} + norm2 = {} + for mu in mus: + S = padd(S, mu) + norm2 = padd(norm2, pmul(mu, mu)) + shift = Fr(1, 4 * (m + 3)) if MUT != "shift-variance" or m != 3 else Fr(1, 20) + quad = padd(pscale(norm2, Fr(1, 4)), pscale(pmul(S, S), -shift)) + return poly, quad + + +# ---------------- Taylor-model quadrature with Cauchy remainders +class Factor: + """An entire factor g(s): point Taylor coefficients at a rational centre s0 (list of K+1 intervals) and a bound M on + |g(z)| for |z - s0| <= rho.""" + + def coefficients(self, s0, K): + raise NotImplementedError + + def disc_bound(self, s0, rho): + raise NotImplementedError + + +class Gauss(Factor): + """exp(-gamma s^2).""" + + def __init__(self, gamma): + self.gamma = Fr(gamma) + + def coefficients(self, s0, K): + g = Iv.frac(self.gamma) + S0 = Iv.frac(s0) + y = [(-(g * S0.sq())).exp()] + prev = Iv(ZERO) + for n in range(K): + nxt = -(2 * g) * (S0 * y[n] + prev) / (n + 1) # (n+1) y_{n+1} = -2 gamma (s0 y_n + y_{n-1}) + prev = y[n] + y.append(nxt) + return y + + def disc_bound(self, s0, rho): + g = Iv.frac(self.gamma) + inner = max(Fr(0), abs(Fr(s0)) - Fr(rho)) + return (g * Iv.frac(rho * rho) - g * Iv.frac(inner * inner)).exp() # |e^{-g z^2}| = e^{-g(x^2 - y^2)} + + +class ErfLin(Factor): + """erf(kappa s) (sign +1) or erfc(kappa s) (sign -1); kappa an interval.""" + + def __init__(self, kappa, sign): + self.kappa, self.sign = kappa, sign + + def coefficients(self, s0, K): + k = self.kappa + S0 = Iv.frac(s0) + e = erf_iv(k * S0) + c0 = e if self.sign > 0 else 1 - e + # derivative of erf(kappa s) is (2 kappa/sqrt pi) exp(-kappa^2 s^2): coefficients of the Gaussian in s + k2 = k.sq() + y = [(-(k2 * S0.sq())).exp()] + prev = Iv(ZERO) + for n in range(K): + nxt = -(2 * k2) * (S0 * y[n] + prev) / (n + 1) + prev = y[n] + y.append(nxt) + pref = 2 * k / SQRTPI * self.sign + return [c0] + [pref * y[n] / (n + 1) for n in range(K)] + + def disc_bound(self, s0, rho): + k = self.kappa + seg = 2 * k / SQRTPI * Iv.frac(abs(Fr(s0)) + Fr(rho)) * (k.sq() * Iv.frac(Fr(rho) ** 2)).exp() + return seg if self.sign > 0 else 1 + seg + + +class Poly(Factor): + """Polynomial with rational coefficients (list by power).""" + + def __init__(self, coeffs): + self.c = list(coeffs) + + def coefficients(self, s0, K): + # Taylor coefficients at s0: sum_j c_j C(j, n) s0^(j-n) + out = [] + for n in range(K + 1): + v = Fr(0) + for j, cj in enumerate(self.c): + if j >= n: + v += cj * math.comb(j, n) * Fr(s0) ** (j - n) + out.append(Iv.frac(v)) + return out + + def disc_bound(self, s0, rho): + r = abs(Fr(s0)) + Fr(rho) + return Iv.frac(sum(abs(cj) * r ** j for j, cj in enumerate(self.c))) + + +class Jk(Factor): + """J_k(s) = int_0^s q^k e^{-q^2/8} dq: J_0 = sqrt(2 pi) erf(s/(2 sqrt2)), J_1 = 4 (1 - e^{-s^2/8}), + J_{k+2}(s) = -4 s^{k+1} e^{-s^2/8} + 4 (k+1) J_k(s).""" + + def __init__(self, k): + self.k = k + + _cache = {} + + def coefficients(self, s0, K): + key = (s0, K) + chain = Jk._cache.get(key) + if chain is None or len(chain) <= self.k: + E3 = Gauss(Fr(1, 8)).coefficients(s0, K) + J0 = [c * SQRT2PI for c in ErfLin(1 / (2 * SQRT2), +1).coefficients(s0, K)] + J1 = [Iv.frac(4) - 4 * E3[0]] + [-4 * c for c in E3[1:]] + chain = [J0, J1] + for k in range(max(self.k, 9) - 1): + sp = Poly([Fr(0)] * (k + 1) + [Fr(1)]).coefficients(s0, K) # s^{k+1} + chain.append(t_add(t_scale(t_mul(sp, E3), Fr(-4)), t_scale(chain[k], Fr(4 * (k + 1))))) + Jk._cache = {key: chain} + return chain[self.k] + + def infinity(self): + """J_k(inf) = (1/2) 8^{(k+1)/2} Gamma((k+1)/2), as an interval.""" + k = self.k + eight = Iv.frac(8 ** ((k + 1) // 2)) * (SQRT2 * 2 if (k + 1) % 2 else Iv(ONE)) + return eight * gamma_half_integer(k + 1) / 2 + + def disc_bound(self, s0, rho): + r = Iv.frac(abs(Fr(s0)) + Fr(rho)) + return ipow(r, self.k + 1) * (Iv.frac(Fr(rho) ** 2 / 8)).exp() / (self.k + 1) + + +def t_mul(u, v): + K = len(u) - 1 + out = [] + for n in range(K + 1): + acc = Iv(ZERO) + for i in range(n + 1): + acc = acc + u[i] * v[n - i] + out.append(acc) + return out + + +def t_add(u, v): + return [x + y for x, y in zip(u, v)] + + +def t_scale(u, c): + return [x * c for x in u] + + +def enclose_integral(branches, K=None, cell=None, T=None, rho=None, remainder=True): + """int_0^T of sum over branches of const * prod(prefactors) * sum_pairs poly * J, plus the tail bound. A branch is + (const Iv, [prefactor factors], [(Poly, Jk or None)], tail_poly, gamma): for s >= T the branch is bounded in absolute + value by |const| tail_poly(s) e^{-gamma s^2} (tail_poly a rational polynomial by power). On each cell the point Taylor + coefficients of the product are integrated exactly; the tail sum_{n > K} c_n x^n of the (everywhere convergent) Taylor + series at the centre is bounded through Cauchy's estimate at the centre, |c_n| <= M_F(rho)/rho^n, M_F the product/sum of + the factors' disc bounds, and the geometric series in |x|/rho <= h/(2 rho) (v1.1: the v1.0 text applied the estimate + off-centre without shrinking the radius; Codex 4148001628).""" + K = K_ORDER if K is None else K + cell = CELL if cell is None else cell + T = T_CUT if T is None else T + rho = RHO if rho is None else rho + ncell = int(T / cell) + if ncell * cell != T: + raise ValueError("cell does not divide T") + half = cell / 2 + total = Iv(ZERO) + even_moments = [Iv.frac(2 * half ** (n + 1) / (n + 1)) if n % 2 == 0 else None for n in range(K + 1)] + # int_{-h/2}^{h/2} sum_{n>K} |c_n| |x|^n dx <= M_F sum_{n>K} rho^-n 2 (h/2)^(n+1)/(n+1) + # <= M_F rho^-(K+1) 2 (h/2)^(K+2)/(K+2) . 1/(1 - h/(2 rho)) + if Fr(half) >= Fr(rho): + raise ValueError("cell half-width must be below the Cauchy radius") + rem_moment = Iv.frac(2 * half ** (K + 2) / (K + 2) / (1 - Fr(half) / Fr(rho))) + rho_pow = Iv.frac(Fr(rho) ** (K + 1)) + for c in range(ncell): + s0 = (2 * c + 1) * half + cell_val = Iv(ZERO) + M_F = Iv(ZERO) + for const, prefactors, pairs, _, _ in branches: + acc, M_acc = None, Iv(ZERO) + for poly, jk in pairs: + pc = poly.coefficients(s0, K) + M_pair = poly.disc_bound(s0, rho) + if jk is not None: + pc = t_mul(pc, jk.coefficients(s0, K)) + M_pair = M_pair * jk.disc_bound(s0, rho) + acc = pc if acc is None else t_add(acc, pc) + M_acc = M_acc + M_pair + M = M_acc + for f in prefactors: + acc = t_mul(acc, f.coefficients(s0, K)) + M = M * f.disc_bound(s0, rho) + for n in range(0, K + 1, 2): + cell_val = cell_val + const * acc[n] * even_moments[n] + M_F = M_F + Iv(const.mag()) * M + if remainder: + R = (M_F / rho_pow * rem_moment).hi + cell_val = cell_val + Iv(-R, R) + total = total + cell_val + # tail: int_T^inf s^j e^{-gamma s^2} ds = (1/2) gamma^{-(j+1)/2} Gamma((j+1)/2, gamma T^2), + # Gamma(a, x) <= x^(a-1) e^-x / (1 - (a-1)/x) for a > 1, <= x^(a-1) e^-x for a <= 1 + tail = Iv(ZERO) + for const, _, _, tail_poly, gamma in branches: + g = Iv.frac(gamma) + x = g * Iv.frac(T * T) + for j, cj in enumerate(tail_poly): + if cj == 0: + continue + a = Fr(j + 1, 2) + if a.denominator == 2: + xa1 = ipow(x, int(a - Fr(1, 2))) * x.sqrt() / x + gpow = ipow(g, int(a)) * g.sqrt() + else: + xa1 = ipow(x, int(a) - 1) + gpow = ipow(g, int(a)) + up = xa1 * (-x).exp() + if a > 1: + if not x.lo > Iv.frac(a - 1).hi: + raise ValueError("tail bound needs gamma T^2 > a - 1") + up = up / (1 - Iv.frac(a - 1) / x) + tail = tail + Iv(const.mag()) * Iv.frac(abs(cj)) * up / gpow / 2 + return total + Iv(-tail.hi, tail.hi) + + +# ---------------- the reductions +def d1_terms(): + """D_1 = (sqrt3/2) (1/(2 sqrt pi)) int_0^inf mu^2 e^{-3 mu^2/16} dmu (m = 1: no a-integration needed).""" + const = SQRT3 / 2 / (2 * SQRTPI) + return [(const, [Gauss(Fr(3, 16))], [(Poly([Fr(0), Fr(0), Fr(1)]), None)], [Fr(0), Fr(0), Fr(1)], Fr(3, 16))] + + +def a_integrated_terms(m): + """For m = 2 and m = 3: integrate the sector polynomial times e^{-Q} over a in (0, inf) exactly; returns the erfc-branch + and polynomial-branch integrands as polynomials in the remaining variables with their (decoupled) exponents.""" + poly, quad = reduce_sector(m) + nv = m + # quad = alpha a^2 + a * beta(p) + rest(p) + alpha = quad.get((2,) + (0,) * (nv - 1), Fr(0)) + beta = {k[1:]: c for k, c in quad.items() if k[0] == 1} + rest = {k[1:]: c for k, c in quad.items() if k[0] == 0} + # G_n = P_n(beta) G_0 + R_n(beta), G_0 = (1/2) sqrt(pi/alpha) e^{beta^2/(4 alpha)} erfc(beta/(2 sqrt alpha)) + nvr = nv - 1 + P = [pconst(1, nvr), pscale(beta, Fr(-1, 2) / alpha)] + R = [{}, pconst(Fr(1, 2) / alpha, nvr)] + deg_a = max(k[0] for k in poly) + for n in range(1, deg_a): + P.append(pscale(padd(pscale(P[n - 1], Fr(n)), pscale(pmul(beta, P[n]), Fr(-1))), Fr(1, 2) / alpha)) + R.append(pscale(padd(pscale(R[n - 1], Fr(n)), pscale(pmul(beta, R[n]), Fr(-1))), Fr(1, 2) / alpha)) + if MUT == "recurrence": + P[3] = pscale(P[3], Fr(101, 100)) + PE, PR = {}, {} + for k, c in poly.items(): + n = k[0] + w = {k[1:]: c} + PE = padd(PE, pmul(w, P[n])) + PR = padd(PR, pmul(w, R[n])) + # exponents: erfc branch -(rest - beta^2/(4 alpha)), polynomial branch -rest + expE = padd(rest, pscale(pmul(beta, beta), Fr(-1, 4) / alpha)) + return alpha, beta, PE, PR, expE, rest + + +def d2_terms(): + """m = 2: one remaining variable p; erfc branch e^{-p^2/8} erfc(p sqrt30/20), polynomial branch e^{-p^2/5}.""" + alpha, beta, PE, PR, expE, rest = a_integrated_terms(2) + if set(expE) - {(2,)} or set(rest) - {(2,)} or set(beta) - {(1,)}: + raise ValueError("m = 2 structure") + gE, gR, b1 = expE[(2,)], rest[(2,)], beta[(1,)] + kappa = Iv.frac(b1) / (2 * Iv.frac(alpha).sqrt()) # erfc argument slope + constE = Iv.frac(Fr(1, 2)) * (PI / Iv.frac(alpha)).sqrt() + pref = (Iv.frac(Fr(3, 5))).sqrt() / mehta_Z(2) * 2 # sqrt(3/(m+3)) . (m! sectors) / Z_m + pe, pr = univariate(PE), univariate(PR) + tE = (pref * constE, [Gauss(gE), ErfLin(kappa, -1)], [(Poly(pe), None)], [2 * abs(c) for c in pe], gE) + tR = (pref, [Gauss(gR)], [(Poly(pr), None)], [abs(c) for c in pr], gR) + return [tE, tR] + + +def d3_terms(): + """m = 3: remaining (p, q); substitute p = (s - q)/2 (Jacobian 1/2), integrate q over (0, s) exactly through J_k.""" + alpha, beta, PE, PR, expE, rest = a_integrated_terms(3) + s, q = pvar(0, 2), pvar(1, 2) + p_img = pscale(padd(s, pscale(q, Fr(-1))), Fr(1, 2)) + images = [p_img, q] + PEs, PRs = psubst(PE, images, 2), psubst(PR, images, 2) + eE, eR, bs = psubst(expE, images, 2), psubst(rest, images, 2), psubst(beta, images, 2) + # the decoupling: exponents must be diagonal in (s, q) with q^2 coefficient 1/8, and beta must be a multiple of s + if set(eE) != {(2, 0), (0, 2)} or set(eR) != {(2, 0), (0, 2)} or eE[(0, 2)] != Fr(1, 8) or eR[(0, 2)] != Fr(1, 8) or set(bs) != {(1, 0)}: + raise ValueError("decoupling failed: " + repr((eE, eR, bs))) + gE, gR = eE[(2, 0)], eR[(2, 0)] + kappa = Iv.frac(bs[(1, 0)]) / (2 * Iv.frac(alpha).sqrt()) + constE = Iv.frac(Fr(1, 2)) * (PI / Iv.frac(alpha)).sqrt() + pref = (Iv.frac(Fr(1, 2))).sqrt() / mehta_Z(3) * 6 / 2 # sqrt(3/6) . 3! / Z_3 . Jacobian 1/2 + if MUT == "erfc-branch": + constE = constE * Iv.frac(Fr(1, 2)) + branches = [] + def by_q_power(Pxy): + out = {} + for (i, j), c in Pxy.items(): + out.setdefault(j, {})[(i,)] = c + return out + for branch, Pxy, g, const, extra in (("E", PEs, gE, pref * constE, [ErfLin(kappa, -1)]), ("R", PRs, gR, pref, [])): + pairs, tail_poly = [], {} + for k, coef_poly in sorted(by_q_power(Pxy).items()): + ek = univariate(coef_poly) + jk = Jk(k) + pairs.append((Poly(ek), jk)) + # tail bound for s >= T: |J_k(s)| <= J_k(inf), erfc <= 2; rational upper bound of J_k(inf) from its interval + jk_bound = Fr(jk.infinity().hi) * (2 if branch == "E" else 1) + for j, c in enumerate(ek): + tail_poly[j] = tail_poly.get(j, Fr(0)) + abs(c) * jk_bound + tp = [tail_poly.get(j, Fr(0)) for j in range(max(tail_poly) + 1)] + branches.append((const, [Gauss(g)] + extra, pairs, tp, g)) + return branches + + +# ---------------- floating controls +def float_orthant(m, T=14.0, n=None): + n = n or {1: 20000, 2: 400, 3: 90}[m] + coef = math.sqrt(3 / (m + 3)) / float(mehta_Z_float(m)) + h = T / n + grid = [-(i + 0.5) * h for i in range(n)] + total = 0.0 + for mu in itertools.product(grid, repeat=m): + S = sum(mu) + w = math.prod(x * x for x in mu) + for i in range(m): + for j in range(i + 1, m): + w *= abs(mu[i] - mu[j]) + total += w * math.exp(-sum(x * x for x in mu) / 4 + S * S / (4 * (m + 3))) + return coef * total * h ** m + + +def mehta_Z_float(m): + return 2 ** (m / 2 + m * (m - 1) / 4) * (2 * math.pi) ** (m / 2) * math.prod(math.gamma(1 + j / 2) for j in range(1, m + 1)) / math.gamma(1.5) ** m + + +def float_Z3(T=14.0, n=100): + h = 2 * T / n + grid = [-T + (i + 0.5) * h for i in range(n)] + tot = 0.0 + for l1 in grid: + for l2 in grid: + for l3 in grid: + tot += abs((l1 - l2) * (l1 - l3) * (l2 - l3)) * math.exp(-(l1 * l1 + l2 * l2 + l3 * l3) / 4) + return tot * h ** 3 + + +def float_D3_pq(T=30.0, n=1200): + """D_3 by numerical (p, q) quadrature after the exact a-integration only: independent of the s-substitution, the + q-integration and J_k; midpoint rule, error O(h^2).""" + alpha = 3 / 8 + def G_list(beta, nmax=6): + G0 = 0.5 * math.sqrt(math.pi / alpha) * math.exp(beta * beta / (4 * alpha)) * math.erfc(beta / (2 * math.sqrt(alpha))) + G = [G0, (1 - beta * G0) / (2 * alpha)] + for k in range(1, nmax): + G.append((k * G[k - 1] - beta * G[k]) / (2 * alpha)) + return G + def pm(u, v): + out = [0.0] * (len(u) + len(v) - 1) + for i, x in enumerate(u): + for j, y in enumerate(v): + out[i + j] += x * y + return out + h = T / n + tot = 0.0 + for i in range(n): + p = (i + 0.5) * h + for j in range(n): + q = (j + 0.5) * h + G = G_list(p / 2 + q / 4) + full = pm(pm([0, 0, 1.0], [p * p, 2 * p, 1.0]), [(p + q) ** 2, 2 * (p + q), 1.0]) + tot += sum(full[k] * G[k] for k in range(7)) * p * q * (p + q) * math.exp(-(p * p / 3 + 5 * q * q / 24 + p * q / 3)) + return tot * h * h / (16 * math.pi) + + +# ---------------- main +def main(): + global MUT + ap = argparse.ArgumentParser() + ap.add_argument("--mutant", choices=MUTANTS) + args = ap.parse_args() + MUT = args.mutant + t0 = time.time() + checks = {} + + D1 = enclose_integral(d1_terms(), T=Fr(40)) + D2 = enclose_integral(d2_terms(), T=Fr(40)) + T3 = d3_terms() + D3 = enclose_integral(T3) + D3_coarse = enclose_integral(T3, K=K_ORDER // 2, cell=Fr(1), remainder=MUT != "remainder-dropped") + if MUT == "remainder-dropped": + D3 = enclose_integral(T3, K=6, cell=Fr(1), remainder=False) + D2_exact = Iv.frac(Fr(29, 6)) - SQRT6 + + def c_ref(d, D): + return GAMMA76 * CBRT32 * sphere_area(d) * D / (SQRT3 * SQRTPI * ipow(2 * PI, d)) + c2, c3, c4 = c_ref(2, D1), c_ref(3, D2), c_ref(4, D3) + c4_closed = GAMMA76 * CBRT32 * D3 / (8 * SQRT3 * PI * PI * SQRTPI) + + # torus transfer for d = 4, L = 24 (NOTE section 4): image bound, covariance sandwich, density and cone-moment ratios + shells = 5120 if MUT != "image-shells" else 2560 + E4 = Iv.frac(2 * shells * (76 * 24 ** 6 + 15)) * (-Iv.frac(288)).exp() + n_jet = 4 + 1 + 10 # dim of (G, t, svec H) in d = 4 + eps = Iv.frac(6 * n_jet) * E4 # spectral norm <= n . 2E; C_ref >= I/3 + b_exp = 7 if MUT != "sandwich-exponent" else 5 # d + n_A/2 with n_A = 6; a = 20/3 <= 7 + # delta = (1+eps)^7/(1-eps)^b - 1 = [sum_k (C(7,k) - (-1)^k C(b,k)) eps^k] / (1-eps)^b, evaluated without cancellation; + # the ratio c_{4,24}/c_{4,ref} lies in [1/(1+delta), 1+delta] (the lower bound (1-eps)^7/(1+eps)^b exceeds 1/(1+delta)) + num = Iv(ZERO) + for k in range(1, 8): + num = num + Iv.frac(math.comb(7, k) - (-1) ** k * math.comb(b_exp, k)) * ipow(eps, k) + delta = num / ipow(1 - eps, b_exp) + ratio = Iv((1 / (1 + delta)).lo, (1 + delta).hi) + c4_torus = c4 * ratio + rel_bound = delta.hi + + # ---- rules + f1, f2, f3, f3pq = float_orthant(1), float_orthant(2), float_orthant(3), float_D3_pq() + def near(iv, fv, tol): + f, t = Decimal(repr(fv)), Decimal(repr(tol)) + return CE.subtract(iv.lo, f) <= t and CE.subtract(f, iv.hi) <= t + z3f = float_Z3() + checks["FLOAT_INSIDE"] = bool(near(D1, f1, 1e-6) and near(D2, f2, 2e-3) and near(D3, f3, 2e-2) and near(D3, f3pq, 5e-4) + and near(mehta_Z(3), z3f, 1.5)) + checks["D1_EXACT"] = bool(D1.contains(Fr(4, 3))) + checks["D2_EXACT"] = bool(D2.intersects(D2_exact)) + checks["WIDTHS"] = bool(D1.width() < Decimal("1e-40") and D2.width() < Decimal("1e-40") and D3.width() < Decimal("1e-30") + and c4.width() < Decimal("1e-30")) + side = {d: Iv(Decimal(a), Decimal(b)) for d, (a, b) in SIDE24.items()} + checks["SIDE24_CONSISTENT"] = bool(all(c.intersects(side[d] + Iv(-side[d].hi * Decimal("1e-106"), side[d].hi * Decimal("1e-106"))) + for d, c in ((2, c2), (3, c3)))) + checks["CLOSED_FORM_D4"] = bool(c4.intersects(c4_closed)) + checks["TRUNCATION_NESTING"] = bool(D3_coarse.intersects(D3) and D3_coarse.width() >= D3.width()) + checks["MONOTONE"] = bool(D1.hi < D2.lo < D2.hi < D3.lo and c2.lo > c3.hi > c3.lo > c4.hi and c4.lo > 0) + checks["IMAGE_BOUND"] = bool(Decimal("1.2e-111") < E4.lo and E4.hi < Decimal("1.3e-111")) + # proving direction (v1.2, Codex 4148778791 on Math-#200): 13 eps <= delta <= 15 eps for every represented value; the + # enclosures of delta and 14 eps intersect (delta - 14 eps ~ 98 eps^2 is far below the interval resolution: recorded, not certified) + checks["TRANSFER_BOUND"] = bool(delta.lo >= CC.multiply(Decimal(13), eps.hi) and delta.hi <= CF.multiply(Decimal(15), eps.lo) + and delta.hi >= CF.multiply(Decimal(14), eps.lo) and delta.lo <= CC.multiply(Decimal(14), eps.hi) + and c4_torus.contains(c4)) + # Mehta's constant at m = 3 by the same decoupling (NOTE section 1): Z_3 = (3/2) sqrt(pi/3) int_0^inf e^{-s^2/24} (4 s^2 - 32 + 32 e^{-s^2/8}) ds = 48 sqrt2 pi + zc = Iv.frac(Fr(3, 2)) * (PI / 3).sqrt() + z3_quad = enclose_integral([(zc, [Gauss(Fr(1, 24))], [(Poly([Fr(-32), Fr(0), Fr(4)]), None)], [Fr(32), Fr(0), Fr(4)], Fr(1, 24)), + (zc * 32, [Gauss(Fr(1, 6))], [(Poly([Fr(1)]), None)], [Fr(1)], Fr(1, 6))]) + checks["MEHTA_Z3"] = bool(z3_quad.intersects(mehta_Z(3)) and z3_quad.intersects(48 * SQRT2 * PI) and z3_quad.width() < Decimal("1e-40")) + checks["LIBRARY_EXACT"] = bool(library_exact()) + pinned = {"D_3": D3, "c_4_ref": c4, "pi": PI} + checks["PINNED"] = bool(all(str(pinned[k].lo).startswith(p) and str(pinned[k].hi).startswith(p) for k, p in PINNED.items())) + if len(checks) != 13: + raise ValueError("rule count changed") + passed = all(checks.values()) + + out = { + "object": "CL-SIDE24-D4-COEFFICIENT-20260930-v1", + "scientific_effect": "NONE", + "certified": True, + "mutant": MUT, + "passed": passed, + "checks": checks, + "cone_moments": {"D_1": D1.pair(), "D_2": D2.pair(), "D_2_exact_29/6-sqrt6": D2_exact.pair(), "D_3": D3.pair(), + "D_3_coarse": D3_coarse.pair()}, + "coefficients": {"c_2_ref": c2.pair(), "c_3_ref": c3.pair(), "c_4_ref": c4.pair(), "c_4_ref_closed": c4_closed.pair(), + "c_4_24": c4_torus.pair()}, + "transfer": {"E_4": E4.pair(), "epsilon": eps.pair(), "delta": delta.pair(), "n_jet": n_jet, + "exponents": {"a": "20/3 (<= 7 used)", "b": 7}}, + "constants": {"pi": PI.pair(), "Gamma(7/6)": GAMMA76.pair(), "(3/2)^(1/3)": CBRT32.pair(), + "Z_1": mehta_Z(1).pair(), "Z_2": mehta_Z(2).pair(), "Z_3": mehta_Z(3).pair(), "Z_3_quadrature": z3_quad.pair(), + "|S^3|": sphere_area(4).pair()}, + "quadrature": {"K": K_ORDER, "cell": str(CELL), "T": str(T_CUT), "rho": str(RHO)}, + "float_reference_10sig": {"D_1_orthant": "%.10g" % f1, "D_2_orthant": "%.10g" % f2, "D_3_orthant": "%.10g" % f3, + "D_3_pq_quadrature": "%.10g" % f3pq, "Z_3_cube": "%.10g" % z3f}, + "exact_relations": { + "D_m": "sqrt(3/(m+3)) Z_m^-1 int_{mu<0} prod mu_i^2 |Delta(mu)| exp(-|mu|^2/4 + (sum mu)^2/(4(m+3))) dmu", + "Z_m": "2^(m/2 + m(m-1)/4) (2 pi)^(m/2) prod_{j<=m} Gamma(1 + j/2) / Gamma(3/2)^m; Z_3 = (3/2) sqrt(pi/3) int_0^inf e^{-s^2/24} (4 s^2 - 32 + 32 e^{-s^2/8}) ds = 48 sqrt2 pi", + "m=3 sector": "mu_3 = -a, mu_2 = -a-p, mu_1 = -a-p-q; Q = 3a^2/8 + p^2/3 + 5q^2/24 + ap/2 + aq/4 + pq/3", + "a-integral": "G_0 = (1/2) sqrt(8 pi/3) exp(p^2/6 + pq/6 + q^2/24) erfc((2p+q)/(2 sqrt6)); G_1 = (1 - beta G_0)/(2 alpha); G_{n+1} = (n G_{n-1} - beta G_n)/(2 alpha)", + "decoupling": "s = 2p + q: erfc branch exp(-s^2/24 - q^2/8), polynomial branch exp(-s^2/12 - q^2/8); dp dq = ds dq/2, 0 < q < s", + "q-integral": "J_0 = sqrt(2 pi) erf(s/(2 sqrt2)), J_1 = 4(1 - e^{-s^2/8}), J_{k+2} = -4 s^{k+1} e^{-s^2/8} + 4(k+1) J_k", + "D_3": "I_3 / (16 pi)", + "c_d_ref": "Gamma(7/6) (3/2)^(1/3) |S^{d-1}| D_{d-1} / (sqrt3 sqrt pi (2 pi)^d); d = 4: Gamma(7/6) (3/2)^(1/3) D_3 / (8 sqrt3 pi^(5/2))", + "transfer": "E_4 = 10240 (76 24^6 + 15) e^-288; eps = 90 E_4; c_{4,24}/c_{4,ref} in [(1-eps)^7/(1+eps)^7, (1+eps)^7/(1-eps)^7]", + }, + } + sys.stdout.write(json.dumps(out, indent=1) + "\n") + sys.stderr.write("seconds %.1f\n" % (time.time() - t0)) + sys.exit(0 if passed else 1) + + +def library_exact(): + ok = True + for q in (Fr(3), Fr(120), Fr(1, 2), Fr(288), Fr(5, 7), Fr(48) - Fr(51035039649098175696, 10 ** 29)): + g = -q / 64 + parts, term, ssum = [], Fr(1), Fr(0) + for k in range(0, 60): + ssum += term + parts.append(ssum) + term = term * g / (k + 1) + lo, hi = min(parts[-1], parts[-2]), max(parts[-1], parts[-2]) + exact = Iv(Iv.frac(lo ** 64).lo, Iv.frac(hi ** 64).hi) + X = Iv.frac(q) + got = (-X).exp() + ok &= got.intersects(exact) and CC.divide(got.width(), got.lo) < Decimal(10) ** (-PREC + 3) + r = root_frac(q, 2) + ok &= X.sqrt().intersects(r) and CC.divide(X.sqrt().width(), r.lo) < Decimal(10) ** (-PREC + 3) + ok &= (-X).lo == X.hi.copy_negate() and (-X).hi == X.lo.copy_negate() + return ok + + +if __name__ == "__main__": + main()