diff --git a/.github/workflows/c2-exact.yml b/.github/workflows/c2-exact.yml new file mode 100644 index 000000000..b8574323d --- /dev/null +++ b/.github/workflows/c2-exact.yml @@ -0,0 +1,79 @@ +name: c2 exact closed forms (replay) + +'on': + pull_request: + paths: + - 'frontiers/c2_exact_20261001/**' + - '.github/workflows/c2-exact.yml' + - 'frontiers/three_fronts_20260924/LIFETIME_REMAINDER.md' + - 'coefficients/side24_v1/PROOF.md' + workflow_dispatch: + +permissions: + contents: read + +jobs: + exact-replay: + runs-on: ubuntu-latest + timeout-minutes: 30 + 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/c2_exact_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: Exact derivation, checks and replay against RESULTS.json, mutants, controls + env: + PYTHONDONTWRITEBYTECODE: '1' + FLAGS: ${{ matrix.flags }} + run: | + set -euo pipefail + cd frontiers/c2_exact_20261001 + python $FLAGS certificate.py --check + for m in target finite-part cone-sign truncation sphere; 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 + python $FLAGS controls.py + cd "$GITHUB_WORKSPACE" + git diff --exit-code + test -z "$(git status --porcelain)" diff --git a/frontiers/c2_exact_20261001/NOTE.md b/frontiers/c2_exact_20261001/NOTE.md new file mode 100644 index 000000000..a2798b876 --- /dev/null +++ b/frontiers/c2_exact_20261001/NOTE.md @@ -0,0 +1,228 @@ +# The third-order coefficient `c2` in closed form: Gaussian kernel, `d = 1, 2, 3` + +**Object:** CL-C2-EXACT-20261001-v1. +**Author:** Anthropic Claude (Claude Code session `session_015wNj8LPTKXsaT68G3DgPPh`), 1 October 2026. +**Disposition:** exact closed forms, by an exact rational derivation replayed with the standard library, and certified +interval enclosures of the coefficient `c2` that Math-#216 defines and Math-#218 (0.1) and Math-#220 use, together with +the leading coefficient `c` from the same pipeline. Author-side. Nonauthor reads at `b43d36a`: OpenAI Slice B ACCEPT +(review 5379878251) and Slice A PASS (review 5379979090) with amendment OA-223-A-01, applied in v1.1 (§§0–2: the analytic +fixed-cone surrogate is written `A~_r`); xAI Slice C arithmetic check (5932294410). Scope claim: Math-#216 comment +5931466669. +**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-#214, #216, #218 or #220 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-#216 ([D], head `4f0927f`) defines, from the two-point kernel `A_r(b, k, u)` of [R] and its fold expansion +`A_r = A_0 + r^2 A_2 + O(r^3)`, the coefficients below. Here `A_0`, `A_2` are the Taylor coefficients of the fixed-cone +surrogate `A~_r` of §1, which is real-analytic in `r` (§2.3); with the consumed comparison `A_r = A~_r + O(r^3)` (§1) they +are the coefficients of the fold expansion of `A_r` itself. No analyticity of the typed `A_r` is claimed. + + c = (1/3) |S^(d-1)| integral_R db integral_0^oo A_0(b, k) k^(-2/3) dk, + c2 = (1/3) |S^(d-1)| integral_R db f.p. integral_0^oo A_2(b, k) k^(-4/3) dk + = (1/3) |S^(d-1)| integral_R db integral_0^oo [A_2(b, k) - A_2(b, 0)] k^(-4/3) dk. + +`c` is the leading coefficient of the short-lifetime law. `c2` is the coefficient of `ell^(1/3)`: +- Math-#218 (Theorem T, candidate status) proves it for the candidate density; +- Math-#220 (Theorem E3) proves it for the elder density, conditional on Math-#187; +- [D] derives it formally and by Monte Carlo, and lists a certified value as not claimed. + +**Theorem (closed forms).** For the Gaussian kernel `K(z) = exp(-|z|^2/2)` on `R^d`: + +| `d` | `c` | `c2` | +|---|---|---| +| `1` | `Gamma(1/6) 12^(-1/6) / (6 pi^(3/2))` | `(3/4) 12^(1/6) Gamma(5/6) / pi^(3/2)` | +| `2` | `Gamma(1/6) 12^(-1/6) / (9 pi^(3/2))` | `(13/18) 12^(1/6) Gamma(5/6) / pi^(3/2)` | +| `3` | `(29/72 - sqrt6/12) Gamma(1/6) 12^(-1/6) / pi^(5/2)` | `(5/48)(33 - 7 sqrt6) 12^(1/6) Gamma(5/6) / pi^(5/2)` | + +In every dimension the ratio is an algebraic multiple of one transcendental factor: + + c2/c = kappa_d Gamma(5/6) / (12^(2/3) Gamma(1/6)), kappa_1 = 54, kappa_2 = 78, kappa_3 = 18 (141 - sqrt6)/25. + +Here `12^(1/6) = 2^(1/3) 3^(1/6)`. The certified values (interval arithmetic, outward to 16 digits, `RESULTS.json`) are: + +| `d` | `c` | `c2` | `c2 / c` | +|---|---|---|---| +| `1` | `[0.1101103789590410, 0.1101103789590625]` | `[0.2300445802661503, 0.2300445802661998]` | `[2.089217950577392, 2.089217950578207]` | +| `2` | `[0.07340691930602738, 0.07340691930604167]` | `[0.2215244106266632, 0.2215244106267110]` | `[3.017759261945122, 3.017759261946299]` | +| `3` | `[0.04177593184059426, 0.04177593184060272]` | `[0.1612340491269447, 0.1612340491269810]` | `[3.859496174547095, 3.859496174548665]` | + +At 30 digits (mpmath, a control only): `c2 = 0.230044580266174952088473542454`, `0.221524410626686990900011559`, +`0.161234049126962737547589515632`. + +**Consequences.** +- **[D]'s eight digits are confirmed and made exact.** [D] gives `c2 = 0.23004458`, `0.22152441`, `0.16123405` and + `c2/c = 2.089218`, `3.017759`, `3.859496` from an `r`-fit. Each is the correct rounding of the closed form (§7). +- **The leading coefficient comes out exactly.** In `d = 2, 3` it equals the closed form `c_(d,ref)` of + `coefficients/side24_v1` (1). In `d = 1` it equals Math-#214's `C_0`, and `c2` equals #214's `2 B_2`. These are + identities of exact numbers, checked by the certificate (§4). +- **The three-term law in closed form.** With [CU]'s `c1`, the relative correction is + `nu/(c ell^(-1/3)) - 1 = (c1/c) ell^(7/12) + (c2/c) ell^(2/3) + ...`. In `d = 3` it is + `-5.071062... ell^(7/12) + 3.859496174548... ell^(2/3)`, where `c1/c` is certified in Math-#219 and `c2/c` here. + +## 1. What is evaluated + +The kernel is that of [R] §§2–4 and [D] §0. Pins sit at `M = -r e1/2` and `S = r e1/2` with zero gradients and heights +`b`, `b - k r^3`. The symmetric pin vector `U_r` has target `v_r = (b - k r^3/2, -k r^2, 0, 12k, 0, ..., 0)`, with rows +as in [D]'s `c2_check.py` (`pseries.pin_rows`), and + + A_r(b, k) = 12 pi_r(v_r) E[ |det H_M| |det H_S| 1{typed} | U_r = v_r ] / r^2, u = e1 (isotropy). + +As in [D] §2.2, the typed indicator is: +- nothing in `d = 1`, where the typed event has probability `1 - O(e^(-c k^2/r^2))`; +- `1{f_yy(0) < 0}` in `d = 2`; +- `1{D_y^2 f(0) < 0}` in `d = 3`. + +Write `A~_r` for the kernel with the typed indicator replaced in this way: the **fixed-cone surrogate**. Each +replacement changes the kernel by `O(r^3)` at fixed admissible `(b, k)`, so + + A_r(b, k) = A~_r(b, k) + O(r^3), A~_r = A_0 + r^2 A_2 + O(r^3) (§2.3). + +The first relation is consumed from [D] §2.2. Everything computed below is `A~_r`; analyticity is proved for `A~_r` only, +and no claim is made that the typed `A_r` is analytic (OA-223-A-01). [T] Lemma F proves the fold expansion and identifies +`A_2(b, 0) = -12 pi_0 E_0[Y^2 1{A < 0} | b]`, which is (0.1). Both are consumed, not reviewed. On the typed set, +`|det H_M| |det H_S| = -det H_M det H_S`. + +## 2. Exact derivation (`pseries.py`, `exact.py`) + +**2.1 Series.** +- Every quantity is a linear form in the jets at 0 with coefficients in `Q[[r]]`: + `d^beta f(x e1) = sum_j x^j/j! d^(beta + j e1) f(0)`, with `x = -+r/2`. +- Jet covariances are the integers `(-1)^|b| d^(a+b) K(0)`. +- The pin rows divide by `r`, `r` and `r^3`. Each division is exact, because the low coefficients vanish identically; + the code asserts this. +- Keeping `r^0 .. r^10`, every derived series is exact through `r^7`. The integrand needs `r^6`. + +**2.2 Conditioning.** +- The `r^0` pin covariance is invertible (`det = 12` in `d = 1, 2`), so the conditional means `G v_r` and covariances + are computed exactly in `Q[[r]]`. +- In `d >= 2` the Hessian entries are conditioned further on the transverse block `T = D_y^2 f(0)`: on `a = x1` in + `d = 2`, and on `T = diag(x1, x2)` in `d = 3`. +- `d = 3` uses transverse rotation invariance. The pins lie on `e1`, and transverse rotations preserve the law and + `det H_M det H_S`. With the Jacobian of the symmetric `2 x 2` eigenvalue map, + `integral_(T<0) g(T) p(T) dT = pi integral_(x2 b - L/(2 a_b)` leaves Gaussian moments `integral b^n e^(-a_b b^2) db` and the cone exponent `R - L^2/(4 a_b)`. + +**3.3 The cone.** +- **`d = 2`:** half-line moments `integral_(-oo)^0 a^n e^(-g a^2) da`, with `g = 3/16`. +- **`d = 3`:** with `x1 = -s` and `x2 = -s - t`, the exponent is `(3/10) s^2 + (3/10) s t + (1/5) t^2` and + `D = 4 al ga - be^2 = 3/20`. The quadrant moments `M(i, j)` follow from two integration-by-parts identities, in `s` + and in `t`. These form a `2 x 2` system with determinant `D`, starting from + `M(0, 0) = theta = (pi/2 - arctan(be/sqrt D))/sqrt D`. +- Every `M(i+1, j)` with `j >= 1` is reached twice, and the two values agree exactly (asserted). +- **The angle `theta` cancels exactly in both `c` and `c2`.** The results lie in `Q(sqrt6)` times powers of `pi`. + +Numbers are carried as `sum q pi^(h/2) sqrt(m) theta^e`, with `q` rational and `m` squarefree. + +## 4. Exact checks (`certificate.py`) + +1. **`theta` does not occur** in any result. +2. **`d = 2, 3`, `c` equals `coefficients/side24_v1` (1).** Since `Gamma(7/6) = Gamma(1/6)/6` and + `(3/2)^(1/3) = sqrt3 12^(-1/6)`, the closed form `Gamma(7/6)(3/2)^(1/3) D_(d-1)/(2 sqrt3 pi^(d-1) sqrt pi)` equals + `Gamma(1/6) 12^(-1/6) D_(d-1)/(12 pi^(d-1/2))`, with `D_1 = 4/3` and `D_2 = 29/6 - sqrt6`. Equality is exact. +3. **`d = 1`, `c = C_0` and `c2 = 2 B_2` of Math-#214 (D1.2).** For the Gaussian spectral moments `1, 3, 15, 105`: + - `D = 6`, `sigma_3 = sqrt6`, `p_12 = 1/(2 pi sqrt3)` and `Q/(120 lambda_2 lambda_4 D) = 3/4`; + - `72^(-1/6) 6^(2/3) = sqrt6 12^(-1/6)` and `3^(1/3) 6^(1/3) = 12 sqrt3 12^(-5/6)` give `C_0 = atom/(6 pi^(3/2))` and + `2 B_2 = 9 atom'/pi^(3/2)`; + - both are exact equalities, and the two reductions are also checked against #214's formulas in floating point. +4. **Truncation stability.** The whole derivation is repeated with series order 12 and polynomial order 8, and every + result is identical. +5. **The ratios.** `kappa_d Gamma(5/6)/(12^(2/3) Gamma(1/6))` must meet the interval quotient of the enclosures. + +## 5. Evaluation + +- `Gamma(1/6)` and `Gamma(5/6)` come from Stirling's series at `q + 10`, with the remainder bounded by the first + omitted term, and the exact product (`gam.py`: the code of Math-#219's `consts.py`). +- `12^(rational)` is `exp(log)`; `pi` and `sqrt` come from `ia.py` and `elem.py` (Math-#219, byte-identical; `ia.py` is + also Math-#217's). +- Outward one-ulp rounding is used throughout. Publication widens by `1e-14` relative and rounds outward to 16 digits. +- `--check` re-derives everything and requires the regenerated document to equal `RESULTS.json`. + +## 6. Verification + +- `python3 -B -S certificate.py --check` takes about one second. It runs the exact derivation, the checks of §4 and the + evaluation, then compares with `RESULTS.json`. +- Five seeded defects must each exit 1: + - `target`: the pin target `11 k` instead of `12 k`; + - `finite-part`: the `l = 0` term without the subtraction of `A_2(b, 0)`; + - `cone-sign`: the sign of the cone cross term; + - `truncation`: polynomial order 3; + - `sphere`: `|S^1| = pi`. +- The workflow `.github/workflows/c2-exact.yml` runs, for `-B -S` and `-B -O -S`: the manifest and pin check, + `--check`, the mutants, the controls, and a clean-tree check. + +## 7. Floating controls (`controls.py`; not part of the certificate) + +1. **A finite-`r` kernel.** + - The fixed-cone surrogate `A~_r(b, k)` is computed at finite `r` directly from the closed-form kernel derivatives at + the displaced points, + `D^g K(z) = prod (-1)^(g_i) He_(g_i)(z_i) e^(-z_i^2/2)`, with 60-digit conditioning and no Taylor jets. + - The Richardson estimate `2 D(r) - D(2r)`, with `D(r) = (A~_r - A_0)/r^2`, converges to the exact `A_2(b, k)` at the + rate `O(r^2)`, at five points in `d = 1, 2, 3`. + - The relative errors at `r = 0.01, 0.005, 0.0025` are `1.3e-4`, `3.2e-5`, `8.0e-6` (`d = 1`), with ratios `3.75` to + `4.17` per halving everywhere. +2. **[D]'s floating values.** All nine (`c`, `c2`, `c2/c` for `d = 1, 2, 3`) are correct roundings of the certified + midpoints. The largest difference, `4.0e-13`, is for `c(d = 3)`, against a half-unit of `5e-13`. + +## 8. What this does not do + +- **Gaussian kernel only.** The periodized covariance changes the jets by `1 + O(e^(-L^2/8))` ([D] §0). No torus + transfer of `c2` is made. The finite part and the typed indicator do not fit the homogeneity argument of Math-#219's + Lemma S. +- **No `d >= 4`.** The method extends; the cone then needs the 3D eigenvalue reduction. +- **The definition is consumed.** + - The identification of `c2` as the `ell^(1/3)` coefficient is [D]'s (formal), [T]'s (candidate density, candidate + status) and Math-#220's (elder density, conditional on #187), and is not reviewed here. + - So is the `O(r^3)` replacement of the typed indicator, `A_r = A~_r + O(r^3)` (§1). Analyticity is proved for the + surrogate `A~_r` only. +- Not a review of any source. Same GitHub account as every lane; zero organizational-independence credit. + **I will not merge.** + +## 9. Provenance + +- **Sources** (`SOURCE_MAP.json`): + - [R] `frontiers/three_fronts_20260924/LIFETIME_REMAINDER.md` (the two-point kernel); + - [SIDE24] `coefficients/side24_v1/PROOF.md` (1). + + Both are on `main` and verified by the workflow. +- **Unmerged, recorded and not checked on this tree:** + - [D] Math-#216 `NOTE.md` and `c2_check.py` at head `4f0927f`; + - [T] Math-#218 `PROOF.md` at head `7f09d34`; + - [D1] Math-#214 `PROOF.md` at head `f4a58df`. +- **Files:** + - the certificate: `certificate.py`, `exact.py`, `pseries.py`, `gam.py`, `ia.py`, `elem.py`; + - `controls.py`; + - `RESULTS.json`, `SOURCE_MAP.json`, `SOURCE_FILES.json`. diff --git a/frontiers/c2_exact_20261001/README.md b/frontiers/c2_exact_20261001/README.md new file mode 100644 index 000000000..9fb11719d --- /dev/null +++ b/frontiers/c2_exact_20261001/README.md @@ -0,0 +1,14 @@ +# The third-order coefficient `c2` in closed form (Math-#216 / #218 / #220), Gaussian kernel + +`CL-C2-EXACT-20261001-v1`. Author-side exact closed forms with certified enclosures; scientific effect NONE. See `NOTE.md`. + +| `d` | `c2` | value | +|---|---|---| +| `1` | `(3/4) 12^(1/6) Gamma(5/6) / pi^(3/2)` (= Math-#214's `2 B_2`) | `0.2300445802661...` | +| `2` | `(13/18) 12^(1/6) Gamma(5/6) / pi^(3/2)` | `0.2215244106266...` | +| `3` | `(5/48)(33 - 7 sqrt6) 12^(1/6) Gamma(5/6) / pi^(5/2)` | `0.1612340491269...` | + +`c2/c = kappa_d Gamma(5/6)/(12^(2/3) Gamma(1/6))`, `kappa = 54, 78, 18(141 - sqrt6)/25`. + +Replay: `python3 -B -S certificate.py --check` (exact derivation, exact checks, interval evaluation; about one second). +Controls: `python3 -B -S controls.py`. Standard library only. diff --git a/frontiers/c2_exact_20261001/RESULTS.json b/frontiers/c2_exact_20261001/RESULTS.json new file mode 100644 index 000000000..3077d9b29 --- /dev/null +++ b/frontiers/c2_exact_20261001/RESULTS.json @@ -0,0 +1,85 @@ +{ + "object": "CL-C2-EXACT-20261001-v1", + "scientific_effect": "NONE", + "statement": "closed forms and interval enclosures (outward decimals) of the leading coefficient c and the third-order coefficient c2 of the short-lifetime law, Gaussian kernel, d = 1, 2, 3 (Math-#216, #218 (0.1)); see NOTE.md", + "exact": { + "d=1": { + "c": "[1/6 pi^(-3/2)] * Gamma(1/6) 12^(-1/6)", + "c2": "[9 pi^(-3/2)] * Gamma(5/6) 12^(-5/6)", + "c2 / c": "54 * Gamma(5/6)/(12^(2/3) Gamma(1/6))" + }, + "d=2": { + "c": "[1/9 pi^(-3/2)] * Gamma(1/6) 12^(-1/6)", + "c2": "[26/3 pi^(-3/2)] * Gamma(5/6) 12^(-5/6)", + "c2 / c": "78 * Gamma(5/6)/(12^(2/3) Gamma(1/6))" + }, + "d=3": { + "c": "[29/72 pi^(-5/2) + -1/12 sqrt(6) pi^(-5/2)] * Gamma(1/6) 12^(-1/6)", + "c2": "[165/4 pi^(-5/2) + -35/4 sqrt(6) pi^(-5/2)] * Gamma(5/6) 12^(-5/6)", + "c2 / c": "18 (141 - sqrt6)/25 * Gamma(5/6)/(12^(2/3) Gamma(1/6))" + } + }, + "values": { + "d=1": { + "c": [ + "0.1101103789590410", + "0.1101103789590625" + ], + "c2": [ + "0.2300445802661503", + "0.2300445802661998" + ], + "c2 / c": [ + "2.089217950577392", + "2.089217950578207" + ] + }, + "d=2": { + "c": [ + "0.07340691930602738", + "0.07340691930604167" + ], + "c2": [ + "0.2215244106266632", + "0.2215244106267110" + ], + "c2 / c": [ + "3.017759261945122", + "3.017759261946299" + ] + }, + "d=3": { + "c": [ + "0.04177593184059426", + "0.04177593184060272" + ], + "c2": [ + "0.1612340491269447", + "0.1612340491269810" + ], + "c2 / c": [ + "3.859496174547095", + "3.859496174548665" + ] + } + }, + "derivation": { + "series_order_R": 10, + "polynomial_order_NR": 6, + "stability_orders": [ + 12, + 8 + ], + "d=3 cone exponent [alpha, beta, gamma, D]": [ + "3/10", + "3/10", + "1/5", + "3/20" + ], + "d=2 cone exponent gamma_a": "3/16" + }, + "publication": { + "margin_relative": 1e-14, + "digits": 16 + } +} diff --git a/frontiers/c2_exact_20261001/SOURCE_FILES.json b/frontiers/c2_exact_20261001/SOURCE_FILES.json new file mode 100644 index 000000000..cf4012408 --- /dev/null +++ b/frontiers/c2_exact_20261001/SOURCE_FILES.json @@ -0,0 +1,71 @@ +{ + "object": "CL-C2-EXACT-20261001-v1", + "files": [ + { + "path": "NOTE.md", + "bytes": 14480, + "sha256": "598fc86a7c6b21bad438ef087e5aa7521f0fae08e25fa4f40f9b2ba635846dbb", + "git_blob": "a2798b87674ed5fc0d2144697e77ccb9b185b4ab" + }, + { + "path": "README.md", + "bytes": 786, + "sha256": "b5fd88c0abf4123c9539a8141b4c77f0123b446385190c65f651ace29b3bb6a5", + "git_blob": "9fb11719d968fb9e88aaa76d6b28dc89301ebb71" + }, + { + "path": "RESULTS.json", + "bytes": 1891, + "sha256": "94aa67f11211881d60e93ee1782e9b9d829d246d8f75d73eb05b2212278cb831", + "git_blob": "3077d9b29bc2393fbae1dd2c8420ca1bcc3707dd" + }, + { + "path": "SOURCE_MAP.json", + "bytes": 3811, + "sha256": "930dc2fcb6d93ab341d8d3f6b07a160a423cf33935438ae52296ced39f46bae2", + "git_blob": "b10ad7300409f7e7478bdd4eb9e13cf925da495b" + }, + { + "path": "certificate.py", + "bytes": 10102, + "sha256": "a76b0d441c8a80bbc77be45f1445fe6905d19b9b38d48fe42c81d4a373bb60f1", + "git_blob": "96f1913f909118d921ea16161a1ff0885ba6aaef" + }, + { + "path": "controls.py", + "bytes": 12100, + "sha256": "0183606269cf39d8da034d71d37b9ac18a892d0276ad28d5be140531721c1f5b", + "git_blob": "1d2b05386f83c20debd145328c4887d9641d0ef6" + }, + { + "path": "elem.py", + "bytes": 9825, + "sha256": "65490f314fd731ced8d893d137ed5d7a516cf3cf36ba68cf77d7845805c55fb0", + "git_blob": "99ad3332b40daa1eee50106dbaa8b4b726cd5ae6" + }, + { + "path": "exact.py", + "bytes": 16099, + "sha256": "974099926619e7c85aea9fbccb5285e54862909cd079895f16e785583c884d6c", + "git_blob": "b06819802d1b90813cd172dc372aaac771ff6720" + }, + { + "path": "gam.py", + "bytes": 1596, + "sha256": "ac1e9214b6a1332e4ba9ca76ed602126d12c753e3e809dd6cf108e34eaac5060", + "git_blob": "e7707bf14aad66ff8aa40ec5d40bf7ab426d4139" + }, + { + "path": "ia.py", + "bytes": 8930, + "sha256": "6d048d58c6ab18f21c0609fc888ca73497a60e234217bb4f7a309a1502fccc9e", + "git_blob": "cf4bbfa16f577fe65f4ab5fb3621943fa63cf9fa" + }, + { + "path": "pseries.py", + "bytes": 6782, + "sha256": "d1faa96ca7e7df0051fbfe78318a17d0c78d608eb9f1ed2e3a3b83a6598a9994", + "git_blob": "fac3c562b33e626d477aefecf2b1e680b0ae5f98" + } + ] +} diff --git a/frontiers/c2_exact_20261001/SOURCE_MAP.json b/frontiers/c2_exact_20261001/SOURCE_MAP.json new file mode 100644 index 000000000..b10ad7300 --- /dev/null +++ b/frontiers/c2_exact_20261001/SOURCE_MAP.json @@ -0,0 +1,67 @@ +{ + "object": "CL-C2-EXACT-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": "[R]", + "path": "frontiers/three_fronts_20260924/LIFETIME_REMAINDER.md", + "ref": "main 124c37d", + "consumed": "sections 2-4: the symmetric pin vector U_r, the target v_r and the two-point kernel A_r(b,k,u) = 12 pi_r(v_r) Z_r / r^2", + "current_required": true, + "bytes": 17734, + "sha256": "380b7d0abdb0fe2de5a6564565af9560d3f1b1ce22fd5538927db0c52f4f3a4a", + "git_blob": "247b3ecf80bfbe896948d5d489b2d5842a81c481" + }, + { + "label": "[SIDE24]", + "path": "coefficients/side24_v1/PROOF.md", + "ref": "main 124c37d", + "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 (an exact check of the pipeline)", + "current_required": true, + "bytes": 10272, + "sha256": "c06daccc4ba4b9168522b9888b76a7d599934fc3b91bd753ee5d492262917769", + "git_blob": "44b66f04f89fcd87383b3603fa69f1feb64cdddd" + }, + { + "label": "[D]", + "path": "frontiers/third_order_coefficient_20261001/NOTE.md", + "ref": "Math-#216 head 4f0927f (unmerged)", + "consumed": "section 0: the definition of c and c2 (fold expansion, finite part); section 2.2: the typed indicator replaced by 1{T < 0} at O(r^3); section 0 table: the floating values (controls only)", + "current_required": false, + "bytes": 24186, + "sha256": "ae15c0ee173b9d01d258536a11780cc4b128c57811d1fb5db9db391ffc759ee5", + "git_blob": "aa078a9cbfc3aa859c089596892e34a8293831eb" + }, + { + "label": "[D-code]", + "path": "frontiers/third_order_coefficient_20261001/c2_check.py", + "ref": "Math-#216 head 4f0927f (unmerged)", + "consumed": "Pinned.pin_rows and Pinned.target: the exact pin rows and target reproduced in pseries.pin_rows", + "current_required": false, + "bytes": 25031, + "sha256": "0f8a742958ce9bfa94382d73748fe93769cba43079570d67180bfa9082a6e68a", + "git_blob": "41f745275db81f6332e791759db5858eec8c5365" + }, + { + "label": "[T]", + "path": "frontiers/candidate_third_order_20261001/PROOF.md", + "ref": "Math-#218 head 7f09d34 (unmerged)", + "consumed": "(0.1) and Lemma F: the fold expansion A_r = A_0 + r^2 A_2 + O(r^3) and the identification of A_2(b, 0) (cited)", + "current_required": false, + "bytes": 37220, + "sha256": "6c1279164499d4fe2fbeee244b40dafa2159758d20f04826cd67d7602ea97fe4", + "git_blob": "70ca57ef05199e75417d7d57928c1e1c93bf944e" + }, + { + "label": "[D1]", + "path": "frontiers/d1_third_order_law_20261001/PROOF.md", + "ref": "Math-#214 head f4a58df (unmerged)", + "consumed": "(D1.2): C_0 and B_2 in closed form for spectral moments (an exact check of the d = 1 pipeline)", + "current_required": false, + "bytes": 44378, + "sha256": "279ed88ff2debd570e6426aec53cc93aabc7a298731cd9b25c45667a2c11077b", + "git_blob": "873532b955dbf7684edc1a53ba169c7216f778ab" + } + ], + "note": "current_required entries are verified byte-for-byte on the checked-out tree by the workflow; the unmerged pins record the exact bytes read and are not checked against the tree." +} diff --git a/frontiers/c2_exact_20261001/certificate.py b/frontiers/c2_exact_20261001/certificate.py new file mode 100644 index 000000000..96f1913f9 --- /dev/null +++ b/frontiers/c2_exact_20261001/certificate.py @@ -0,0 +1,249 @@ +"""Certificate for CL-C2-EXACT-20261001: the third-order coefficient c2 of Math-#216 / #218 (0.1) for the Gaussian kernel, +d = 1, 2, 3, in closed form and as interval enclosures; also the leading coefficient c by the same pipeline. + + python3 -B -S certificate.py --write derive, check, write RESULTS.json + python3 -B -S certificate.py --check derive, check, compare with RESULTS.json (must be identical) + python3 -B -S certificate.py --mutant N a seeded defect; the exact checks must reject it (exit 1) + +Standard library only. The derivation is exact (Fractions; exact.py, pseries.py); only the final evaluation of the closed +forms uses outward-rounded binary64 intervals (ia.py, elem.py, gam.py).""" +import argparse +import json +import math +import os +import sys +from fractions import Fraction as Fr + +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +sys.dont_write_bytecode = True + +import pseries as P +import exact as X +from exact import Num +from ia import dn, up, mul, div, isqrt, from_frac, PI +import gam + +PACKET = 'CL-C2-EXACT-20261001-v1' +MUTANTS = ('target', 'finite-part', 'cone-sign', 'truncation', 'sphere') +DIGITS = 16 +MARGIN = 1e-14 + + +def install_mutant(name): + if name in ('target', 'finite-part', 'cone-sign'): + X.MUTANT = name + elif name == 'truncation': + X.NR = 3 + elif name == 'sphere': + X.SPHERE[2] = Num.pi_half(2) + elif name is not None: + raise SystemExit('unknown mutant: %s' % name) + + +# ------------------------------------------------------------------------------------------- derivation +def derive(): + out = {} + for d in (1, 2, 3): + data = X.kernel_data(d) + for which, nm in ((0, 'c'), (2, 'c2')): + val, atom, info = X.coefficient(d, which, data) + out[(d, nm)] = (val, atom, info) + return out + + +def side24_cref(d): + """c_(d,ref) of coefficients/side24_v1 (1) in the atom Gamma(1/6) 12^(-1/6): + 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. + With Gamma(7/6) = Gamma(1/6)/6 and (3/2)^(1/3) = sqrt3 12^(-1/6): c_(d,ref) = atom * D_(d-1) / (12 pi^(d - 1/2)).""" + D = {2: Num.rat(Fr(4, 3)), 3: Num.rat(Fr(29, 6)) - Num.sqrt(6)}[d] + return D * Fr(1, 12) * Num.pi_half(-(2 * d - 1)) + + +def d1_closed(): + """Math-#214 (D1.2) for the Gaussian kernel: lambda_2, 4, 6, 8 = 1, 3, 15, 105, D = lambda_2 lambda_6 - lambda_4^2 = 6, + sigma_3 = sqrt(D/lambda_2) = sqrt6, p_12 = 1/(2 pi sqrt(lambda_2 lambda_4)) = 1/(2 pi sqrt3), + Q/(120 lambda_2 lambda_4 D) = 1620/2160 = 3/4. + C_0 = 2 . 72^(-1/6) Gamma(7/6) (2 pi)^(-1/2) p_12 sigma_3^(4/3); 72^(-1/6) 6^(2/3) = sqrt6 . 12^(-1/6), so + C_0 = atom_c . 2 (1/6) sqrt6 (2 pi)^(-1/2) / (2 pi sqrt3) = atom_c . (1/6) pi^(-3/2). + 2 B_2 = 2 . 2^(1/2) 3^(1/3) Gamma(5/6) (2 pi)^(-1/2) p_12 sigma_3^(2/3) . 3/4; 3^(1/3) 6^(1/3) = 12 sqrt3 . 12^(-5/6), so + 2 B_2 = atom_c2 . 2 sqrt2 12 sqrt3 (3/4) (2 pi)^(-1/2) / (2 pi sqrt3) = atom_c2 . 9 pi^(-3/2).""" + return Num.rat(Fr(1, 6)) * Num.pi_half(-3), Num.rat(9) * Num.pi_half(-3) + + +def d1_float(): + """the same two numbers directly from #214's formulas in floating point (a check of the reductions above)""" + l2, l4, l6, l8 = 1.0, 3.0, 15.0, 105.0 + D = l2 * l6 - l4 * l4 + s3 = math.sqrt(D / l2) + p12 = 1 / (2 * math.pi * math.sqrt(l2 * l4)) + C0 = 2 * 72 ** (-1 / 6) * math.gamma(7 / 6) * (2 * math.pi) ** -0.5 * p12 * s3 ** (4 / 3) + Q = 4 * l2 ** 2 * l4 * l8 + 26 * l2 * l4 ** 2 * l6 - 5 * l2 ** 2 * l6 ** 2 - 25 * l4 ** 4 + B2 = 2 ** 0.5 * 3 ** (1 / 3) * math.gamma(5 / 6) * (2 * math.pi) ** -0.5 * p12 * s3 ** (2 / 3) * Q / (120 * l2 * l4 * D) + return C0, 2 * B2 + + +def checks(res): + """Exact checks. Returns a list of failure strings.""" + fails = [] + for key, (val, atom, info) in res.items(): + if any(e for (_, _, e) in val.t): + fails.append('%s: the cone atom theta does not cancel' % (key,)) + for d in (2, 3): + if not res[(d, 'c')][0] == side24_cref(d): + fails.append('d = %d: c differs from side24_v1 (1): %s' % (d, res[(d, 'c')][0].text())) + C0, B2x2 = d1_closed() + if not res[(1, 'c')][0] == C0: + fails.append('d = 1: c differs from Math-#214 C_0') + if not res[(1, 'c2')][0] == B2x2: + fails.append('d = 1: c2 differs from Math-#214 2 B_2') + fc, fb = d1_float() + if not (abs(fc - num_float(C0, 'c')) < 1e-14 and abs(fb - num_float(B2x2, 'c2')) < 1e-14): + fails.append('d = 1: the reductions of #214 closed forms disagree in floating point') + return fails + + +def stability(res): + """recompute with higher truncation orders; the exact results must be identical""" + oldR, oldNR = P.R, X.NR + try: + P.R, X.NR = 12, 8 + res2 = derive() + finally: + P.R, X.NR = oldR, oldNR + return [('%s changes with the truncation order' % (key,)) for key in res if not res[key][0] == res2[key][0]] + + +# ------------------------------------------------------------------------------------------- evaluation +ATOMS = {'c': (Fr(1, 6), Fr(-1, 6)), 'c2': (Fr(5, 6), Fr(-5, 6))} + + +def num_float(n, which): + g, e = ATOMS[which] + a = math.gamma(float(g)) * 12.0 ** float(e) + return a * math.fsum(float(v) * math.pi ** (h / 2) * math.sqrt(m) for (h, m, _), v in n.t.items()) + + +def ipow(x, n): + out = (1.0, 1.0) + for _ in range(abs(n)): + out = mul(out, x) + return out if n >= 0 else div((1.0, 1.0), out) + + +def num_iv(n, which): + g, e = ATOMS[which] + atom = mul(gam.gamma(g), gam.rpow((12.0, 12.0), e)) + tot = (0.0, 0.0) + for (h, m, _), v in sorted(n.t.items()): + term = from_frac(v) + term = mul(term, ipow(PI, h // 2)) + if h % 2: + term = mul(term, isqrt(PI)) + if m != 1: + term = mul(term, isqrt(from_frac(Fr(m)))) + tot = (dn(tot[0] + term[0]), up(tot[1] + term[1])) + return mul(atom, tot) + + +def _round_out(x, upward): + from decimal import Decimal, ROUND_FLOOR, ROUND_CEILING + dd = Decimal(x) + q = Decimal(1).scaleb(dd.adjusted() - DIGITS + 1) + return str(dd.quantize(q, rounding=ROUND_CEILING if upward else ROUND_FLOOR)) + + +def pub(iv): + lo, hi = iv[0] - MARGIN * abs(iv[0]), iv[1] + MARGIN * abs(iv[1]) + return [_round_out(dn(lo), False), _round_out(up(hi), True)] + + +RATIO_TEXT = {1: '54', 2: '78', 3: '18 (141 - sqrt6)/25'} + + +def document(res): + vals, exact = {}, {} + for d in (1, 2, 3): + c, c2 = num_iv(res[(d, 'c')][0], 'c'), num_iv(res[(d, 'c2')][0], 'c2') + vals['d=%d' % d] = {'c': pub(c), 'c2': pub(c2), 'c2 / c': pub(div(c2, c))} + exact['d=%d' % d] = {'c': '[%s] * %s' % (res[(d, 'c')][0].text(), res[(d, 'c')][1]), + 'c2': '[%s] * %s' % (res[(d, 'c2')][0].text(), res[(d, 'c2')][1]), + 'c2 / c': '%s * Gamma(5/6)/(12^(2/3) Gamma(1/6))' % RATIO_TEXT[d]} + info3 = res[(3, 'c2')][2] + return { + 'object': PACKET, + 'scientific_effect': 'NONE', + 'statement': 'closed forms and interval enclosures (outward decimals) of the leading coefficient c and the ' + 'third-order coefficient c2 of the short-lifetime law, Gaussian kernel, d = 1, 2, 3 (Math-#216, ' + '#218 (0.1)); see NOTE.md', + 'exact': exact, + 'values': vals, + 'derivation': {'series_order_R': P.R, 'polynomial_order_NR': X.NR, 'stability_orders': [12, 8], + 'd=3 cone exponent [alpha, beta, gamma, D]': [str(info3[k]) for k in ('alpha', 'beta', 'gamma', + 'D')], + 'd=2 cone exponent gamma_a': str(res[(2, 'c2')][2]['gamma_a'])}, + 'publication': {'margin_relative': MARGIN, 'digits': DIGITS}, + } + + +def ratio_check(res): + """c2/c in closed form against the interval quotient""" + fails = [] + base = div(mul(gam.gamma(Fr(5, 6)), (1.0, 1.0)), mul(gam.rpow((12.0, 12.0), Fr(2, 3)), gam.gamma(Fr(1, 6)))) + k = {1: from_frac(Fr(54)), 2: from_frac(Fr(78)), + 3: mul(from_frac(Fr(18, 25)), (dn(141 - isqrt((6.0, 6.0))[1]), up(141 - isqrt((6.0, 6.0))[0])))} + for d in (1, 2, 3): + a = mul(k[d], base) + q = div(num_iv(res[(d, 'c2')][0], 'c2'), num_iv(res[(d, 'c')][0], 'c')) + if a[1] < q[0] or q[1] < a[0]: + fails.append('d = %d: closed-form ratio %r disjoint from the quotient %r' % (d, a, q)) + return fails + + +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) + 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) + try: + res = derive() + except (ArithmeticError, ValueError) as e: + print('DERIVATION FAILURE', e) + return 1 + fails = checks(res) + ratio_check(res) + if not fails and a.mutant is None: + fails += stability(res) + print('exact checks: %d failures' % len(fails)) + for f in fails: + print('CHECK FAILURE', f) + if fails: + return 1 + for d in (1, 2, 3): + print('d = %d: c = %s ; c2 = %s' % (d, res[(d, 'c')][0].text() + ' * ' + res[(d, 'c')][1], + res[(d, 'c2')][0].text() + ' * ' + res[(d, 'c2')][1])) + if not (a.write or a.check): + return 0 + doc = document(res) + path = os.path.join(HERE, 'RESULTS.json') + if a.write: + with open(path, 'w') as fh: + json.dump(doc, fh, indent=1) + fh.write('\n') + print('wrote', path) + return 0 + with open(path) as fh: + old = json.load(fh) + if old != json.loads(json.dumps(doc)): + print('MISMATCH: the recomputed document differs from RESULTS.json') + return 1 + print(json.dumps({'object': PACKET, 'check': 'passed', 'scientific_effect': 'NONE'})) + return 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/frontiers/c2_exact_20261001/controls.py b/frontiers/c2_exact_20261001/controls.py new file mode 100644 index 000000000..1d2b05386 --- /dev/null +++ b/frontiers/c2_exact_20261001/controls.py @@ -0,0 +1,301 @@ +"""Floating controls for CL-C2-EXACT-20261001. These are not part of the certificate. + +1. Finite-r kernel, pointwise. The fixed-cone surrogate A~_r(b, k) (NOTE.md section 1) is computed at finite r + directly from the closed-form kernel derivatives at the displaced points M = -r e1/2, S = r e1/2, 0, without Taylor + jets or series: D^g K(z) = prod_i (-1)^(g_i) He_(g_i)(z_i) e^(-z_i^2/2), conditioning in Decimal (60 digits), then + floating Isserlis and quadrature. The Richardson combination + 2 D(r) - D(2r), D(r) = (A~_r - A_0)/r^2, must approach the exact A_2(b, k) of exact.py at the rate O(r^2) + (error ratio 4 per halving of r). +2. Math-#216's floating values (head 4f0927f, NOTE.md section 0 table and README) against the closed forms. + +Run: python3 -B controls.py (pure Python, a few seconds).""" +import json +import math +import os +import sys +from decimal import Decimal, getcontext + +HERE = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, HERE) +sys.dont_write_bytecode = True +import exact as X + +getcontext().prec = 60 +DEC = Decimal + + +# ------------------------------------------------------------------------------------------- finite-r kernel +def he(n, x): + a, b = DEC(1), x + if n == 0: + return a + for k in range(1, n): + a, b = b, x * b - k * a + return b + + +def dK(g, z): + out = DEC(1) + for gi, zi in zip(g, z): + out *= (-1) ** gi * he(gi, zi) * (-(zi * zi) / 2).exp() + return out + + +def cov(u, w): + """u, w: lists of (coef, point, alpha)""" + tot = DEC(0) + for cu, xu, au in u: + for cw, xw, aw in w: + z = tuple(p - q for p, q in zip(xu, xw)) + tot += cu * cw * (-1) ** sum(aw) * dK(tuple(p + q for p, q in zip(au, aw)), z) + return tot + + +def solve(A, B): + """A^-1 B (Decimal Gauss-Jordan), B a list of columns""" + n = len(A) + M = [A[i][:] + [col[i] for col in B] for i in range(n)] + for c in range(n): + p = max(range(c, n), key=lambda i: abs(M[i][c])) + M[c], M[p] = M[p], M[c] + piv = M[c][c] + M[c] = [x / piv for x in M[c]] + for i in range(n): + if i != c and M[i][c] != 0: + f = M[i][c] + M[i] = [x - f * y for x, y in zip(M[i], M[c])] + return [[M[i][n + j] for i in range(n)] for j in range(len(B))] + + +def det(A): + n = len(A) + M = [row[:] for row in A] + out = DEC(1) + for c in range(n): + p = max(range(c, n), key=lambda i: abs(M[i][c])) + if p != c: + M[c], M[p] = M[p], M[c] + out = -out + out *= M[c][c] + for i in range(c + 1, n): + f = M[i][c] / M[c][c] + M[i] = [x - f * y for x, y in zip(M[i], M[c])] + return out + + +def unit(d, *idx): + e = [0] * d + for i in idx: + e[i] += 1 + return tuple(e) + + +def structure(d, r): + r = DEC(r) + Mp = (-r / 2,) + (DEC(0),) * (d - 1) + Sp = (r / 2,) + (DEC(0),) * (d - 1) + O = (DEC(0),) * d + e0, ex = (0,) * d, unit(d, 0) + rows = [[(DEC('0.5'), Mp, e0), (DEC('0.5'), Sp, e0)], + [(-1 / r, Mp, e0), (1 / r, Sp, e0)], + [(-1 / r, Mp, ex), (1 / r, Sp, ex)], + [(6 / r ** 2, Mp, ex), (6 / r ** 2, Sp, ex), (12 / r ** 3, Mp, e0), (-12 / r ** 3, Sp, e0)]] + for t in range(1, d): + ey = unit(d, t) + rows += [[(DEC('0.5'), Mp, ey), (DEC('0.5'), Sp, ey)], [(-1 / r, Mp, ey), (1 / r, Sp, ey)]] + forms = {} + for nm, pt in (('M', Mp), ('S', Sp)): + for i in range(d): + for j in range(i, d): + forms[(nm, i, j)] = [(DEC(1), pt, unit(d, i, j))] + trans = [] + for i in range(1, d): + for j in range(i, d): + forms[('T', i, j)] = [(DEC(1), O, unit(d, i, j))] + trans.append(('T', i, j)) + p = len(rows) + Sm = [[cov(a, b) for b in rows] for a in rows] + names = list(forms) + CFP = {u: [cov(forms[u], row) for row in rows] for u in names} + G = {u: col for u, col in zip(names, solve(Sm, [CFP[u] for u in names]))} # S^-1 C_PF -> mean coefficients + C = {(u, w): cov(forms[u], forms[w]) - sum(G[u][q] * CFP[w][q] for q in range(p)) for u in names for w in names} + Sinv_cols = solve(Sm, [[DEC(int(i == j)) for i in range(p)] for j in range(p)]) + return {'p': p, 'G': G, 'C': C, 'Sinv': Sinv_cols, 'detS': det(Sm), 'names': names, 'trans': trans, 'r': r} + + +def gl(n, a, b): + """Gauss-Legendre nodes (floating, Newton) on [a, b]""" + xs, ws = [], [] + for i in range(1, n + 1): + x = math.cos(math.pi * (i - 0.25) / (n + 0.5)) + for _ in range(100): + p0, p1 = 1.0, x + for k in range(2, n + 1): + p0, p1 = p1, ((2 * k - 1) * x * p1 - (k - 1) * p0) / k + dp = n * (x * p1 - p0) / (x * x - 1) + dx = p1 / dp + x -= dx + if abs(dx) < 1e-16: + break + xs.append(0.5 * (a + b) + 0.5 * (b - a) * x) + ws.append(0.5 * (b - a) * 2 / ((1 - x * x) * dp * dp)) + return list(zip(xs, ws)) + + +# small float polynomials in the cone variables: dict (i, j) -> coefficient +def pmul(a, b): + out = {} + for e1, c1 in a.items(): + for e2, c2 in b.items(): + e = (e1[0] + e2[0], e1[1] + e2[1]) + out[e] = out.get(e, 0.0) + c1 * c2 + return out + + +def detdet_poly(d, aff, cp): + """E[det H_M det H_S | T] as a float polynomial in the cone variables, by Isserlis""" + import itertools + tot = {} + for pp in itertools.permutations(range(d)): + for qq in itertools.permutations(range(d)): + fac = [('M', min(i, pp[i]), max(i, pp[i])) for i in range(d)] + \ + [('S', min(i, qq[i]), max(i, qq[i])) for i in range(d)] + sg = X.perm_sign(pp) * X.perm_sign(qq) + for mt in X.matchings(list(range(2 * d))): + used = set(i for pr in mt for i in pr) + c = float(sg) + for (i, j) in mt: + c *= cp[(fac[i], fac[j])] + if c == 0.0: + continue + poly = {(0, 0): c} + for i in range(2 * d): + if i not in used: + poly = pmul(poly, aff[fac[i]]) + for e, v in poly.items(): + tot[e] = tot.get(e, 0.0) + v + return tot + + +def A_finite(d, st, b, k, nq=40): + """A~_r(b, k) = -12 pi_r(v) E[det H_M det H_S 1{T < 0} | U_r = v] / r^2 at finite r (the fixed-cone surrogate)""" + r = st['r'] + p = st['p'] + v = [DEC(b) - DEC(k) * r ** 3 / 2, -DEC(k) * r * r, DEC(0), 12 * DEC(k)] + [DEC(0)] * (p - 4) + q = sum(v[i] * st['Sinv'][j][i] * v[j] for i in range(p) for j in range(p)) + pi_r = float((2 * DEC(math.pi)) ** (-DEC(p) / 2) / st['detS'].sqrt() * (-q / 2).exp()) + names, trans = st['names'], st['trans'] + mean = {u: sum(st['G'][u][j] * v[j] for j in range(p)) for u in names} + hn = [u for u in names if u not in trans] + if d == 1: + cp = {(u, w): float(st['C'][(u, w)]) for u in hn for w in hn} + aff = {u: {(0, 0): float(mean[u])} for u in hn} + Z = detdet_poly(d, aff, cp)[(0, 0)] + return -12 * pi_r * Z / float(r) ** 2 + CTT = [[st['C'][(x, y)] for y in trans] for x in trans] + nt = len(trans) + Bc = solve(CTT, [[st['C'][(u, trans[i])] for i in range(nt)] for u in hn]) # columns: C_TT^-1 C_Tu + B = {u: col for u, col in zip(hn, Bc)} + cp = {(u, w): float(st['C'][(u, w)] - sum(B[u][i] * st['C'][(trans[i], w)] for i in range(nt))) for u in hn for w in hn} + mT = [mean[t] for t in trans] + # cone variables (x1, x2): d = 2: a = x1; d = 3: diag(x1, x2) + aff = {} + for u in hn: + base = mean[u] - sum(B[u][i] * mT[i] for i in range(nt)) + if d == 2: + aff[u] = {(0, 0): float(base), (1, 0): float(B[u][0])} + else: + aff[u] = {(0, 0): float(base), (1, 0): float(B[u][0]), (0, 1): float(B[u][2])} + poly = detdet_poly(d, aff, cp) + CTTinv = solve(CTT, [[DEC(int(i == j)) for i in range(nt)] for j in range(nt)]) + Ci = [[float(CTTinv[j][i]) for j in range(nt)] for i in range(nt)] + dT = float(det(CTT)) + mTf = [float(x) for x in mT] + norm = (2 * math.pi) ** (-nt / 2) / math.sqrt(dT) + + def dens(z): + qq = sum(z[i] * Ci[i][j] * z[j] for i in range(nt) for j in range(nt)) + return norm * math.exp(-qq / 2) + tot = 0.0 + if d == 2: + L = 14.0 + abs(b) + for x, w in gl(nq, -L, 0.0): + val = sum(c * x ** e[0] for e, c in poly.items()) + tot += w * val * dens([x - mTf[0]]) + else: + L = 14.0 + abs(b) + nodes = gl(nq, 0.0, L) + for s, ws in nodes: + for t, wt in nodes: + x1, x2 = -s, -s - t + val = sum(c * x1 ** e[0] * x2 ** e[1] for e, c in poly.items()) + tot += ws * wt * val * dens([x1 - mTf[0], -mTf[1], x2 - mTf[2]]) * math.pi * t + return -12 * pi_r * tot / float(r) ** 2 + + +def A_exact(d, data, which, b, k, nq=40): + """A_which(b, k) from exact.py's polynomial data, evaluated by floating quadrature over the cone""" + H = data['G'].rcoef(which + 2) + E = (data['q0'] + data['Q0']) * X.Fr(1, 2) + const = sum(float(v) * math.pi ** (h / 2) * math.sqrt(m) for (h, m, _), v in data['const'].t.items()) + + def ev(poly, x1, x2): + return sum(float(c) * b ** e[1] * k ** e[2] * x1 ** e[3] * x2 ** e[4] for e, c in poly.t.items()) + if d == 1: + return const * math.exp(-ev(E, 0, 0)) * ev(H, 0, 0) + tot = 0.0 + L = 14.0 + abs(b) + if d == 2: + for x, w in gl(nq, -L, 0.0): + tot += w * math.exp(-ev(E, x, 0)) * ev(H, x, 0) + else: + nodes = gl(nq, 0.0, L) + for s, ws in nodes: + for t, wt in nodes: + x1, x2 = -s, -s - t + tot += ws * wt * math.exp(-ev(E, x1, x2)) * ev(H, x1, x2) * t + return const * tot + + +QUOTED_216 = { # Math-#216 head 4f0927f, NOTE.md section 0 table + 1: {'c': '0.110110378959', 'c2': '0.23004458', 'c2/c': '2.089218'}, + 2: {'c': '0.073406919306', 'c2': '0.22152441', 'c2/c': '3.017759'}, + 3: {'c': '0.041775931841', 'c2': '0.16123405', 'c2/c': '3.859496'}, +} + + +def main(): + res = json.load(open(os.path.join(HERE, 'RESULTS.json')))['values'] + print('1. finite-r kernel against the exact A_2: relative error of 2 D(r) - D(2r) at r = 0.01, 0.005, 0.0025') + out, bad = [], [] + for d, pts in ((1, [(0.3, 0.7)]), (2, [(0.5, 0.4), (-0.8, 0.9)]), (3, [(0.6, 0.5), (-0.4, 0.3)])): + data = X.kernel_data(d) + for b, k in pts: + A0, A2 = A_exact(d, data, 0, b, k), A_exact(d, data, 2, b, k) + errs = [] + for r in ('0.01', '0.005', '0.0025'): + D1 = (A_finite(d, structure(d, r), b, k) - A0) / float(r) ** 2 + D2 = (A_finite(d, structure(d, str(2 * float(r))), b, k) - A0) / (2 * float(r)) ** 2 + errs.append(abs(2 * D1 - D2 - A2) / abs(A2)) + print(' d = %d, (b, k) = (%4.1f, %3.1f): A_2 = %+.10f; errors %.1e, %.1e, %.1e (ratios %.2f, %.2f)' + % (d, b, k, A2, errs[0], errs[1], errs[2], errs[0] / errs[1], errs[1] / errs[2])) + out.append(errs) + if not (errs[2] < 1e-4 and all(3.0 < errs[i] / errs[i + 1] < 5.0 for i in range(2))): + bad.append('finite-r control, d = %d, (b, k) = (%g, %g)' % (d, b, k)) + print('2. Math-#216 floating values against the certified enclosures:') + for d in (1, 2, 3): + for key, rk in (('c', 'c'), ('c2', 'c2'), ('c2/c', 'c2 / c')): + q = QUOTED_216[d][key] + lo, hi = (float(x) for x in res['d=%d' % d][rk]) + mid = 0.5 * (lo + hi) + half = 0.5 * 10 ** -len(q.split('.')[1]) + print(' d = %d %-5s quoted %-15s certified mid %.13f |difference| %.1e (half-unit %.0e)%s' + % (d, key, q, mid, abs(float(q) - mid), half, '' if abs(float(q) - mid) < half else ' NOT A ROUNDING')) + if not abs(float(q) - mid) < half: + bad.append('#216 value d = %d %s' % (d, key)) + print('controls: %d failures' % len(bad)) + return 1 if bad else 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/frontiers/c2_exact_20261001/elem.py b/frontiers/c2_exact_20261001/elem.py new file mode 100644 index 000000000..99ad3332b --- /dev/null +++ b/frontiers/c2_exact_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/c2_exact_20261001/exact.py b/frontiers/c2_exact_20261001/exact.py new file mode 100644 index 000000000..b06819802 --- /dev/null +++ b/frontiers/c2_exact_20261001/exact.py @@ -0,0 +1,432 @@ +"""Exact derivation of c and c2 for the Gaussian kernel, d = 1, 2, 3 (NOTE.md sections 2-3). Standard library only. + +Polynomials (Poly) in (r, b, k, x1, x2) over Q, truncated at r^NR; closed-form numbers (Num) + sum q * pi^(h/2) * sqrt(m) * theta^e, q rational, h integer, m squarefree, e in {0, 1}, +where theta = M(0, 0) is the one cone atom of d = 3 (section 3.3). Each coefficient is returned as a Num times one +Gamma atom: c = Num * Gamma(1/6) 12^(-1/6), c2 = Num * Gamma(5/6) 12^(-5/6).""" +from fractions import Fraction as Fr +import itertools +import math + +import pseries as P + +NR = 6 +VARS = ('r', 'b', 'k', 'x1', 'x2') +NV = len(VARS) +ONE = (0,) * NV +MUTANT = None + + +# ------------------------------------------------------------------------------------------- multivariate polynomials +class Poly: + __slots__ = ('t',) + + def __init__(self, t=None): + self.t = {e: c for e, c in (t or {}).items() if c} + + @staticmethod + def var(name, power=1): + e = [0] * NV + e[VARS.index(name)] = power + return Poly({tuple(e): Fr(1)}) + + @staticmethod + def const(c): + return Poly({ONE: Fr(c)}) + + def __add__(self, o): + if not isinstance(o, Poly): + o = Poly.const(o) + t = dict(self.t) + for e, c in o.t.items(): + t[e] = t.get(e, 0) + c + return Poly(t) + __radd__ = __add__ + + def __neg__(self): + return Poly({e: -c for e, c in self.t.items()}) + + def __sub__(self, o): + return self + (-o if isinstance(o, Poly) else Poly.const(-o)) + + def __mul__(self, o): + if not isinstance(o, Poly): + o = Fr(o) + return Poly({e: c * o for e, c in self.t.items()}) + t = {} + for e1, c1 in self.t.items(): + for e2, c2 in o.t.items(): + if e1[0] + e2[0] > NR: + continue + e = tuple(a + b for a, b in zip(e1, e2)) + t[e] = t.get(e, 0) + c1 * c2 + return Poly(t) + __rmul__ = __mul__ + + def rcoef(self, n): + return Poly({(0,) + e[1:]: c for e, c in self.t.items() if e[0] == n}) + + def degree(self, name): + i = VARS.index(name) + return max((e[i] for e in self.t), default=0) + + def coef_of(self, name, p): + i = VARS.index(name) + return Poly({e[:i] + (0,) + e[i + 1:]: c for e, c in self.t.items() if e[i] == p}) + + def scalar(self): + if set(self.t) - {ONE}: + raise ValueError('not a constant: %r' % (self.t,)) + return self.t.get(ONE, Fr(0)) + + def subs(self, name, poly): + i = VARS.index(name) + out, powers = Poly(), [Poly.const(1)] + for e, c in self.t.items(): + while len(powers) <= e[i]: + powers.append(powers[-1] * poly) + out = out + Poly({e[:i] + (0,) + e[i + 1:]: c}) * powers[e[i]] + return out + + def is_zero(self): + return not self.t + + +def from_series(s): + return Poly({(i,) + (0,) * (NV - 1): c for i, c in enumerate(s.c[:NR + 1]) if c}) + + +# ------------------------------------------------------------------------------------------- closed-form numbers +def squarefree(n): + q, m, p = 1, 1, 2 + while p * p <= n: + while n % (p * p) == 0: + n //= p * p + q *= p + if n % p == 0: + n //= p + m *= p + p += 1 + return q, m * n + + +def sqrt_frac(x): + x = Fr(x) + if x <= 0: + raise ValueError('sqrt of a non-positive rational') + q, m = squarefree(x.numerator * x.denominator) + return Fr(q, x.denominator), m + + +class Num: + __slots__ = ('t',) + + def __init__(self, t=None): + self.t = {k: v for k, v in (t or {}).items() if v} + + @staticmethod + def rat(q): + return Num({(0, 1, 0): Fr(q)}) + + @staticmethod + def sqrt(x): + q, m = sqrt_frac(x) + return Num({(0, m, 0): q}) + + @staticmethod + def pi_half(h): + return Num({(h, 1, 0): Fr(1)}) + + def __add__(self, o): + t = dict(self.t) + for k, v in o.t.items(): + t[k] = t.get(k, 0) + v + return Num(t) + + def __neg__(self): + return Num({k: -v for k, v in self.t.items()}) + + def __sub__(self, o): + return self + (-o) + + def __eq__(self, o): + return isinstance(o, Num) and not (self - o).t + + def __mul__(self, o): + if not isinstance(o, Num): + return Num({k: v * Fr(o) for k, v in self.t.items()}) + t = {} + for (h1, m1, e1), v1 in self.t.items(): + for (h2, m2, e2), v2 in o.t.items(): + g = math.gcd(m1, m2) + key = (h1 + h2, (m1 // g) * (m2 // g), e1 + e2) + t[key] = t.get(key, 0) + v1 * v2 * g + return Num(t) + __rmul__ = __mul__ + + def text(self): + parts = [] + for (h, m, e), v in sorted(self.t.items()): + s = str(v) + if m != 1: + s += ' sqrt(%d)' % m + if h: + s += ' pi^(%s)' % Fr(h, 2) + if e: + s += ' theta' + ('' if e == 1 else '^%d' % e) + parts.append(s) + return ' + '.join(parts) if parts else '0' + + +# ------------------------------------------------------------------------------------------- the kernel A_r +def forms_for(d): + forms = {} + for nm, sg in (('M', -1), ('S', 1)): + for i in range(d): + for j in range(i, d): + be = [0] * d + be[i] += 1 + be[j] += 1 + forms[(nm, i, j)] = P.deriv_at(tuple(be), sg, d) + trans = [] + for i in range(1, d): + for j in range(i, d): + be = [0] * d + be[i] += 1 + be[j] += 1 + forms[('T', i, j)] = {tuple(be): P.S.const(1)} + trans.append(('T', i, j)) + return forms, trans + + +def matchings(lst): + if not lst: + yield [] + return + f, rest = lst[0], lst[1:] + for mm in matchings(rest): + yield mm + for i, y in enumerate(rest): + for mm in matchings(rest[:i] + rest[i + 1:]): + yield [(f, y)] + mm + + +def perm_sign(p): + return -1 if sum(1 for i in range(len(p)) for j in range(i + 1, len(p)) if p[i] > p[j]) % 2 else 1 + + +def series_pow(x, alpha): + """(1 + x)^alpha for a Poly x without r^0 part""" + out, term, coef = Poly.const(1), Poly.const(1), Fr(1) + for n in range(1, NR + 1): + coef = coef * (Fr(alpha) - (n - 1)) / n + term = term * x + out = out + term * coef + return out + + +def series_exp(x): + out, term = Poly.const(1), Poly.const(1) + for n in range(1, NR + 1): + term = term * x * Fr(1, n) + out = out + term + return out + + +def kernel_data(d): + """A~_r(b, k) = const * [cone integral of exp(-(q0 + Q0)/2) G w] / r^2 (NOTE.md 2.3; the fixed-cone surrogate), with + w = 1 (d = 1), da over a < 0 (d = 2), (x1 - x2) dx over x2 < x1 < 0 (d = 3; the factor pi is in const).""" + forms, trans = forms_for(d) + pin = P.Pinned(d, forms) + p, nt = pin.p, len(trans) + r, b, k = Poly.var('r'), Poly.var('b'), Poly.var('k') + tgt = 11 if MUTANT == 'target' else 12 + v = [b - k * r * r * r * Fr(1, 2), -(k * r * r), Poly(), k * tgt] + [Poly()] * (p - 4) + mean = {u: sum((from_series(pin.G[u][j]) * v[j] for j in range(p)), Poly()) for u in pin.names} + hn = [u for u in pin.names if u not in trans] + if nt: + CTT = [[pin.C[(x, y)] for y in trans] for x in trans] + CTTinv = P.mat_inv(CTT) + detT = P.mat_det(CTT) + B = {u: [sum((pin.C[(u, trans[j])] * CTTinv[j][i] for j in range(nt)), P.S()) for i in range(nt)] for u in hn} + cp = {(u, w): from_series(pin.C[(u, w)] - sum((B[u][i] * pin.C[(trans[i], w)] for i in range(nt)), P.S())) + for u in hn for w in hn} + tval = {('T', 1, 1): Poly.var('x1'), ('T', 1, 2): Poly(), ('T', 2, 2): Poly.var('x2')} + zT = [tval[nm] - mean[nm] for nm in trans] + aff = {u: mean[u] + sum((from_series(B[u][i]) * zT[i] for i in range(nt)), Poly()) for u in hn} + else: + CTTinv, detT, zT = [], P.S.const(1), [] + cp = {(u, w): from_series(pin.C[(u, w)]) for u in hn for w in hn} + aff = dict((u, mean[u]) for u in hn) + F = Poly() # E[det H_M det H_S | T, U_r = v_r] by Isserlis + for pp in itertools.permutations(range(d)): + for qq in itertools.permutations(range(d)): + fac = [('M', min(i, pp[i]), max(i, pp[i])) for i in range(d)] + \ + [('S', min(i, qq[i]), max(i, qq[i])) for i in range(d)] + sg = perm_sign(pp) * perm_sign(qq) + for mt in matchings(list(range(2 * d))): + used = set(i for pr in mt for i in pr) + c = Poly.const(sg) + for (i, j) in mt: + c = c * cp[(fac[i], fac[j])] + for i in range(2 * d): + if i not in used: + c = c * aff[fac[i]] + F = F + c + q = sum((v[i] * from_series(pin.Sinv[i][j]) * v[j] for i in range(p) for j in range(p)), Poly()) + Q = sum((zT[i] * from_series(CTTinv[i][j]) * zT[j] for i in range(nt) for j in range(nt)), Poly()) + q0, Q0 = q.rcoef(0), Q.rcoef(0) + detS0, detT0 = pin.detS.c[0], detT.c[0] + ratio = from_series(pin.detS) * (1 / detS0) * from_series(detT) * (1 / detT0) - Poly.const(1) + G = series_pow(ratio, Fr(-1, 2)) * series_exp(((q - q0) + (Q - Q0)) * Fr(-1, 2)) * F + nn = p + nt + const = Num.rat(-12) * Num.sqrt(Fr(1, 2 ** nn)) * Num.pi_half(-nn) * Num.sqrt(Fr(1) / (detS0 * detT0)) + if d == 3: + const = const * Num.pi_half(2) + return {'G': G, 'q0': q0, 'Q0': Q0, 'const': const, 'p': p, 'detS0': detS0, 'detT0': detT0, 'F': F} + + +# ------------------------------------------------------------------------------------------- the integrations +SPHERE = {1: Num.rat(2), 2: Num.pi_half(2) * 2, 3: Num.pi_half(2) * 4} + + +def half_moment(j, c): + """int_0^oo t^j e^(-c t^2) dt = Gamma((j+1)/2) c^(-(j+1)/2)/2 (c > 0 rational)""" + if j % 2 == 0: + g = Fr(1) + for i in range(j // 2): + g *= Fr(2 * i + 1, 2) + return Num.pi_half(1) * Num.sqrt(Fr(1) / c) * (g / c ** (j // 2) / 2) + m = (j + 1) // 2 + return Num.rat(Fr(math.factorial(m - 1)) / c ** m / 2) + + +def quadrant_moments(al, be, ga, nmax): + """M(i, j) = int int_(s, t > 0) s^i t^j exp(-(al s^2 + be s t + ga t^2)) ds dt, i + j <= nmax, from + (i) 2 al M(i+1, j) + be M(i, j+1) = i M(i-1, j) + [i = 0] H_j(ga), + (ii) be M(i+1, j) + 2 ga M(i, j+1) = j M(i, j-1) + [j = 0] H_i(al) (H the half-line moments), + integration by parts in s and in t. M(0, 0) = theta = (pi/2 - arctan(be/sqrt(D)))/sqrt(D), D = 4 al ga - be^2. + Every M(i+1, j) with j >= 1 is reached twice (from (i, j) and from (i+1, j-1)); the two values must agree.""" + D = 4 * al * ga - be * be + if not (D > 0 and al > 0 and ga > 0): + raise ValueError('cone exponent not positive definite') + M = {(0, 0): Num({(0, 1, 1): Fr(1)})} + for n in range(nmax): + for i in range(n + 1): + j = n - i + R1 = (M[(i - 1, j)] * i if i else half_moment(j, ga)) + R2 = (M[(i, j - 1)] * j if j else half_moment(i, al)) + vals = (((i + 1, j), (R1 * (2 * ga) - R2 * be) * (Fr(1) / D)), + ((i, j + 1), (R2 * (2 * al) - R1 * be) * (Fr(1) / D))) + for key, val in vals: + if key in M: + if not M[key] == val: + raise ArithmeticError('quadrant recurrence inconsistent at %r' % (key,)) + else: + M[key] = val + return M, D + + +def gamma_ratio(a, a0): + """Gamma(a)/Gamma(a0) for a - a0 an integer (exact)""" + out, x = Fr(1), Fr(a0) + while x < a: + out *= x + x += 1 + while x > a: + x -= 1 + out /= x + return out + + +def coefficient(d, which, data=None): + """which = 0: c = (|S^(d-1)|/3) int db int_0^oo A_0 k^(-2/3) dk; + which = 2: c2 = (|S^(d-1)|/3) int db f.p. int_0^oo A_2 k^(-4/3) dk (the finite part = the subtraction of A_2(b, 0)). + Returns (Num, Gamma atom, info).""" + D_ = data or kernel_data(d) + G = D_['G'] + for n in (0, 1): + if not G.rcoef(n).is_zero(): + raise ArithmeticError('integrand does not vanish to order r^2') + if not G.rcoef(3).is_zero(): + raise ArithmeticError('A_r has an r^1 term') + H = G.rcoef(which + 2) + E = (D_['q0'] + D_['Q0']) * Fr(1, 2) + # --- k: E = ak k^2 + (k-free part) + ak = E.coef_of('k', 2).scalar() + if not (E.coef_of('k', 1).is_zero() and E.degree('k') == 2 and ak == 12): + raise ArithmeticError('unexpected k-dependence of the exponent') + a0 = Fr(1, 6) if which == 0 else Fr(5, 6) + expo = Fr(-2, 3) if which == 0 else Fr(-4, 3) + Hk = Poly() + for j in range(H.degree('k') + 1): + cj = H.coef_of('k', j) + if cj.is_zero(): + continue + if j % 2: + raise ArithmeticError('odd power of k') + # int_0^oo k^(j + expo) e^(-12 k^2) dk (finite part) = Gamma(a) 12^(-a)/2, a = (j + expo + 1)/2 + a = (j + expo + 1) / 2 + if MUTANT == 'finite-part' and j == 0 and which == 2: + a = a0 # drops the subtraction of A_2(b, 0): Gamma(-1/6) replaced by Gamma(5/6) + shift = a - a0 + Hk = Hk + cj * (gamma_ratio(a, a0) * Fr(12) ** (-int(shift)) / 2) + katom = 'Gamma(1/6) 12^(-1/6)' if which == 0 else 'Gamma(5/6) 12^(-5/6)' + # --- b: k-free exponent = ab b^2 + b L(x) + R(x); complete the square + Ek = E.coef_of('k', 0) + ab = Ek.coef_of('b', 2).scalar() + if Ek.degree('b') != 2: + raise ArithmeticError('exponent not quadratic in b') + Lx, Rx = Ek.coef_of('b', 1), Ek.coef_of('b', 0) + shifted = Hk.subs('b', Poly.var('b') - Lx * (Fr(1) / (2 * ab))) + Hb = {} + for e, c in shifted.t.items(): + if e[1] % 2: + continue + key = e[3:] + Hb[key] = Hb.get(key, Num()) + half_moment(e[1], ab) * (2 * c) + if not Hb: + raise ArithmeticError('empty integrand: the r-truncation is too low') + Ex = Rx - Lx * Lx * (Fr(1) / (4 * ab)) + total, info = Num(), {} + if d == 1: + if not Ex.is_zero(): + raise ArithmeticError('d = 1: residual exponent') + for val in Hb.values(): + total = total + val + elif d == 2: + g = Ex.coef_of('x1', 2).scalar() + if not (Ex.degree('x1') == 2 and Ex.coef_of('x1', 1).is_zero() and Ex.coef_of('x1', 0).is_zero()): + raise ArithmeticError('d = 2: cone exponent') + for key, val in Hb.items(): + total = total + val * half_moment(key[0], g) * (-1) ** key[0] # int over a < 0 + info['gamma_a'] = g + else: + X1, X2 = -Poly.var('x1'), -Poly.var('x1') - Poly.var('x2') # x1 = -s, x2 = -s - t (s -> x1, t -> x2 slots) + + def subst(poly): + out = Poly() + for e, c in poly.t.items(): + term = Poly({e[:3] + (0, 0): c}) + for _ in range(e[3]): + term = term * X1 + for _ in range(e[4]): + term = term * X2 + out = out + term + return out + Est = subst(Ex) + al = Est.coef_of('x1', 2).coef_of('x2', 0).scalar() + ga = Est.coef_of('x2', 2).coef_of('x1', 0).scalar() + be = Est.coef_of('x1', 1).coef_of('x2', 1).scalar() + chk = Poly.var('x1', 2) * al + Poly.var('x1') * Poly.var('x2') * be + Poly.var('x2', 2) * ga + if not (Est - chk).is_zero(): + raise ArithmeticError('d = 3: cone exponent not a quadratic form') + nmax = max(key[0] + key[1] for key in Hb) + 2 + M, Dc = quadrant_moments(al, -be if MUTANT == 'cone-sign' else be, ga, nmax) + for key, val in Hb.items(): + mono = subst(Poly({(0, 0, 0) + key: Fr(1)})) * Poly.var('x2') # weight x1 - x2 = t + for e, c in mono.t.items(): + total = total + val * M[(e[3], e[4])] * c + info.update({'alpha': al, 'beta': be, 'gamma': ga, 'D': Dc}) + value = SPHERE[d] * Fr(1, 3) * D_['const'] * total + return value, katom, info diff --git a/frontiers/c2_exact_20261001/gam.py b/frontiers/c2_exact_20261001/gam.py new file mode 100644 index 000000000..e7707bf14 --- /dev/null +++ b/frontiers/c2_exact_20261001/gam.py @@ -0,0 +1,44 @@ +"""Gamma at rationals and rational powers, as outward-rounded intervals. These are the functions lngamma, gamma and rpow +of frontiers/cusp_coefficient_certified_20261001/consts.py (Math-#219); their code lines are unchanged. + +log Gamma(y) for y = q + N >= 10 uses Stirling's series with the remainder bounded by the first omitted term (real y > 0), +and log Gamma(q) = log Gamma(q + N) - log prod_(j= 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/c2_exact_20261001/pseries.py b/frontiers/c2_exact_20261001/pseries.py new file mode 100644 index 000000000..fac3c562b --- /dev/null +++ b/frontiers/c2_exact_20261001/pseries.py @@ -0,0 +1,203 @@ +"""Exact truncated power series in the separation r, and the pinned two-point structure (NOTE.md section 2). + +Field: centred Gaussian on R^d with covariance K(z) = exp(-|z|^2/2), so Cov(d^a f(0), d^b f(0)) = (-1)^|b| d^(a+b) K(0), +d^g K(0) = prod_i (-1)^(g_i/2) (g_i - 1)!! for even g and 0 otherwise. Points M = -r e1/2, S = +r e1/2. + +Every quantity is a linear form in the jets at 0 whose coefficients are power series in r: + d^beta f(x e1) = sum_j x^j/j! d^(beta + j e1) f(0), x = -+ r/2. +A Series keeps the coefficients of r^0..r^R exactly (Fractions). The pin rows divide by r, r and r^3; each division is +exact (the low coefficients vanish identically, which `shift` asserts), and costs at most three orders, so every derived +series is exact through r^(R - 3) = r^7. The certificate needs r^6 (exact.NR).""" +from fractions import Fraction as Fr + +R = 10 + + +class S: + """a power series c_0 + c_1 r + ... + c_R r^R over Q (higher orders dropped)""" + __slots__ = ('c',) + + def __init__(self, c=None): + c = list(c or []) + self.c = (c + [Fr(0)] * (R + 1))[:R + 1] + + @staticmethod + def const(x): + return S([Fr(x)]) + + def __add__(self, o): + if not isinstance(o, S): + o = S.const(o) + return S([a + b for a, b in zip(self.c, o.c)]) + __radd__ = __add__ + + def __neg__(self): + return S([-a for a in self.c]) + + def __sub__(self, o): + return self + (-o if isinstance(o, S) else S.const(-o)) + + def __rsub__(self, o): + return (-self) + o + + def __mul__(self, o): + if not isinstance(o, S): + return S([a * o for a in self.c]) + out = [Fr(0)] * (R + 1) + for i, a in enumerate(self.c): + if a: + for j in range(R + 1 - i): + if o.c[j]: + out[i + j] += a * o.c[j] + return S(out) + __rmul__ = __mul__ + + def inv(self): + a0 = self.c[0] + if a0 == 0: + raise ZeroDivisionError('series with zero constant term') + out = [Fr(0)] * (R + 1) + out[0] = 1 / a0 + for n in range(1, R + 1): + out[n] = -sum(self.c[k] * out[n - k] for k in range(1, n + 1)) / a0 + return S(out) + + def is_zero(self): + return not any(self.c) + + +def dfact(n): + out = 1 + while n > 1: + out *= n + n -= 2 + return out + + +def rho_der(g): + """d^g K(0) for K = exp(-|z|^2/2)""" + if any(x % 2 for x in g): + return 0 + out = 1 + for x in g: + out *= (-1) ** (x // 2) * dfact(x - 1) + return out + + +def cov(u, w): + """Cov of two linear forms (dict jet -> S)""" + tot = S() + for a, sa in u.items(): + for b, sb in w.items(): + v = rho_der(tuple(x + y for x, y in zip(a, b))) + if v: + tot = tot + sa * sb * ((-1) ** sum(b) * v) + return tot + + +def deriv_at(beta, sign, d): + """d^beta f at x = sign r/2 along e1, as jets: sum_(j <= R) (sign/2)^j r^j/j! d^(beta + j e1) f(0)""" + out = {} + fact = 1 + for j in range(R + 1): + if j: + fact *= j + c = [Fr(0)] * (R + 1) + c[j] = Fr(sign, 2) ** j / fact + out[(beta[0] + j,) + tuple(beta[1:])] = S(c) + return out + + +def lin(*terms): + """sum of coef * form""" + out = {} + for coef, form in terms: + for g, s in form.items(): + out[g] = out.get(g, S()) + s * coef + return {g: s for g, s in out.items() if not s.is_zero()} + + +def shift(form, p): + """divide by r^p; the coefficients of r^0..r^(p-1) must vanish identically""" + out = {} + for g, s in form.items(): + if any(x != 0 for x in s.c[:p]): + raise ArithmeticError('inexact division by r^%d at jet %r' % (p, g)) + out[g] = S(s.c[p:]) + return out + + +def pin_rows(d): + """The symmetric pin vector U_r of [R] section 2 (Math-#216 c2_check.py Pinned.pin_rows): + (f(M) + f(S))/2, (f(S) - f(M))/r, (f_x(S) - f_x(M))/r, 6/r^2 (f_x(M) + f_x(S) - 2 (f(S) - f(M))/r), + and for each transverse direction y: (f_y(M) + f_y(S))/2, (f_y(S) - f_y(M))/r. + Target v_r = (b - k r^3/2, -k r^2, 0, 12 k, 0, ..., 0).""" + e0 = (0,) * d + ex = (1,) + (0,) * (d - 1) + fa, fc = deriv_at(e0, -1, d), deriv_at(e0, 1, d) + fxa, fxc = deriv_at(ex, -1, d), deriv_at(ex, 1, d) + rows = [lin((Fr(1, 2), fa), (Fr(1, 2), fc)), + shift(lin((1, fc), (-1, fa)), 1), + shift(lin((1, fxc), (-1, fxa)), 1)] + num = lin((1, fxa), (1, fxc), (-2, shift(lin((1, fc), (-1, fa)), 1))) + rows.append({g: s * 6 for g, s in shift(num, 2).items()}) + for t in range(1, d): + ey = tuple(1 if i == t else 0 for i in range(d)) + ya, yc = deriv_at(ey, -1, d), deriv_at(ey, 1, d) + rows += [lin((Fr(1, 2), ya), (Fr(1, 2), yc)), shift(lin((1, yc), (-1, ya)), 1)] + return rows + + +def mat_inv(M): + """inverse of a square matrix of series whose r^0 part is invertible (Gauss-Jordan over Q[[r]])""" + n = len(M) + A = [row[:] + [S.const(int(i == j)) for j in range(n)] for i, row in enumerate(M)] + for c in range(n): + p = next(i for i in range(c, n) if A[i][c].c[0] != 0) + A[c], A[p] = A[p], A[c] + iv = A[c][c].inv() + A[c] = [x * iv for x in A[c]] + for i in range(n): + if i != c and not A[i][c].is_zero(): + f = A[i][c] + A[i] = [x - f * y for x, y in zip(A[i], A[c])] + return [row[n:] for row in A] + + +def mat_det(M): + n = len(M) + A = [row[:] for row in M] + det = S.const(1) + for c in range(n): + p = next(i for i in range(c, n) if A[i][c].c[0] != 0) + if p != c: + A[c], A[p] = A[p], A[c] + det = -det + det = det * A[c][c] + iv = A[c][c].inv() + for i in range(c + 1, n): + if not A[i][c].is_zero(): + f = A[i][c] * iv + A[i] = [x - f * y for x, y in zip(A[i], A[c])] + return det + + +class Pinned: + """conditional law of linear forms given U_r: mean G v (G[u][j] series), covariance C[(u, w)] (series)""" + + def __init__(self, d, forms): + self.d = d + self.rows = pin_rows(d) + p = len(self.rows) + self.p = p + Sm = [[cov(a, b) for b in self.rows] for a in self.rows] + self.Sinv = mat_inv(Sm) + self.detS = mat_det(Sm) + self.names = list(forms) + CF = {u: [cov(forms[u], row) for row in self.rows] for u in self.names} + self.G = {u: [sum((CF[u][q] * self.Sinv[q][j] for q in range(p)), S()) for j in range(p)] for u in self.names} + self.C = {} + for i, u in enumerate(self.names): + for w in self.names[i:]: + val = cov(forms[u], forms[w]) - sum((self.G[u][q] * CF[w][q] for q in range(p)), S()) + self.C[(u, w)] = self.C[(w, u)] = val