diff --git a/.github/sympy-parity/.gitignore b/.github/sympy-parity/.gitignore
new file mode 100644
index 000000000..6b58567d4
--- /dev/null
+++ b/.github/sympy-parity/.gitignore
@@ -0,0 +1,2 @@
+# Written by check.py on every run; the committed snapshot is sympy-baseline.json.
+sympy-current.json
diff --git a/.github/sympy-parity/check.py b/.github/sympy-parity/check.py
new file mode 100644
index 000000000..ab40d8787
--- /dev/null
+++ b/.github/sympy-parity/check.py
@@ -0,0 +1,136 @@
+# Re-runs the SymPy half of the parity survey against whatever SymPy pip installs today and
+# compares it with the snapshot in sympy-baseline.json.
+#
+# Only SymPy runs here. The survey's other half asks AngouriMath the same questions, and that
+# half is not in this repository, so nothing below compares the two libraries or restates a
+# verdict about either. What it answers is the narrower question that the published snapshot
+# cannot answer about itself: is SymPy still the version, and still giving the answers, that
+# the snapshot was taken from.
+#
+# Exit code 1 means the snapshot no longer describes SymPy, and a human has to re-take it.
+# https://github.com/asc-community/AngouriMath/issues/746 item 73.
+
+import json
+import os
+import signal
+import sys
+
+import sympy
+from sympy import *                                    # noqa: F401,F403  the probes are written against the flat namespace
+from sympy.abc import x, y, k, n, s, a, b, c           # noqa: F401
+
+# Every import below serves at least one probe. A release that moves or removes one of these is
+# itself a finding, so an import failure is recorded and left to surface as the probes that need
+# it failing, rather than taking the whole run down before anything is measured.
+FAILED_IMPORTS = []
+for statement in (
+    "import sympy.ntheory.modular",
+    "from sympy.stats import Die, P",
+    "from sympy.logic.boolalg import to_cnf",
+    "from sympy.logic.inference import satisfiable",
+    "from sympy.combinatorics import PermutationGroup",
+    "from sympy.vector import CoordSys3D",
+    "from sympy.tensor.tensor import TensorIndexType",
+    "from sympy.diffgeom import Manifold",
+    "from sympy.holonomic import expr_to_holonomic",
+    "from sympy.physics.units import convert_to",
+    "from sympy.printing.mathml import mathml",
+    "from sympy.discrete import convolution, fft",
+    "from sympy.simplify.sqrtdenest import sqrtdenest",
+    "from sympy.polys.numberfields import minimal_polynomial",
+    "from sympy.ntheory.continued_fraction import continued_fraction_iterator",
+    "from sympy.sets.conditionset import ConditionSet",
+    "from sympy.integrals.transforms import laplace_transform, fourier_transform",
+):
+    try:
+        exec(statement)
+    except Exception as e:                             # noqa: BLE001  any import failure is data
+        FAILED_IMPORTS.append(f"{statement} -- {type(e).__name__}: {e}")
+
+
+class Timeout(Exception):
+    pass
+
+
+def _alarm(sig, frm):
+    raise Timeout()
+
+
+signal.signal(signal.SIGALRM, _alarm)
+
+here = os.path.dirname(os.path.abspath(__file__))
+probes = json.load(open(os.path.join(here, "probes.json")))
+baseline = json.load(open(os.path.join(here, "sympy-baseline.json")))
+
+answers = []
+for p in probes:
+    signal.alarm(20)
+    try:
+        value = str(eval(p["sympy"]))
+    except Timeout:
+        value = "*** timeout (>20s) ***"
+    except Exception as e:                             # noqa: BLE001  an exception is the answer
+        value = f"!!! {type(e).__name__}: {e}".split("\n")[0]
+    finally:
+        signal.alarm(0)
+    answers.append({"area": p["area"], "what": p["what"], "sympy": value})
+
+current = {"sympy": sympy.__version__, "answers": answers}
+with open(os.path.join(here, "sympy-current.json"), "w") as f:
+    json.dump(current, f, indent=1)
+
+was = {(r["area"], r["what"]): r["sympy"] for r in baseline["answers"]}
+now = {(r["area"], r["what"]): r["sympy"] for r in answers}
+
+moved = [(k2, was[k2], now[k2]) for k2 in now if k2 in was and was[k2] != now[k2]]
+added = sorted(k2 for k2 in now if k2 not in was)
+gone = sorted(k2 for k2 in was if k2 not in now)
+version_moved = baseline["sympy"] != sympy.__version__
+
+
+def cell(t, at=70):
+    t = t.replace("|", "\\|").replace("\n", " ")
+    return t if len(t) <= at else t[: at - 1] + "…"
+
+
+lines = [
+    "## SymPy parity watch",
+    "",
+    f"- baseline: **SymPy {baseline['sympy']}**, taken {baseline['measured']}",
+    f"- installed today: **SymPy {sympy.__version__}**",
+    f"- {len(answers)} probes run, **{len(moved)}** answer{'' if len(moved) == 1 else 's'} moved",
+    "",
+]
+if FAILED_IMPORTS:
+    lines += ["Imports that no longer resolve:", ""] + [f"- `{i}`" for i in FAILED_IMPORTS] + [""]
+if moved:
+    lines += ["| area | capability | was | is now |", "|---|---|---|---|"]
+    lines += [
+        f"| {area} | {what} | `{cell(old)}` | `{cell(new)}` |" for (area, what), old, new in sorted(moved)
+    ]
+    lines.append("")
+if added or gone:
+    lines += [f"Probes added since the baseline: {added}", f"Probes gone since the baseline: {gone}", ""]
+
+summary = "\n".join(lines)
+print(summary)
+
+step_summary = os.environ.get("GITHUB_STEP_SUMMARY")
+if step_summary:
+    with open(step_summary, "a") as f:
+        f.write(summary + "\n")
+
+if not (version_moved or moved or added or gone):
+    print(f"The snapshot still describes SymPy {sympy.__version__}.")
+    sys.exit(0)
+
+# A version bump on its own is a failure because the published comparison names the version it
+# measured; a version it does not name is a version it did not measure, whether or not any of
+# these 80 answers happens to have moved with it.
+print(
+    "\nThe snapshot no longer describes SymPy. Re-run work/sympyparity end to end against "
+    f"SymPy {sympy.__version__}, republish the comparison, and replace sympy-baseline.json "
+    "with the sympy-current.json this run uploaded.",
+    file=sys.stderr,
+)
+sys.exit(1)
diff --git a/.github/sympy-parity/probes.json b/.github/sympy-parity/probes.json
new file mode 100644
index 000000000..dda40a5c1
--- /dev/null
+++ b/.github/sympy-parity/probes.json
@@ -0,0 +1,402 @@
+[
+ {
+  "area": "simplify",
+  "what": "trigonometric simplification",
+  "sympy": "simplify(sin(x)**2 + cos(x)**2)"
+ },
+ {
+  "area": "simplify",
+  "what": "rational function cancellation",
+  "sympy": "cancel((x**2-1)/(x-1))"
+ },
+ {
+  "area": "simplify",
+  "what": "radical denesting",
+  "sympy": "sqrtdenest(sqrt(5+2*sqrt(6)))"
+ },
+ {
+  "area": "simplify",
+  "what": "combine logarithms",
+  "sympy": "logcombine(log(x)+log(y), force=True)"
+ },
+ {
+  "area": "simplify",
+  "what": "hypergeometric simplification",
+  "sympy": "simplify(gamma(x+1)/gamma(x))"
+ },
+ {
+  "area": "simplify",
+  "what": "trigsimp of a product",
+  "sympy": "trigsimp(2*sin(x)*cos(x))"
+ },
+ {
+  "area": "polys",
+  "what": "factor over the rationals",
+  "sympy": "factor(x**4-1)"
+ },
+ {
+  "area": "polys",
+  "what": "multivariate gcd",
+  "sympy": "gcd(x**2*y+x*y**2, x*y)"
+ },
+ {
+  "area": "polys",
+  "what": "resultant",
+  "sympy": "resultant(x**2-1, x**2-4, x)"
+ },
+ {
+  "area": "polys",
+  "what": "square-free decomposition",
+  "sympy": "sqf(x**4-2*x**2+1)"
+ },
+ {
+  "area": "polys",
+  "what": "Groebner basis",
+  "sympy": "str(groebner([x**2+y**2-1, x-y], x, y))"
+ },
+ {
+  "area": "polys",
+  "what": "factor over a finite field",
+  "sympy": "factor(x**2+1, modulus=2)"
+ },
+ {
+  "area": "polys",
+  "what": "partial fractions",
+  "sympy": "apart(1/(x**2-1))"
+ },
+ {
+  "area": "polys",
+  "what": "minimal polynomial of an algebraic number",
+  "sympy": "minimal_polynomial(sqrt(2)+sqrt(3), x)"
+ },
+ {
+  "area": "solvers",
+  "what": "polynomial equation",
+  "sympy": "solveset(x**2-4, x)"
+ },
+ {
+  "area": "solvers",
+  "what": "transcendental equation",
+  "sympy": "solveset(exp(x)-2, x, domain=S.Reals)"
+ },
+ {
+  "area": "solvers",
+  "what": "linear system",
+  "sympy": "linsolve([x+y-3, x-y-1], [x, y])"
+ },
+ {
+  "area": "solvers",
+  "what": "nonlinear system",
+  "sympy": "nonlinsolve([x**2+y**2-1, x-y], [x, y])"
+ },
+ {
+  "area": "solvers",
+  "what": "ODE",
+  "sympy": "str(dsolve(Derivative(Function('f')(x),x)+Function('f')(x)-x, Function('f')(x)))"
+ },
+ {
+  "area": "solvers",
+  "what": "PDE",
+  "sympy": "'pdsolve exists: ' + str(callable(pdsolve))"
+ },
+ {
+  "area": "solvers",
+  "what": "diophantine equation",
+  "sympy": "diophantine(x**2+y**2-25)"
+ },
+ {
+  "area": "solvers",
+  "what": "recurrence relation",
+  "sympy": "rsolve(Function('y')(n+1)-Function('y')(n)-1, Function('y')(n))"
+ },
+ {
+  "area": "solvers",
+  "what": "inequality",
+  "sympy": "solveset(x**2-4>0, x, domain=S.Reals)"
+ },
+ {
+  "area": "solvers",
+  "what": "numeric root",
+  "sympy": "nsolve(cos(x)-x, x, 1)"
+ },
+ {
+  "area": "integrals",
+  "what": "indefinite integral",
+  "sympy": "integrate(x**2, x)"
+ },
+ {
+  "area": "integrals",
+  "what": "definite integral",
+  "sympy": "integrate(x**2, (x, 0, 1))"
+ },
+ {
+  "area": "integrals",
+  "what": "integral needing Risch",
+  "sympy": "integrate(exp(-x**2), x)"
+ },
+ {
+  "area": "integrals",
+  "what": "integration by parts",
+  "sympy": "integrate(x*exp(x), x)"
+ },
+ {
+  "area": "integrals",
+  "what": "trigonometric substitution",
+  "sympy": "integrate(1/sqrt(1-x**2), x)"
+ },
+ {
+  "area": "integrals",
+  "what": "Laplace transform",
+  "sympy": "str(laplace_transform(exp(-x), x, s))"
+ },
+ {
+  "area": "integrals",
+  "what": "Fourier transform",
+  "sympy": "str(fourier_transform(exp(-x**2), x, k))"
+ },
+ {
+  "area": "integrals",
+  "what": "multiple integral",
+  "sympy": "integrate(x*y, (x, 0, 1), (y, 0, 1))"
+ },
+ {
+  "area": "series",
+  "what": "limit",
+  "sympy": "limit(sin(x)/x, x, 0)"
+ },
+ {
+  "area": "series",
+  "what": "one-sided limit",
+  "sympy": "limit(1/x, x, 0, '+')"
+ },
+ {
+  "area": "series",
+  "what": "Taylor series",
+  "sympy": "series(sin(x), x, 0, 6)"
+ },
+ {
+  "area": "series",
+  "what": "asymptotic expansion",
+  "sympy": "series(exp(x)/x, x, oo, 3)"
+ },
+ {
+  "area": "series",
+  "what": "residue",
+  "sympy": "residue(1/x, x, 0)"
+ },
+ {
+  "area": "concrete",
+  "what": "symbolic summation",
+  "sympy": "summation(k, (k, 1, n))"
+ },
+ {
+  "area": "concrete",
+  "what": "infinite series",
+  "sympy": "summation(1/k**2, (k, 1, oo))"
+ },
+ {
+  "area": "concrete",
+  "what": "symbolic product",
+  "sympy": "product(k, (k, 1, n))"
+ },
+ {
+  "area": "matrices",
+  "what": "determinant",
+  "sympy": "Matrix([[x,1],[2,y]]).det()"
+ },
+ {
+  "area": "matrices",
+  "what": "inverse",
+  "sympy": "Matrix([[1,2],[3,4]]).inv()"
+ },
+ {
+  "area": "matrices",
+  "what": "eigenvalues",
+  "sympy": "str(Matrix([[1,2],[2,1]]).eigenvals())"
+ },
+ {
+  "area": "matrices",
+  "what": "row echelon form",
+  "sympy": "str(Matrix([[1,2],[2,4]]).rref())"
+ },
+ {
+  "area": "matrices",
+  "what": "characteristic polynomial",
+  "sympy": "Matrix([[1,2],[2,1]]).charpoly().as_expr()"
+ },
+ {
+  "area": "matrices",
+  "what": "Jordan form",
+  "sympy": "str(Matrix([[1,1],[0,1]]).jordan_form()[1])"
+ },
+ {
+  "area": "ntheory",
+  "what": "integer factorisation",
+  "sympy": "str(factorint(360))"
+ },
+ {
+  "area": "ntheory",
+  "what": "primality",
+  "sympy": "isprime(104729)"
+ },
+ {
+  "area": "ntheory",
+  "what": "Euler totient",
+  "sympy": "totient(36)"
+ },
+ {
+  "area": "ntheory",
+  "what": "Moebius",
+  "sympy": "mobius(30)"
+ },
+ {
+  "area": "ntheory",
+  "what": "Chinese remainder",
+  "sympy": "str(sympy.ntheory.modular.crt([3,5,7],[2,3,2]))"
+ },
+ {
+  "area": "ntheory",
+  "what": "continued fraction",
+  "sympy": "str(list(continued_fraction_iterator(Rational(415,93))))"
+ },
+ {
+  "area": "logic",
+  "what": "boolean simplification",
+  "sympy": "simplify_logic((a & b) | (a & ~b))"
+ },
+ {
+  "area": "logic",
+  "what": "satisfiability",
+  "sympy": "str(satisfiable(a & ~a))"
+ },
+ {
+  "area": "logic",
+  "what": "conjunctive normal form",
+  "sympy": "to_cnf((a & b) | c)"
+ },
+ {
+  "area": "sets",
+  "what": "interval intersection",
+  "sympy": "Interval(0,2).intersect(Interval(1,3))"
+ },
+ {
+  "area": "sets",
+  "what": "set-builder",
+  "sympy": "str(ConditionSet(x, x>0, S.Reals))"
+ },
+ {
+  "area": "sets",
+  "what": "image of a set",
+  "sympy": "str(imageset(Lambda(x, x**2), Interval(0,2)))"
+ },
+ {
+  "area": "functions",
+  "what": "gamma",
+  "sympy": "gamma(Rational(1,2))"
+ },
+ {
+  "area": "functions",
+  "what": "zeta",
+  "sympy": "zeta(2)"
+ },
+ {
+  "area": "functions",
+  "what": "error function",
+  "sympy": "erf(0)"
+ },
+ {
+  "area": "functions",
+  "what": "Bessel",
+  "sympy": "besselj(0,0)"
+ },
+ {
+  "area": "functions",
+  "what": "Lambert W",
+  "sympy": "LambertW(0)"
+ },
+ {
+  "area": "functions",
+  "what": "polylogarithm",
+  "sympy": "polylog(2,1)"
+ },
+ {
+  "area": "assumptions",
+  "what": "assumption-driven simplification",
+  "sympy": "sqrt(Symbol('p', positive=True)**2)"
+ },
+ {
+  "area": "assumptions",
+  "what": "ask a predicate",
+  "sympy": "str(ask(Q.positive(Symbol('p', positive=True))))"
+ },
+ {
+  "area": "printing",
+  "what": "LaTeX output",
+  "sympy": "latex(Integral(x**2, x))"
+ },
+ {
+  "area": "printing",
+  "what": "C code generation",
+  "sympy": "ccode(x**2+sin(x))"
+ },
+ {
+  "area": "printing",
+  "what": "MathML output",
+  "sympy": "'mathml exists: ' + str(callable(mathml))"
+ },
+ {
+  "area": "discrete",
+  "what": "convolution",
+  "sympy": "str(convolution([1,2],[3,4]))"
+ },
+ {
+  "area": "discrete",
+  "what": "fast Fourier transform",
+  "sympy": "str(fft([1,2,3,4]))"
+ },
+ {
+  "area": "combinatorics",
+  "what": "permutation group",
+  "sympy": "'PermutationGroup exists: ' + str(callable(PermutationGroup))"
+ },
+ {
+  "area": "stats",
+  "what": "probability of an event",
+  "sympy": "str(P(Die('d',6) > 3))"
+ },
+ {
+  "area": "geometry",
+  "what": "line intersection",
+  "sympy": "str(Line((0,0),(1,1)).intersection(Line((0,1),(1,0))))"
+ },
+ {
+  "area": "physics",
+  "what": "units",
+  "sympy": "'units exists: ' + str(callable(convert_to))"
+ },
+ {
+  "area": "vector",
+  "what": "vector calculus",
+  "sympy": "'CoordSys3D exists: ' + str(callable(CoordSys3D))"
+ },
+ {
+  "area": "tensor",
+  "what": "tensor algebra",
+  "sympy": "'tensor exists: ' + str(callable(TensorIndexType))"
+ },
+ {
+  "area": "diffgeom",
+  "what": "differential geometry",
+  "sympy": "'Manifold exists: ' + str(callable(Manifold))"
+ },
+ {
+  "area": "holonomic",
+  "what": "holonomic functions",
+  "sympy": "'expr_to_holonomic exists: ' + str(callable(expr_to_holonomic))"
+ },
+ {
+  "area": "codegen",
+  "what": "compile to a callable",
+  "sympy": "str(lambdify(x, x**2)(3))"
+ }
+]
\ No newline at end of file
diff --git a/.github/sympy-parity/sympy-baseline.json b/.github/sympy-parity/sympy-baseline.json
new file mode 100644
index 000000000..d66ee67f4
--- /dev/null
+++ b/.github/sympy-parity/sympy-baseline.json
@@ -0,0 +1,406 @@
+{
+ "sympy": "1.14.0",
+ "measured": "2026-08-23",
+ "answers": [
+  {
+   "area": "simplify",
+   "what": "trigonometric simplification",
+   "sympy": "1"
+  },
+  {
+   "area": "simplify",
+   "what": "rational function cancellation",
+   "sympy": "x + 1"
+  },
+  {
+   "area": "simplify",
+   "what": "radical denesting",
+   "sympy": "sqrt(2) + sqrt(3)"
+  },
+  {
+   "area": "simplify",
+   "what": "combine logarithms",
+   "sympy": "log(x*y)"
+  },
+  {
+   "area": "simplify",
+   "what": "hypergeometric simplification",
+   "sympy": "x"
+  },
+  {
+   "area": "simplify",
+   "what": "trigsimp of a product",
+   "sympy": "sin(2*x)"
+  },
+  {
+   "area": "polys",
+   "what": "factor over the rationals",
+   "sympy": "(x - 1)*(x + 1)*(x**2 + 1)"
+  },
+  {
+   "area": "polys",
+   "what": "multivariate gcd",
+   "sympy": "x*y"
+  },
+  {
+   "area": "polys",
+   "what": "resultant",
+   "sympy": "9"
+  },
+  {
+   "area": "polys",
+   "what": "square-free decomposition",
+   "sympy": "(x**2 - 1)**2"
+  },
+  {
+   "area": "polys",
+   "what": "Groebner basis",
+   "sympy": "GroebnerBasis([x - y, 2*y**2 - 1], x, y, domain='ZZ', order='lex')"
+  },
+  {
+   "area": "polys",
+   "what": "factor over a finite field",
+   "sympy": "(x + 1)**2"
+  },
+  {
+   "area": "polys",
+   "what": "partial fractions",
+   "sympy": "-1/(2*(x + 1)) + 1/(2*(x - 1))"
+  },
+  {
+   "area": "polys",
+   "what": "minimal polynomial of an algebraic number",
+   "sympy": "x**4 - 10*x**2 + 1"
+  },
+  {
+   "area": "solvers",
+   "what": "polynomial equation",
+   "sympy": "{-2, 2}"
+  },
+  {
+   "area": "solvers",
+   "what": "transcendental equation",
+   "sympy": "{log(2)}"
+  },
+  {
+   "area": "solvers",
+   "what": "linear system",
+   "sympy": "{(2, 1)}"
+  },
+  {
+   "area": "solvers",
+   "what": "nonlinear system",
+   "sympy": "{(-sqrt(2)/2, -sqrt(2)/2), (sqrt(2)/2, sqrt(2)/2)}"
+  },
+  {
+   "area": "solvers",
+   "what": "ODE",
+   "sympy": "Eq(f(x), C1*exp(-x) + x - 1)"
+  },
+  {
+   "area": "solvers",
+   "what": "PDE",
+   "sympy": "pdsolve exists: True"
+  },
+  {
+   "area": "solvers",
+   "what": "diophantine equation",
+   "sympy": "{(-3, 4), (3, 4), (4, -3), (4, 3), (-4, -3), (-5, 0), (-4, 3), (5, 0), (0, 5), (-3, -4), (0, -5), (3, -4)}"
+  },
+  {
+   "area": "solvers",
+   "what": "recurrence relation",
+   "sympy": "C0 + n"
+  },
+  {
+   "area": "solvers",
+   "what": "inequality",
+   "sympy": "Union(Interval.open(-oo, -2), Interval.open(2, oo))"
+  },
+  {
+   "area": "solvers",
+   "what": "numeric root",
+   "sympy": "0.739085133215161"
+  },
+  {
+   "area": "integrals",
+   "what": "indefinite integral",
+   "sympy": "x**3/3"
+  },
+  {
+   "area": "integrals",
+   "what": "definite integral",
+   "sympy": "1/3"
+  },
+  {
+   "area": "integrals",
+   "what": "integral needing Risch",
+   "sympy": "sqrt(pi)*erf(x)/2"
+  },
+  {
+   "area": "integrals",
+   "what": "integration by parts",
+   "sympy": "(x - 1)*exp(x)"
+  },
+  {
+   "area": "integrals",
+   "what": "trigonometric substitution",
+   "sympy": "asin(x)"
+  },
+  {
+   "area": "integrals",
+   "what": "Laplace transform",
+   "sympy": "(1/(s + 1), -1, True)"
+  },
+  {
+   "area": "integrals",
+   "what": "Fourier transform",
+   "sympy": "sqrt(pi)*exp(-pi**2*k**2)"
+  },
+  {
+   "area": "integrals",
+   "what": "multiple integral",
+   "sympy": "1/4"
+  },
+  {
+   "area": "series",
+   "what": "limit",
+   "sympy": "1"
+  },
+  {
+   "area": "series",
+   "what": "one-sided limit",
+   "sympy": "oo"
+  },
+  {
+   "area": "series",
+   "what": "Taylor series",
+   "sympy": "x - x**3/6 + x**5/120 + O(x**6)"
+  },
+  {
+   "area": "series",
+   "what": "asymptotic expansion",
+   "sympy": "exp(x)/x"
+  },
+  {
+   "area": "series",
+   "what": "residue",
+   "sympy": "1"
+  },
+  {
+   "area": "concrete",
+   "what": "symbolic summation",
+   "sympy": "n**2/2 + n/2"
+  },
+  {
+   "area": "concrete",
+   "what": "infinite series",
+   "sympy": "pi**2/6"
+  },
+  {
+   "area": "concrete",
+   "what": "symbolic product",
+   "sympy": "factorial(n)"
+  },
+  {
+   "area": "matrices",
+   "what": "determinant",
+   "sympy": "x*y - 2"
+  },
+  {
+   "area": "matrices",
+   "what": "inverse",
+   "sympy": "Matrix([[-2, 1], [3/2, -1/2]])"
+  },
+  {
+   "area": "matrices",
+   "what": "eigenvalues",
+   "sympy": "{3: 1, -1: 1}"
+  },
+  {
+   "area": "matrices",
+   "what": "row echelon form",
+   "sympy": "(Matrix([\n[1, 2],\n[0, 0]]), (0,))"
+  },
+  {
+   "area": "matrices",
+   "what": "characteristic polynomial",
+   "sympy": "lambda**2 - 2*lambda - 3"
+  },
+  {
+   "area": "matrices",
+   "what": "Jordan form",
+   "sympy": "Matrix([[1, 1], [0, 1]])"
+  },
+  {
+   "area": "ntheory",
+   "what": "integer factorisation",
+   "sympy": "{2: 3, 3: 2, 5: 1}"
+  },
+  {
+   "area": "ntheory",
+   "what": "primality",
+   "sympy": "True"
+  },
+  {
+   "area": "ntheory",
+   "what": "Euler totient",
+   "sympy": "12"
+  },
+  {
+   "area": "ntheory",
+   "what": "Moebius",
+   "sympy": "-1"
+  },
+  {
+   "area": "ntheory",
+   "what": "Chinese remainder",
+   "sympy": "(23, 105)"
+  },
+  {
+   "area": "ntheory",
+   "what": "continued fraction",
+   "sympy": "[4, 2, 6, 7]"
+  },
+  {
+   "area": "logic",
+   "what": "boolean simplification",
+   "sympy": "a"
+  },
+  {
+   "area": "logic",
+   "what": "satisfiability",
+   "sympy": "False"
+  },
+  {
+   "area": "logic",
+   "what": "conjunctive normal form",
+   "sympy": "(a | c) & (b | c)"
+  },
+  {
+   "area": "sets",
+   "what": "interval intersection",
+   "sympy": "Interval(1, 2)"
+  },
+  {
+   "area": "sets",
+   "what": "set-builder",
+   "sympy": "ConditionSet(x, x > 0, Reals)"
+  },
+  {
+   "area": "sets",
+   "what": "image of a set",
+   "sympy": "Interval(0, 4)"
+  },
+  {
+   "area": "functions",
+   "what": "gamma",
+   "sympy": "sqrt(pi)"
+  },
+  {
+   "area": "functions",
+   "what": "zeta",
+   "sympy": "pi**2/6"
+  },
+  {
+   "area": "functions",
+   "what": "error function",
+   "sympy": "0"
+  },
+  {
+   "area": "functions",
+   "what": "Bessel",
+   "sympy": "1"
+  },
+  {
+   "area": "functions",
+   "what": "Lambert W",
+   "sympy": "0"
+  },
+  {
+   "area": "functions",
+   "what": "polylogarithm",
+   "sympy": "pi**2/6"
+  },
+  {
+   "area": "assumptions",
+   "what": "assumption-driven simplification",
+   "sympy": "p"
+  },
+  {
+   "area": "assumptions",
+   "what": "ask a predicate",
+   "sympy": "True"
+  },
+  {
+   "area": "printing",
+   "what": "LaTeX output",
+   "sympy": "\\int x^{2}\\, dx"
+  },
+  {
+   "area": "printing",
+   "what": "C code generation",
+   "sympy": "pow(x, 2) + sin(x)"
+  },
+  {
+   "area": "printing",
+   "what": "MathML output",
+   "sympy": "mathml exists: True"
+  },
+  {
+   "area": "discrete",
+   "what": "convolution",
+   "sympy": "[3, 10, 8]"
+  },
+  {
+   "area": "discrete",
+   "what": "fast Fourier transform",
+   "sympy": "[10, -2 - 2*I, -2, -2 + 2*I]"
+  },
+  {
+   "area": "combinatorics",
+   "what": "permutation group",
+   "sympy": "PermutationGroup exists: True"
+  },
+  {
+   "area": "stats",
+   "what": "probability of an event",
+   "sympy": "1/2"
+  },
+  {
+   "area": "geometry",
+   "what": "line intersection",
+   "sympy": "[Point2D(1/2, 1/2)]"
+  },
+  {
+   "area": "physics",
+   "what": "units",
+   "sympy": "units exists: True"
+  },
+  {
+   "area": "vector",
+   "what": "vector calculus",
+   "sympy": "CoordSys3D exists: True"
+  },
+  {
+   "area": "tensor",
+   "what": "tensor algebra",
+   "sympy": "tensor exists: True"
+  },
+  {
+   "area": "diffgeom",
+   "what": "differential geometry",
+   "sympy": "Manifold exists: True"
+  },
+  {
+   "area": "holonomic",
+   "what": "holonomic functions",
+   "sympy": "expr_to_holonomic exists: True"
+  },
+  {
+   "area": "codegen",
+   "what": "compile to a callable",
+   "sympy": "9"
+  }
+ ]
+}
\ No newline at end of file
diff --git a/.github/workflows/CPPTest.yml b/.github/workflows/CPPTest.yml
index 47820ae58..6ea49ca29 100644
--- a/.github/workflows/CPPTest.yml
+++ b/.github/workflows/CPPTest.yml
@@ -1,5 +1,18 @@
 name: 'C++ Test'
 
+# Coverage for this leg is https://github.com/asc-community/AngouriMath/issues/746 item 35 and
+# https://github.com/asc-community/AngouriMath/issues/397, and it is not the same measurement the
+# other three legs make. C#, F# and Interactive each run `dotnet test` and collect coverage of
+# managed assemblies with the VSTest XPlat collector. Nothing in this job can: the library reaches
+# the tests as a Native AOT shared object, the published artefact contains no managed assembly at
+# all, and ctest runs a native executable rather than a test host the collector can attach to.
+#
+# What this leg does have, and what nothing measured until now, is the C++ of the wrapper itself:
+# Sources/Wrappers/AngouriMath.CPP.Importing holds the handle cache, the error-code translation
+# and the exception boundary that every C++ caller goes through. That is ordinary C++, and gcov
+# covers it, so the Linux job builds with --coverage and uploads Cobertura under the `cpp` flag.
+# The number is coverage of the C++ wrapper, not of AngouriMath through it; see
+# Sources/AngouriMath/Docs/Contributing/Coverage.md.
 on:
   push:
     branches:
@@ -38,8 +51,8 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         sudo apt-get update
-        sudo apt-get install -y zlib1g-dev libkrb5-dev
-        
+        sudo apt-get install -y zlib1g-dev libkrb5-dev gcovr
+
     - name: 'Building the library into native for Windows'
       if: ${{ matrix.os == 'windows-latest' }}
       shell: cmd
@@ -57,11 +70,21 @@ jobs:
         dotnet publish -p:NativeLib=Shared -p:SelfContained=true -r ${{ matrix.flag }} -c release        
 
     - name: 'Preparing tests'
+      if: ${{ matrix.os != 'ubuntu-latest' }}
       run: |
         cd ${{ env.tests_path }}
         cmake -S . -B build
         mkdir "build/Debug"
 
+    # Same configure, plus gcov instrumentation. -O0 so that a line in the report is a line in the
+    # source; the tests take under a second either way.
+    - name: 'Preparing tests, instrumented for coverage'
+      if: ${{ matrix.os == 'ubuntu-latest' }}
+      run: |
+        cd ${{ env.tests_path }}
+        cmake -S . -B build -DCMAKE_BUILD_TYPE=Debug -DCMAKE_CXX_FLAGS="--coverage -O0 -g" -DCMAKE_EXE_LINKER_FLAGS="--coverage"
+        mkdir "build/Debug"
+
     - name: 'Copying the library for Windows'
       if: ${{ matrix.os == 'windows-latest' }}
       run: |
@@ -88,6 +111,34 @@ jobs:
         cd build
         ctest
       
+    # The filter is what keeps this honest: without it the report is dominated by the vendored
+    # googletest sources, which nobody here maintains, and by RunTests.cpp, which is 100% covered
+    # by construction because it is the thing being run.
+    - name: 'Coverage of the C++ wrapper'
+      if: ${{ matrix.os == 'ubuntu-latest' }}
+      run: |
+        cd ${{ env.tests_path }}/build
+        gcovr --root ../../../../.. --filter '.*AngouriMath\.CPP\.Importing.*' \
+              --print-summary --cobertura coverage-cpp.xml
+        # gcovr writes a well-formed report with nothing in it when the filter matches no file, and
+        # codecov accepts it. Fail here instead, so that a flag claiming this leg is covered cannot
+        # go green on an empty report.
+        grep -q 'lines-valid="[1-9]' coverage-cpp.xml
+
+    # `files`, not `file`: v5 of the action does not take `file`, and given it, it warns about an
+    # unexpected input and then searches the workspace itself. Measured here, that search found 96
+    # reports -- the whole of googletest and of /usr/include/c++/13 -- so the flag would have
+    # carried everything the filter above exists to exclude. `disable_search` and `plugins: noop`
+    # keep it to the one report gcovr wrote.
+    - name: 'Send to codecov'
+      if: ${{ matrix.os == 'ubuntu-latest' }}
+      uses: codecov/codecov-action@v5
+      with:
+        files: ./Sources/Tests/CPPWrapperUnitTests/tests/build/coverage-cpp.xml
+        disable_search: true
+        plugins: noop
+        flags: cpp
+
     - name: 'Upload logs'
       if: failure()
       uses: actions/upload-artifact@v6
diff --git a/.github/workflows/CSharpTest.yml b/.github/workflows/CSharpTest.yml
index f825535b6..b73e31615 100644
--- a/.github/workflows/CSharpTest.yml
+++ b/.github/workflows/CSharpTest.yml
@@ -84,9 +84,12 @@ jobs:
         if-no-files-found: warn
         retention-days: 90
 
+    # `files`, not `file`. v5 of the action does not take `file`: it warns about an unexpected
+    # input, ignores the path, and searches the workspace instead, so the named report was never
+    # what made this flag.
     - name: 'Send to codecov'
       if: ${{ matrix.os == 'windows-latest' }}
       uses: codecov/codecov-action@v5
       with:
-        file: ./Sources/Tests/UnitTests/coverage.opencover.xml
+        files: ./Sources/Tests/UnitTests/coverage.opencover.xml
         flags: csharp
diff --git a/.github/workflows/SymPyParity.yml b/.github/workflows/SymPyParity.yml
new file mode 100644
index 000000000..91b1e16b6
--- /dev/null
+++ b/.github/workflows/SymPyParity.yml
@@ -0,0 +1,76 @@
+name: 'SymPy parity watch'
+
+# The comparison against SymPy is a snapshot: 80 probes over 23 areas, taken by hand, against one
+# SymPy version, and published in Sources/AngouriMath/Docs/Usage/Comparison.md. A snapshot cannot
+# tell you it has gone out of date, and this is
+# https://github.com/asc-community/AngouriMath/issues/746 item 73 -- "the one SymPy snapshot goes
+# stale on its own -- it wants a scheduled re-run rather than another hand-taken snapshot".
+#
+# What runs here is the SymPy half only. The other half asks this library the same 80 questions
+# and lives in the analysis workspace, not in this repository, so this job compares nothing
+# between the two libraries and repeats no verdict about either. It answers the one question the
+# published page cannot answer about itself: is SymPy still the version, and still giving the
+# answers, that the page was written from. When either has moved the job fails, and the run
+# uploads the new answers so re-taking the snapshot is a file copy plus a re-run of the harness.
+#
+# Maxima and SageMath are named alongside SymPy in that roadmap item and are deliberately absent
+# here. Nothing in this repository has ever measured them, and a workflow that installed them
+# without a baseline to compare against would report a green tick for a comparison nobody made.
+#
+# The path filters carry no leading "./" -- a filter is a glob relative to the repository root,
+# and "./x" matches nothing there. That defect silenced Benchmark.yml for ninety-eight merged
+# pull requests and CPPBuild.yml for its whole existence
+# (https://github.com/asc-community/AngouriMath/issues/775). The filters name this file and the
+# probe set so that a change to either demonstrates itself on its own pull request; the weekly
+# run is what actually watches SymPy, and a schedule trigger only fires once the workflow is on
+# the default branch.
+on:
+  schedule:
+    - cron: '17 6 * * 1'
+  workflow_dispatch:
+  push:
+    branches:
+      - master
+    paths:
+      - '.github/sympy-parity/**'
+      - '.github/workflows/SymPyParity.yml'
+  pull_request:
+    branches:
+      - '*'
+    paths:
+      - '.github/sympy-parity/**'
+      - '.github/workflows/SymPyParity.yml'
+
+permissions:
+  contents: read
+
+jobs:
+  Parity:
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v6
+
+    - uses: actions/setup-python@v5
+      with:
+        python-version: '3.x'
+
+    # Unpinned on purpose: the question this job asks is what today's SymPy does. Pinning it
+    # would make the job agree with the baseline for ever and measure nothing.
+    - name: 'Install the newest SymPy'
+      run: |
+        python -m pip install --upgrade pip
+        python -m pip install --upgrade sympy
+        python -c "import sympy; print('SymPy', sympy.__version__)"
+
+    - name: 'Re-run the SymPy probes against the published snapshot'
+      run: python .github/sympy-parity/check.py
+
+    - name: 'Upload what SymPy answers today'
+      if: always()
+      uses: actions/upload-artifact@v6
+      with:
+        name: sympy-current
+        path: .github/sympy-parity/sympy-current.json
+        if-no-files-found: warn
+        retention-days: 90
diff --git a/Sources/AngouriMath/Docs/Contributing/Coverage.md b/Sources/AngouriMath/Docs/Contributing/Coverage.md
new file mode 100644
index 000000000..518e5740d
--- /dev/null
+++ b/Sources/AngouriMath/Docs/Contributing/Coverage.md
@@ -0,0 +1,59 @@
+# Coverage
+
+Four legs of CI run tests, and they do not all measure the same thing. This says what each one
+measures, so that a missing number is read as the boundary it is rather than as a leg nobody got
+round to instrumenting.
+
+| leg | workflow | what it runs | how coverage is collected | codecov flag |
+|---|---|---|---|---|
+| C# | `CSharpTest.yml` | `dotnet test Sources/Tests/UnitTests` | coverlet, both as `/p:CollectCoverage=true` and as the VSTest `XPlat Code Coverage` collector | `csharp` |
+| F# | `FSharpTest.yml` | `dotnet test Sources/Tests/FSharpWrapperUnitTests` | VSTest `XPlat Code Coverage` collector | `fsharp` |
+| Interactive | `InteractiveTest.yml` | `dotnet test` over the Interactive wrapper and the terminal | VSTest `XPlat Code Coverage` collector | `interactive` |
+| C++ | `CPPTest.yml` | `ctest` over a native executable | `gcov`, reported by `gcovr` as Cobertura | `cpp` |
+
+The first three all take the same route: `dotnet test` hosts managed assemblies, and a data
+collector instruments them in that host. The fourth cannot, and the reason is worth writing down
+because it decides what the `cpp` number means.
+
+## The C++ leg does not carry managed code into CI
+
+`AngouriMath.CPP.Exporting` sets `PublishAot`, and `CPPTest.yml` publishes it with
+`-p:NativeLib=Shared -p:SelfContained=true`. What arrives in
+`AngouriMath.CPP.Importing/out-x64` is one shared object and a symbol file; there is no
+`AngouriMath.dll` and no `AngouriMath.CPP.Exporting.dll` in it. The tests are then a C++ binary
+that ctest runs, not a test host.
+
+So neither half of the mechanism the other three legs use is present. There is no `dotnet test` for
+a collector to attach to, and there is no managed assembly in the deployed artefact for a collector
+to instrument. Published Native AOT output is a native binary, and .NET documents the diagnostic
+tooling that follows managed code into it as
+[partial](https://learn.microsoft.com/dotnet/core/deploying/native-aot/diagnostics) — the managed
+debugger does not attach, and managed heap analysis is unsupported.
+
+That is a statement about this route, not a proof that no route exists. What it settles is that
+copying a step from `CSharpTest.yml` into `CPPTest.yml` would upload an empty report under a flag
+claiming the C++ leg was covered, which is worse than the flag being absent.
+
+## What the `cpp` flag does cover
+
+The C++ leg has C++ of its own, and until now nothing measured it:
+`Sources/Wrappers/AngouriMath.CPP.Importing` is the handle cache, the error-code translation and
+the exception boundary that every C++ caller passes through, and a defect there is a defect no
+managed test can see. It is ordinary C++, so `gcov` covers it.
+
+The Linux job in the matrix configures cmake with `--coverage -O0 -g`, runs the same ctest suite as
+the other two, and reports with
+
+```
+gcovr --root <repo root> --filter '.*AngouriMath\.CPP\.Importing.*' --cobertura coverage-cpp.xml
+```
+
+The filter is the honest part. Without it the report is mostly vendored googletest, and it includes
+`RunTests.cpp`, which is fully covered by construction because it is the thing being run.
+
+**So `cpp` is coverage of the C++ wrapper, not of AngouriMath through the C++ wrapper.** A line of
+`Simplify` exercised only from a C++ test appears in no report at all. Read the flag as answering
+"how much of the wrapper do these tests reach", and read `csharp` for the kernel.
+
+To reproduce it locally, publish the exporting project as `CPPTest.yml` does, then configure the
+tests with the two cmake flags above and run `gcovr` from the build directory.