fix(mm): stable tQSSA free-substrate root and subtraction-free partials (#89)

[wshlavacek]

Closes #89.

The bug, and the half of it the issue did not have

The free substrate was sFree = 0.5*(delta + D) with delta = S - Km - E and
D = sqrt(delta^2 + 4*Km*S). For delta < 0 that subtracts two nearly-equal
positive numbers and loses ~2 significant digits per decade of
|delta|/sqrt(4*Km*S) — no correct digit by 1e8. sFree is the positive root of
x^2 - delta*x - Km*S = 0, so the delta < 0 branch is now the conjugate form
2*Km*S/(D - delta). delta >= 0 keeps the textbook expression bit-for-bit.

Applying only that — the fix as written in the issue — is not enough, and the
measurement is what says so. The derivatives were the chain rule through sFree
(dsF/dE = (-1 - delta/D)/2, and similarly for S and Km), and those cancel in
exactly the same regime. With a correct sFree in hand, d(rate)/dE still
sat at 1e+10 relative error.

On the regression fixture the cancellation goes past magnitude and flips the
sign: d(rate)/dE comes out -3.83e-05 where the true value is +4.96e-15. A
negative d(rate)/dE says adding enzyme slows the reaction, and that entry drives
CVODE's Newton step — so it steps confidently the wrong way rather than merely
converging slowly. The C++ test asserts that sign on its own line, so the sharpest
symptom is the first thing to fail.

The forms that ship

Differentiating the symmetric form of the rate collapses each partial to a
single subtraction-free quotient. The tQSSA complex is c = (A - D)/2 with
A = E + S + Km, and that same D is sqrt(A^2 - 4*E*S), so dc/dE = (S - c)/D
with S - c = sFree:

d(rate)/dE  =  kcat*stat*sFree/D
d(rate)/dS  =  kcat*stat*E*Km/((Km + sFree)*D)
d(rate)/dKm = -rate/D

Each is identical to sympy.diff of the shipped rate (asserted in the suite, not
just checked once). Measured against mpmath at 60 digits over the four regime
sweeps from the issue — worst relative error in float64, before → after:

sweep	sFree	rate	d/dE	d/dS	d/dKm
uniform, O(1)	4.5e-14 → 4.1e-16	4.3e-14 → 4.9e-16	5.6e-13 → 6.1e-16	2.6e-15 → 7.3e-16	1.4e-12 → 7.6e-16
log-spread 1e-4..1e3	6.2e-04 → 3.5e-15	6.2e-04 → 5.8e-16	3.0e+03 → 1.8e-15	2.5e-10 → 6.1e-15	2.6e+03 → 2.4e-15
Km >> S,E	2.6e-09 → 5.4e-16	2.6e-09 → 6.5e-16	2.6e-09 → 7.6e-16	7.6e-16 → 4.7e-16	9.4e-09 → 7.9e-16
Km << S,E (deep saturation)	2.1e-02 → 2.4e-13	2.1e-02 → 5.1e-16	2.3e+09 → 4.3e-14	5.3e-05 → 4.7e-13	8.7e+09 → 2.3e-13

Sites

The three lines were open-coded at every site rather than shared, which is why the
issue asks for a grep first. They now come from one helper per language:
mm_tqssa() in include/bngsim/mm_jacobian.hpp, which compute_rxn_rate and
mm_tqssa_derivatives both call, and _mm_sfree_c_lines() in
python/bngsim/_codegen.py, which every emitter calls.

Covered: the interpreted RHS and the SSA propensity emitter (src/model.cpp), the
closed-form analytical Jacobian shared by the dense and sparse CVODE paths, both
compiled-RHS emitters, the analytical Jacobian and its fused J*v, #55's analytic
df/dp, and the JAX RHS (_jax_rhs.py, where the denominator is masked before
the divide so the unselected branch cannot inject a NaN into a tangent).

_CODEGEN_VERSION → 26: MM models' emitted C changes, and a cached v25 .so would
keep serving pre-fix numbers.

Why the fixture is synthetic

The corpus cannot gate this. Both MM models (test_MM, mCaMKII_Ca_Spike) do
reach delta < 0, so they take the new branch — but at
|delta|/sqrt(4*Km*S) <= 4.95 and <= 1.32, where the two forms differ by 1–3
ulp. Largest sFree shift measured over their trajectories: 1.9e-16 and
6.0e-16. generate_rhs_c is byte-identical on 120 sampled non-MM corpus models
and differs only on those two.

So the regression fixture is a deliberately stiff one,
tests/data/mm_tqssa_stiff.net at ratio ~1e7 (large enzyme excess, small Km*S),
shared by the C++ and Python suites.

Verification

• 51/51 C++ (ctest), 2773 passed / 89 skipped Python, pre-commit clean.
• Both halves independently gated: reverting only the root fails
  test_mm_tqssa_stiff_root at rel=0.0199; reverting only the derivatives fails
  it at rel=1.36e10.
• TestNumericalStability pins each rejected form against mpmath — the textbook
  root, the chain rule through sFree, and Analytic sensitivity RHS for Functional rate laws (df/dp), so forward sensitivity leaves the CVODES difference-quotient path #55's already-rejected "simplified"
  d/dKm — so a future tidy-up fails there instead of silently in a trajectory or
  a gradient.
• TestTheEmittedCAtAStiffRatio compiles the emitted RHS and df/dp and checks
  them against mpmath rather than against a finite difference of themselves. That
  matters: a self-FD check is blind to a cancellation both sides share, which is
  how Analytic sensitivity RHS for Functional rate laws (df/dp), so forward sensitivity leaves the CVODES difference-quotient path #55's tests passed over a broken rate.
• SSA propensity library compiled and run at the stiff point: 1.6e-16 relative,
  against 2.1e-2 before.

The SBML exporter is a required companion, not an optional follow-up

_mm_formula in python/bngsim/convert/_sbml_writer.py still exports the
cancelling form, and this PR is what makes that matter: moving the engine to
the stable root while the exporter keeps 0.5*(delta + D) means an exported model
no longer agrees with the engine that exported it. Measured on this branch with
validate_conversion(..., direction="net2sbml"):

fixture	L2 round-trip	L3 numerical
mm_tqssa.net (ratio ~1)	pass, max scale-relative Δ exactly 0.0	pass, 0/303 cells
mm_tqssa_stiff.net (ratio ~1e7)	fail, Δ = 2.17e-02 > 1e-06	fail, 200/303 cells, max rel 1.0, max abs 148

So nothing that exists today regresses — the corpus MM models sit at ratio ~1 and
round-trip bit-for-bit — but a user model in the stiff regime now gets an SBML
export that is 2% off its own source, and the gate that catches it
(RHS_FAITHFUL_TOL = 1e-6 in parity_checks/bng_parity/net_roadrunner.py) fails
silently: it returns None and the caller turns that into an unscored
REFERENCE_FAILED row, so an export that drifts drops models out of parity rather
than failing anything loudly.

The companion change is in progress separately: export the same branch as a
piecewise(0.5*(delta+D), delta >= 0, 2*Km*S/(D-delta)). Two things make that
cheaper than it looks — the emitted MM law already contains max(sFree, 0),
and _validate.py's _NONSMOOTH tuple lists Max and Piecewise together, so
L4 is already inconclusive on every MM model and a new piecewise cannot
regress it; and the writer already targets L3v2 for piecewise/min-max, since every
BNGL if() exports as one.

Separately, and genuinely pre-existing: the compiled RHS has no sFree < 0 clamp
where the interpreted path does, so on a negative-concentration probe the two
backends disagree (at S = -1 the compiled RHS gives dP/dt = -0.169 where the
interpreted one gives 0) — and the emitted Jacobian, which does guard on
sFree > 0, contradicts the emitted RHS beside it in the same artifact. Not
touched here; it needs a decision about which behaviour is right rather than an
edit.