Wire critical symbolic anchors into SageMath CI (C1/C2/M2 + Bach on non-Einstein metric)

CHANGELOG.md

@@ -1,5 +1,22 @@
 # Changelog
 
+## [Unreleased]
+
+### Added
+- `tests/test_pe/test_critical_anchors.py`: SageMath regression anchors that
+  re-prove the audit's critical fixes on discriminating geometries — round S³
+  (Fefferman–Graham `g₄=1/16 g₀` (C1), renormalized volume `v₂=−3/4` (C2),
+  `Q₄=15/8`) and the non-Einstein product S²(1)×S²(2) (Schouten 7/24, −1/3;
+  `J=5/12`; conformal Laplacian `P₂(1)=−5/12` isolating the curvature term (M2);
+  `Bach≠0`). These run in the existing Track A SageMath CI. (167 tests total.)
+
+### Fixed
+- Symbolic curvature operators (`P₂`, Paneitz `P₄`, `Q₄`) formed rational
+  coefficients via Python float division on `cs.dimension` (a Python int),
+  contaminating results with floats (masked when the exact value is
+  binary-representable). Reordered so a Sage object is divided by the integer
+  denominator, keeping outputs exact (e.g. `P₂(1)=−5/12`, not `−0.41666…`).
+
 ## [0.1.1] - 2026-05-31
 
 ### Fixed (mathematical correctness audit)

README.md

@@ -246,7 +246,7 @@ conformal-toolkit/
 │   ├── discrete/               # Curvature, Q, Bach, Willmore, cross-ratios, Yamabe, spectral
 │   ├── features/               # mesh_conformal_features() pipeline
 │   └── benchmarks/             # ShapeNet, SHREC, FAUST evaluation (WIP)
-├── tests/                      # 160 tests across both packages
+├── tests/                      # 167 tests across both packages
 ├── examples/                   # 6 Jupyter notebooks
 └── paper.md                    # JOSS paper draft
 ```
@@ -367,7 +367,13 @@ Every push and pull request runs **both** automatically on GitHub Actions
 ([`.github/workflows/test.yml`](.github/workflows/test.yml)) — so the symbolic
 formulas in this README are re-verified by a real SageMath install in the cloud,
 not just asserted. (You can watch it: the green check on a commit means Track A
-recomputed things like `Q₄(S⁴)=6` from scratch.)
+recomputed things like `Q₄(S⁴)=6` from scratch.) In particular,
+`tests/test_pe/test_critical_anchors.py` re-proves the audit's *critical* fixes
+on geometries chosen to expose them — the round S³ (where the corrected
+Fefferman–Graham `g₄=1/16 g₀` and renormalized volume `v₂=−J/2` differ from the
+old buggy formulas, which only agreed at n=4) and the **non-Einstein** product
+S²(1)×S²(2) (where Bach≠0 and the conformal Laplacian's curvature term is
+nonzero, unlike the Ricci-flat/Einstein metrics that would mask those bugs).
 
 The part worth stealing if you're learning SageMath: **how to get Sage into CI.**
 Sage has no usable pip wheel, but it *is* on conda-forge, so the trick is to

conformal_toolkit/core/gjms.py

@@ -26,7 +26,12 @@ def laplacian_operator(cs, f):
     delta_f = laplacian(cs.connection(), cs.metric, f)
     if n == 2:
         return delta_f
-    return delta_f - ((n - 2) / (4 * (n - 1))) * R * f
+    # Order the arithmetic so the rational coefficient stays EXACT: dividing a
+    # Sage object (R * f) by the integer denominator avoids Python float
+    # division of (n-2)/(4(n-1)), which would contaminate the symbolic result
+    # (cs.dimension is a Python int). Cf. schouten.py, which derives n from
+    # g.domain().dim() (a Sage Integer) and is exact for the same reason.
+    return delta_f - (n - 2) * R * f / (4 * (n - 1))
 
 
 def paneitz_operator(cs, f):
@@ -47,8 +52,9 @@ def paneitz_operator(cs, f):
     term2 = divergence(nabla, g, V_df)
 
     Q4 = cs.q_curvature(order=4)
-    coeff = (n - 4) / 2
-    term3 = coeff * Q4 * f
+    # Exact rational arithmetic (see laplacian_operator): divide the Sage
+    # object by the integer 2 rather than forming the Python float (n-4)/2.
+    term3 = (n - 4) * Q4 * f / 2
 
     return delta2_f + term2 + term3
 

conformal_toolkit/core/q_curvature.py

@@ -42,7 +42,9 @@ def _q4(cs):
     # Then trace the remaining contra/covariant pair
     P_norm_sq = P.contract(0, P_up, 0).trace(0, 1)
 
-    n_half = n / 2
-    Q4 = -delta_J - 2 * P_norm_sq + n_half * J * J
+    # Exact rational arithmetic: form n*J*J/2 (Sage object / integer) rather
+    # than the Python float n/2, which would contaminate the symbolic result
+    # when cs.dimension is a Python int (e.g. (n/2) for odd n).
+    Q4 = -delta_J - 2 * P_norm_sq + n * J * J / 2
 
     return Q4

docs/TOOLING_GAPS.md

@@ -77,18 +77,28 @@ hypersurface count (`M5`) could only be bounded below (`> 4`), not pinned.
 returning an independent generating set (algebraic + derivative invariants).
 Hard, but it is what makes `count_invariants` checkable exactly.
 
-## Gap 5 — No library-execution harness in CI *(P1, highest ROI after backends)*
-
-We could not run the repo's own Sage code, so prose-vs-return mismatches went
-unverified (does `q_curvature(order=4)` literally return 6? does
-`count_invariants` return the stated ints?).
-
-**Proposal — a Sage-enabled CI job** that evaluates each operator on canonical
-anchors — round `Sⁿ`, `S²×S²`, **and a non-conformally-flat metric** — and
-asserts known values (Branson `Q_n=(n−1)!`, conformal invariance / vanishing).
-This catches `M2`, `M8`, `C1`/`C2`, `M16` **empirically, with no symbolic engine
-at all.** The anchor library is already specified in `ERRATA.md → How we caught
-them`; turning it into `tests/test_anchors/` is the concrete first deliverable.
+## Gap 5 — Library-execution harness in CI *(P1 — ✅ LARGELY ADDRESSED)*
+
+Original problem: we could not run the repo's own Sage code, so prose-vs-return
+mismatches went unverified, and the *critical* fixes were witnessed only by
+hand-derivation.
+
+**Done.** The repo's GitHub Actions already provisions SageMath (Track A) and
+runs the symbolic suite on every push. `tests/test_pe/test_critical_anchors.py`
+now turns the audit anchors into automated regression tests on discriminating
+geometries:
+- round **S³** (n≠4, where the bug's wrong n-dependence no longer coincides with
+  the correct value): asserts `v₂ = −3/4` (C2), Fefferman–Graham `g₄ = 1/16 g₀`
+  (C1), `Q₄ = 15/8`;
+- non-Einstein **S²(1)×S²(2)**: asserts Schouten `P[0,0]=7/24, P[2,2]=−1/3`,
+  `J=5/12`, the conformal Laplacian `P₂(1)=−5/12` (isolating the curvature term,
+  M2), and `Bach≠0` — the regime an Einstein/Ricci-flat metric would mask.
+
+So C1, C2, and M2 are now re-proven by CI every push. **Remaining:** add a
+higher even-dimensional sphere (S⁶) cross-check, and anchors for the
+research-level items (`M8` extrinsic Q₄ once implemented, `C1`'s Bach
+differential terms). The discriminating-anchor pattern and the exact closed
+forms are documented in `ERRATA.md → How we caught them`.
 
 ## Gap 6 — Literature retrieval isn't citation-grade *(P3)*
 
@@ -108,13 +118,13 @@ rendered PDF) and indexes numbered equations for verbatim citation lookup.
 |-----|----------|--------|---------|
 | Backend healthcheck + Wolfram key | **P0** | Low | the silent-degradation risk itself |
 | Tensor / abstract-index verifier | **P0** | High | `M2`, `M3`, `M8`, `M12`, `C1`, `M11` |
-| Library-execution CI harness | **P1** | Low–Med | `M2`, `M8`, `C1`, `C2`, `M16` (empirically) |
+| Library-execution CI harness | ✅ **done** | Low–Med | `M2`, `C1`, `C2` re-proven in CI (S³ + non-Einstein S²×S²) |
 | Conformal-weight checker | **P1** | Low | `M4`, `M9`, `m4` |
 | Degenerate-metric / Carroll engine | **P2** | Med | `M9`, `M10`, `m5` |
 | Invariant-basis enumerator | **P2** | High | `M5`, `M6` |
 | arXiv source fetcher | **P3** | Low | `M5`, `M7`, `m3` |
 
-**Do first:** the backend healthcheck (P0) and the library-execution CI harness
-(P1). Together they are low-effort and would have caught the two *critical*
+**Do first:** the backend healthcheck (P0). The library-execution CI harness
+(P1) is now ✅ in place for the critical fixes (C1/C2/M2). Together they are low-effort and would have caught the two *critical*
 errors and most majors automatically — the highest return per hour of the whole
 list.

paper.md

@@ -84,7 +84,7 @@ symbolic conformal geometry to discrete mesh representations usable in PyTorch.
 # Functionality
 
 The software consists of two Python packages, 40 source modules in total, with
-160 tests and 6 example notebooks.
+167 tests and 6 example notebooks.
 
 ## Symbolic Package: `conformal_toolkit`
 

tests/conftest_sage.py

@@ -59,6 +59,54 @@ def _make_flat_rn(n):
     return {'manifold': M, 'metric': g, 'chart': X, 'dim': n}
 
 
+def _make_round_sphere_3():
+    """Unit round S^3: g = dchi^2 + sin^2(chi) dtheta^2 + sin^2(chi) sin^2(theta) dphi^2.
+
+    A cheap n != 4 anchor (R = n(n-1) = 6, P = (1/2)g, J = n/2 = 3/2). Used to
+    discriminate the Fefferman-Graham g_4 (ERRATA C1) and renormalized-volume
+    v_2 (ERRATA C2) bugs, whose wrong n-dependence coincides with the correct
+    value ONLY at n = 4.
+    """
+    from sage.all import Manifold, sin
+    S = Manifold(3, 'S3', structure='Riemannian')
+    X = S.chart(r'chi:(0,pi):\chi theta:(0,pi):\theta phi:(0,2*pi):\phi')
+    chi, theta, phi = X[:]
+    g = S.metric('g')
+    g[0, 0] = 1
+    g[1, 1] = sin(chi)**2
+    g[2, 2] = sin(chi)**2 * sin(theta)**2
+    return {'manifold': S, 'metric': g, 'chart': X, 'coords': (chi, theta, phi), 'dim': 3}
+
+
+def _make_product_s2_s2():
+    """NON-EINSTEIN product S^2(1) x S^2(2): two round 2-spheres of radii 1 and 2.
+
+        g = [dth1^2 + sin^2(th1) dph1^2] + 4*[dth2^2 + sin^2(th2) dph2^2]
+
+    Factor Gaussian curvatures K1 = 1, K2 = 1/4, so Ric has eigenvalue ratios
+    (1, 1, 1/4, 1/4): NOT proportional to g, hence NON-EINSTEIN, and not
+    conformally flat. Verified exact values (independently checked in sympy):
+        R = 5/2,  P[0,0] = 7/24,  P[2,2] = -1/3,  J = trP = 5/12.
+    Unlike Schwarzschild (Ricci-flat) or any Einstein metric -- both of which
+    have Bach = 0 -- this metric has Bach != 0 (its Cotton tensor vanishes, so
+    Bach reduces to the algebraic P^{cd} W_{acbd} part, which is nonzero). It is
+    therefore the anchor that actually exercises the Bach computation and the
+    conformal Laplacian's curvature term in a regime where they don't vanish.
+    """
+    from sage.all import Manifold, sin
+    M = Manifold(4, 'S2xS2', structure='Riemannian')
+    X = M.chart(r'th1:(0,pi):\theta_1 ph1:(0,2*pi):\varphi_1 '
+                r'th2:(0,pi):\theta_2 ph2:(0,2*pi):\varphi_2')
+    th1, ph1, th2, ph2 = X[:]
+    g = M.metric('g')
+    g[0, 0] = 1
+    g[1, 1] = sin(th1)**2
+    g[2, 2] = 4
+    g[3, 3] = 4 * sin(th2)**2
+    return {'manifold': M, 'metric': g, 'chart': X,
+            'coords': (th1, ph1, th2, ph2), 'dim': 4}
+
+
 def _make_hyperbolic_2():
     """Hyperbolic plane H^2 in upper half-plane model: g = (dx^2 + dy^2)/y^2."""
     from sage.all import Manifold

tests/test_pe/test_critical_anchors.py

@@ -0,0 +1,153 @@
+"""Critical symbolic regression anchors for the May-2026 correctness audit.
+
+These pin the corrected closed-form values of the CRITICAL/major numeric fixes
+on geometries that actually DISCRIMINATE the bugs:
+
+  * round S^3 (n != 4): the Fefferman-Graham g_4 (ERRATA C1) and renormalized-
+    volume v_2 (ERRATA C2) bugs have the wrong n-dependence and coincide with
+    the correct value ONLY at n = 4, so a non-4 dimension is required to catch
+    them.
+
+  * non-Einstein S^2(1) x S^2(2): Einstein and Ricci-flat metrics (e.g.
+    Schwarzschild) have Bach = 0 and a Schouten that vanishes or is pure-trace,
+    which masks Bach/Weyl-coefficient and conformal-Laplacian-curvature bugs.
+    This product has Bach != 0 and a genuinely non-proportional Schouten, so it
+    exercises those terms (ERRATA M2 and the Bach code path).
+
+All expected values were derived in closed form and independently re-checked in
+sympy (see ERRATA.md "How we caught them"). Requires SageMath; skipped if
+SageManifolds is unavailable.
+"""
+import pytest
+
+pytest.importorskip("sage.all", reason="SageMath required for symbolic anchors")
+
+from sage.all import Rational
+from tests.conftest_sage import _make_round_sphere_3, _make_product_s2_s2
+
+
+def _frac(a, b):
+    return Rational((a, b))
+
+
+def _scalar(x):
+    """Simplified symbolic expression from a Sage scalar field/expr."""
+    if hasattr(x, "expr"):
+        x = x.expr()
+    return x.simplify_full()
+
+
+def _comp(tensor, frame, i, j):
+    c = tensor[frame, i, j]
+    if hasattr(c, "expr"):
+        c = c.expr()
+    return c.simplify_full()
+
+
+def _is_zero(expr):
+    """Robust symbolic zero test (subtract-then-simplify is safer than ==)."""
+    return bool(expr.simplify_full() == 0)
+
+
+# ----------------------------------------------------------------------------
+# Round S^3 (n != 4): discriminates C1 (g_4) and C2 (v_2).
+# ----------------------------------------------------------------------------
+
+def test_renormalized_volume_v2_s3_is_minus_three_quarters():
+    """ERRATA C2: v_2 = -J/2 (n-independent) = -3/4 on S^3.
+
+    The old -1/(n-2) J would give -3/2 here; the two agree only at n = 4.
+    """
+    g = _make_round_sphere_3()['metric']
+    from conformal_toolkit.poincare_einstein.renormalized_volume import (
+        renormalized_volume_coefficient,
+    )
+    v2 = _scalar(renormalized_volume_coefficient(g, order=2))
+    assert _is_zero(v2 - _frac(-3, 4)), f"v_2 on S^3 should be -3/4, got {v2}"
+
+
+def test_fg_g4_s3_is_one_sixteenth_g0():
+    """ERRATA C1: on the round sphere g_4 = (1/16) g_0 (algebraic piece).
+
+    The old spurious 1/(n-4) prefactor would give a different coefficient at
+    n != 4 (e.g. 7/80 at n = 6); here every diagonal component must be 1/16 g0.
+    """
+    data = _make_round_sphere_3()
+    g = data['metric']
+    frame = data['chart'].frame()
+    from conformal_toolkit.poincare_einstein.fefferman_graham import fg_coefficient_g4
+    g4 = fg_coefficient_g4(g)
+    for i in range(3):
+        gii = _comp(g, frame, i, i)
+        g4ii = _comp(g4, frame, i, i)
+        assert _is_zero(g4ii - gii / 16), \
+            f"g_4[{i},{i}] should be (1/16) g_0[{i},{i}] = {gii / 16}, got {g4ii}"
+
+
+def test_q4_s3_is_fifteen_eighths():
+    """Branson Q_4 (order 4) on S^3 = -n/2 + n^3/8 = 15/8 at n = 3."""
+    data = _make_round_sphere_3()
+    from conformal_toolkit.core.conformal_structure import ConformalStructure
+    cs = ConformalStructure(data['metric'])
+    q4 = _scalar(cs.q_curvature(order=4))
+    assert _is_zero(q4 - _frac(15, 8)), f"Q_4 on S^3 should be 15/8, got {q4}"
+
+
+# ----------------------------------------------------------------------------
+# Non-Einstein S^2(1) x S^2(2): discriminates M2 and exercises Bach != 0.
+# ----------------------------------------------------------------------------
+
+def test_schouten_product_components():
+    """Schouten on the non-Einstein product: P[0,0]=7/24, P[2,2]=-1/3."""
+    data = _make_product_s2_s2()
+    frame = data['chart'].frame()
+    from conformal_toolkit.core.schouten import compute_schouten
+    P = compute_schouten(data['metric'])
+    p00 = _comp(P, frame, 0, 0)
+    p22 = _comp(P, frame, 2, 2)
+    assert _is_zero(p00 - _frac(7, 24)), f"P[0,0] should be 7/24, got {p00}"
+    assert _is_zero(p22 - _frac(-1, 3)), f"P[2,2] should be -1/3, got {p22}"
+
+
+def test_schouten_trace_product_is_five_twelfths():
+    """J = trace(P) = R/(2(n-1)) = 5/12 on S^2(1) x S^2(2)."""
+    g = _make_product_s2_s2()['metric']
+    from conformal_toolkit.core.schouten import schouten_trace
+    J = _scalar(schouten_trace(g))
+    assert _is_zero(J - _frac(5, 12)), f"J on S^2xS^2 should be 5/12, got {J}"
+
+
+def test_conformal_laplacian_has_curvature_term_product():
+    """ERRATA M2: the conformal Laplacian P_2 carries -(n-2)/(4(n-1)) R.
+
+    On a CONSTANT field f = 1 the bare Laplacian gives 0, so P_2(1) isolates the
+    curvature term: -(n-2)/(4(n-1)) R = -(2/12)(5/2) = -5/12 on this metric.
+    The pre-fix bare-Laplacian implementation returned 0 here.
+    """
+    data = _make_product_s2_s2()
+    M = data['manifold']
+    from conformal_toolkit.core.conformal_structure import ConformalStructure
+    cs = ConformalStructure(data['metric'])
+    f = M.scalar_field(1)
+    p2f = _scalar(cs.gjms_operator(f, order=2))
+    assert _is_zero(p2f - _frac(-5, 12)), \
+        f"P_2(1) should be -5/12 (the curvature term), got {p2f}"
+
+
+def test_bach_nonzero_on_non_einstein_product():
+    """The Bach tensor is NOT identically zero on the non-Einstein product.
+
+    Schwarzschild/Einstein metrics have Bach = 0 and would mask Bach-coefficient
+    bugs; here Cotton = 0 but the algebraic P^{cd} W part is nonzero, so a
+    correct Bach implementation must return a nonzero tensor.
+    """
+    data = _make_product_s2_s2()
+    frame = data['chart'].frame()
+    from conformal_toolkit.core.conformal_structure import ConformalStructure
+    cs = ConformalStructure(data['metric'])
+    B = cs.bach()
+    nonzero = any(
+        not _is_zero(_comp(B, frame, i, j))
+        for i in range(4) for j in range(4)
+    )
+    assert nonzero, "Bach should be nonzero on the non-Einstein S^2(1)xS^2(2)"