Fix Bessel numerical instability and DnonDn inversion completeness

[johmathe]

Summary

Re-ran the bug assessment on v0.3.6 with property-based checks (group-invariance and inversion round-trips at tight tolerances). Four real bugs confirmed and fixed, plus small cleanups. All were masked by existing test tolerances or incomplete assertions.

1. bessel_jn forward recurrence diverges for x < n (_bessel.py)

The three-term recurrence was applied unconditionally, but it amplifies the Y_n admixture exponentially when x < n: J_20(0.5) returned -1e+11 instead of 3.7e-31. Since disk harmonics evaluate J_n(lambda*r) with r -> 0, the SO2onDisk basis matrix was rank-degenerate — rank 5 out of 804 at L=32 — making inversion and rotation-invariance meaningless at that size. Fixed by switching to Miller's backward recurrence (with even-sum normalization and overflow rescaling) for |x| < n. After the fix: L=32 basis rank 798/804, all values match mpmath to ~1e-14 relative.

2. _bisect_newton_batch walks away from converged roots (_bessel.py)

When an iterate hit J_n(x) == 0.0 exactly, the sign test fa * fx < 0 became false on every subsequent iteration, so the bisection kept moving the lower bound past the root and converged to the wrong bracket endpoint. Example: j_{6,5} returned 25.430 (which is j_{5,6}) instead of 23.586. This silently corrupted the interlacing brackets for all higher orders — dozens of roots were off by ~2.0. Fixed by pinning the bracket when the root is hit exactly and collapsing brackets whose endpoint is already a root. Added interlacing (j_{n-1,k} < j_{n,k} < j_{n-1,k+1}) and mpmath-reference tests that fail loudly on the old behavior.

3. Newton denominator copysign bug (_bessel.py)

dfx.clamp_min(1e-30).copysign(dfx) maps every negative derivative to -1e-30, producing huge rejected Newton steps and degrading the solver to pure bisection for about half the roots. Now dfx.abs().clamp_min(1e-30).copysign(dfx).

4. DnonDn.invert broken for even n and n=3 (dn_on_dn.py)

forward(invert(beta)) differed from beta with relative errors of order 1 for all even n and n=3 (the docstring promised recovery up to a D_n action). Root cause: the symmetric-sqrt seed for F(rho_1) is only determined up to an O(2) twist, and the fold entries of the extraction chain — rho_02/rho_03 for even n (from rho_1 x rho_{n/2-1}), and the folded rho_1 block at n=3 — are mixed by a rotation R(h*theta) that depends on the unknown twist angle. The old code wrote these gauge-twisted values over exact coefficients.

The fix:

• the chain now only writes forward-chain 2D blocks (rho_{k+1}) and rho_01, never overwriting exact seeds;
• the fold rotation angle is measured from the fold block, the seed is corrected over a small finite candidate set (rotation shifts x reflection), and the candidate whose reconstruction reproduces beta is selected — the winner is certified via forward, so the result is exact rather than approximate;
• odd n >= 5 has no fold and keeps the single-pass path.

Round-trip is now exact to ~1e-11 relative across n=3..32 (5 seeds x batch 64), and even-n reconstructions land in the D_n orbit of the input (~1e-14). The one caveat, documented in the docstring: signals with near-singular F(rho_1) are intrinsically ill-conditioned for this algorithm (worst observed 7e-7 relative on one such sample).

Minor

• Removed dead expressions (r1 - r0, f.shape[0], 2 * l2 + 1) in octa_on_octa.py / _cg.py.
• so3_on_s2.py sparse-cache build used the global RNG for its randperm; now uses a private seeded generator so building the cache no longer perturbs the caller's RNG stream (test added).
• DnonDn.forward now validates input shape and rejects complex input, matching TorusOnTorus. (CnonCn intentionally keeps accepting complex signals — its invert returns complex reconstructions that are re-verified through forward; test documents this.)

Test plan

• New tests: full-bispectrum roundtrip for DnonDn (n=3..16), orbit recovery for even n, mpmath-referenced Bessel values in the x < n regime, high-order reference roots, root interlacing, DHT basis rank, RNG isolation for SO3onS2, input validation.
• Full suite: 606 passed, 22 skipped (CUDA).
• Pre-commit hooks pass (ruff, ruff-format, docformatter, bandit, codespell).

Made with Cursor

[codecov]

Codecov Report

❌ Patch coverage is 96.66667% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 95.46%. Comparing base (183431f) to head (e836eff).

Files with missing lines	Patch %	Lines
src/bispectrum/_bessel.py	91.22%	5 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #49      +/-   ##
==========================================
- Coverage   95.50%   95.46%   -0.04%     
==========================================
  Files          11       11              
  Lines        2135     2227      +92     
==========================================
+ Hits         2039     2126      +87     
- Misses         96      101       +5     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.