Post-0d09740 regression audit: GPU period fix, elliptic123 routing, scipy-free runtime, adversarial-review round - #36
Merged
Conversation
…iptic123 routing, poles, cancellation) Follow-up audit on baseline 0d09740 (docs/specs/post-0d09740-regression-audit.md): - gpu_elliptic12/gpu_ellipj: apply the same quasi-period reduction and 2kK/2kE correction as the serial path -- the issue #35 fix had not reached the GPU kernels, so large phases returned the wrong branch. - elliptic123: retire the stale private elliptic12i/elliptic12ic copies that shadowed the repaired public function and preserved old behavior; all dispatch now routes through matlab/src/elliptic12i.m. - elliptic12 at m=1: detect first-kind pole crossings from the original phase (reduction previously hid them: F(pi,1) returned 0) and restore the full quasi-period contribution to E. - ellipticBD small m: compute D via Carlson RD and B as RF - RD/3 -- (K-E)/m cancelled catastrophically (at m=1e-20: D=0, S wrong by 1e20). - elliptic3 near the third-kind pole: hybrid evaluation, quadrature on regular inputs, Carlson near n*sin^2(phi) = 1; validate pole crossings with abs(phi). - ellipj extreme inputs: no clipping of representable m, reduce by 2K before the Landen descent, reconstruct quasi-period signs; stable dn/sech forms. - elliptic12i: safe internal denominator at the exact m=0 endpoint. - uniquetol grouping: group exact duplicates only -- the 1e-11 tolerance silently substituted one m's result for another's (5e-12 apart -> 1.69e-9 discrepancy). Tests: testRegressionFollowup.m (6 blocks) pins each fix; testGpu extended for the GPU period reduction. Octave: 220/220 blocks across all 16 test files. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ncellation fixes Same audit round as the MATLAB commit (see docs/specs/post-0d09740-regression-audit.md): - Remove every runtime scipy import (theta, nome, inverse, complex functions): a base install failed only when those public functions were first called. Replaced with backend-native theta series, Carlson/AGM forms, and fixed-iteration inverse solvers. - JAX/PyTorch dispatch: replace value-dependent Python branches, NumPy casts, and scalar flattening with masked array expressions and fixed iteration counts so tracing works and results stay on device (complex, inverse, theta, jacobi EDJ, weierstrass, arclength). - elliptic12 m=1 poles, ellipticBD small-m cancellation, elliptic3 via Carlson RF/RJ throughout, ellipj clipping/large-phase reduction, Weierstrass lattice-point detection by ULP-scale tolerance (P(1e-11) is finite again) and explicit rejection of complex input, theta endpoint values (theta(1,v,0) = 0 exactly, no endpoint clipping), Carlson RJ safe under tracing, arclength_ellipse backend-native. Tests: test_regression_followup.py pins each fix; test_backends.py runs the dispatch matrix (NumPy/Torch/JAX). 451 passed + 1 optional-JAX skip in the default environment; 458 with Torch and JAX installed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…e install check Octave CI ran only 6 of 15 test files; discover every test*.m instead. The Torch/JAX jobs installed their backends but never passed those arrays through the library; they now run the dedicated dispatch tests. Add a base-dependency-only job proving the public API works without scipy, and gate publishing on all of the above. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The spec records the full audit: findings table with severity and data impact, numerical/differential verification against mpmath and Carlson references, automated-verification counts, and the deliberate limits (GPU kernels reviewed but not hardware-run; Weierstrass real-only; elliptic3 principal-value continuation not implemented). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An independent adversarial review of this branch (Codex, cross-checked against mpmath at 40-100 digits) produced five material counterexamples. Each is verified, fixed in the affected port(s), and pinned by a regression test (testEdgeCases.m block Q, TestAdversarialRound): - elliptic3 (MATLAB): negative amplitudes formed 0*Inf = NaN whenever the complete integral has a pole (Pi(-1|.5,1), Pi(-1|1,.2)) -- the complete-integral correction is now applied only where a half-period or reflection is actually present. - elliptic12i (both ports): the A&S 17.4.11 decomposition loses ~sqrt(eps/m) digits as m -> 0 -- F(0.2i|1e-20) returned ~0 instead of 0.2i, 9.2e-3 error at m = 1e-14. Added the Maclaurin series through m^2, switched at m*max(1, e^(2|psi|)) < 1e-4; crossover error ~2e-12 against 40-digit mpmath, and the series subsumes the m = 0 override. - weierstrassP/PPrime/Zeta (both ports): pole tolerances (abs(sn) < eps^(1/3) in MATLAB -- a 6e-6 window! -- and 8*eps*max(1,|z|) in python) replaced huge finite near-pole values with Inf: P(1e-16) is 1e32, not a pole (DLMF 23.9.2). A pole is now only the exact lattice point. - inversenomeq (both ports): replaced by the DLMF 20.9.1 closed form m = (theta2(0,q)/theta3(0,q))^4 with q^(1/4) kept outside the ratio. The python bisection had a 2^-64 absolute floor (m(1e-30) came back 2.7e-20 instead of 1.6e-29); the MATLAB interpolation tables were documented-unreliable outside [1e-5, 0.76]. Both are now exact at every scale and the computed q_max endpoint is accepted. - carlsonRC (python): the branch selector compared |y-x| against an ABSOLUTE 1e-14, sending every small-scale input down the degenerate x == y branch -- RC(1e-20,2e-20) was 27% off and contaminated RJ. Selection is now relative; DLMF 19.20 homogeneity (RF,RC ~ l^-1/2, RD,RJ ~ l^-3/2) is tested at lambda = 1e+/-20 in both ports. Also from the review: reversed arc intervals are signed for circles as for ellipses (python applications.py dropped its abs); DLMF citations corrected (19.25.14 for incomplete Pi, 19.25.5/19.25.9 for F/E); docs now state the double-precision phase-reduction bound for ellipj (full precision to |u| ~ 1e12, phase lost by ~1e16 -- shared by every double implementation) and the A&S 17.4.11 branch convention of elliptic12i (Re F = K(m) on Re u = pi/2, divergent as m -> 1, where mpmath/Mathematica may follow a different sheet). Claims from the review NOT adopted, with reasons: the m -> 1 complex values are a branch convention, not an error (now documented); the u = 1e16 Jacobi phase loss is the double-precision wall itself -- fixing it requires K(m) to ~32 digits, recorded as a deliberate limit. New anchors: ellipticBD at m = 0.2, 0.7, 0.999 (mpmath, dps=40) and elliptic12i on both sides of the series crossover. Octave 240/240 blocks across 16 files; pytest 458 + 1 optional-JAX skip. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…s (both ports) An independent adversarial pass fuzzed every public function of both ports against mpmath (dps=40) over parameter endpoints (m -> 0, m -> 1, m = 1), extreme scales, near-pole/near-lattice arguments, exact period multiples +/- ulps and complex arguments across the branch point. Each candidate was classified by evaluating scipy at the same DOUBLE inputs: wherever scipy reaches machine precision the loss is ours. Nine implementation defects, all fixed and pinned (testEdgeCases.m block R, test_edge_cases.py::TestAdversarialRound2): - Delta^2 = 1 - m sin^2 and 1 - n sin^2 formed by subtraction (python elliptic12; elliptic3 both ports, Carlson path and GL integrand): F(pi/2-1e-9 | 1-eps/2) off by 4e-3, Pi(pi/2-1e-6 | m, n=1) off by 3e-5. Now (1-m) + m cos^2 and (1-n) + n cos^2. - F(phi|1) = log(tan(pi/4+phi/2)) (both): F(0|1) = -1.1e-16, wrong sign at 1e-16. Now atanh(sin phi). - (ratio-1)/m in the A&S 17.4.11 decomposition (both): Im F(pi/2+1e-9 i) came back 0 and small m lost sqrt(eps/m) digits. Derived the cancellation-free form: with X1 = cot^2(phi) + Y, Y solves Y^2 + B'Y - C' = 0 and tan^2(mu) = 2 sinh^2 csc^2 / (B' + sqrt(B'^2+4C')) -- m cancels analytically, m = 0 is exact. Verified against the original at 120 digits (the 50-digit reference itself lost digits at cot^2 ~ 1e-33, a cautionary note in its own right). The m^2 series stays for m_eff < 1e-4 where it is still ~30x more accurate than the closed form. - Landen back-substitution asin(c sin/a) near +/-1 as m -> 1 (both, CPU and GPU): cn(9.4 | 1-eps/2) off 5e-10. Now atan2(c sin, sqrt(a^2 cos^2 + b^2 sin^2)) using a^2 - c^2 = b^2 for the AGM triple. - R_C for y << x used arctanh(sqrt(1-y/x)) (both): RC(3,1e-10) lost 8 digits and contaminated R_J and J. Now log1p(((x-y)/(sqrt x + sqrt y) + sqrt(x-y))/sqrt y)/sqrt(x-y) -- the intermediate plain-log form still lost 9 digits at RC(1+1e-13, 1), which is what tripped the RJ(x,y,z,z) == RD identity in testCarlson before the log1p. - R_J duplication cap 30 (both): the ratio exponent (base 4) shrinks by one per step, so ratios beyond ~1e16 never converged -- RJ(1e-20,2e-20,3e-20, .5) was 11% off. python: 60 fixed steps (documented limit ~3e32, JAX- traceable); MATLAB: adaptive break with cap 200 (also RF 20 -> 200, RD 30 -> 200; the break decides). - Two zero Carlson arguments (both): RF(0,0,1) returned 2e6; now Inf (DLMF 19.16). - inverselliptic2 (both): folding a tiny negative z through 2E1 - (z+2E1) lost its digits and the tol-gated Newton froze at an absolute 1e-12 (rel 1e-7 at z = -1e-9 E1). Oddness first; unconditional (python) / relative-stop (MATLAB) Newton. - 1 - m by subtraction (both ports' Weierstrass, MATLAB nomeq via ellipke(1-m), MATLAB weierstrassP/PPrime reducing inside ellipj): q(1e-16) 11% off, q(1e-17) = 0, P(2 omega1 + 1e-9) 40% off on a near-m=1 lattice. 1-m = (e1-e2)/(e1-e3); K' = R_F(0, m, 1) from the exact argument; reduce z by 2 omega1 before ellipj (sn^2 and cn dn/sn^3 are invariant under w -> w + 2K, no sign bookkeeping). Not defects, recorded as deliberate limits in the spec: F(phi|m) for m -> 1 at phi -> pi/2 and the Weierstrass functions within ~1e-9 omega1 of a lattice point now sit exactly at the input's conditioning floor (2 eps |z| / |z - 2k omega1|); the m -> 1 complex values follow the A&S 17.4.11 sheet. Reference-construction lesson written into both test files: anchors are evaluated at the exact double the library receives (mpf(float(x))) -- the decimal input a test author types and its double rounding differ, near singularities, at the 1e-9 level (F(pi/2-1e-9 | 1-eps/2): 19.6599302656 vs 19.6599302792). Octave 16/16 files (246 blocks); pytest 465 + 1 optional-JAX skip. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ckend)
Ran the library on a real NVIDIA L4 (Cloud Run job, europe-west1): PyTorch
CUDA, JAX CUDA, and Octave + ocl/OpenCL. Three things only hardware could
catch, plus one the run's warnings pointed at:
- gpu_ellipj (MATLAB): dn(I) = sqrt((1-m(I)) + m(I).*cn(I).^2) mixed the
column m(I) with the row-indexed cn(I) into a 4005x4005 outer product
("nonconformant arguments") -- the kernel had never executed on a device.
Keep cn as a column local. Every MATLAB GPU kernel is now also exercised
locally through gpuArray/gather identity stubs (worst CPU/GPU-path
difference 3.8e-16 across elliptic12, ellipj, elliptic3, ellipticBDJ,
jacobiThetaEta, weierstrassP/Zeta; matrix shapes preserved).
- gpu_jacobiThetaEta (MATLAB) still carried the retired AGM-product form
with its input perturbation hack, so GPU theta values disagreed with the
CPU q-series by up to 5e-9. Now the same q-series on gpuArray.
- carlsonRF/RJ (python): the two-zero-argument mask used .astype, which
torch tensors do not have -- every torch call into R_F/R_J raised. Pure
boolean algebra now; the whole torch backend was re-verified on 44 outputs.
- F(phi|1): atanh(sin phi) saturates when sin phi rounds to 1
(F(pi/2-1e-9|1) came back inf with a divide-by-zero warning in the run
log); asinh(tan phi) -- the inverse Gudermannian -- is exact at 0, odd,
and finite there (21.416412873533655 vs mpmath 21.416412873533656).
Both ports.
JAX CUDA on the L4: all 44 outputs on device and within 1e-14 of numpy
except Weierstrass P/P'/zeta at 1-3e-12 -- 1/sn^2 and theta1'/theta1 near
lattice points amplify the 1-2 ulp differences between GPU and CPU
transcendental implementations; a cross-backend, not a correctness, gap.
2e6-point jitted elliptic12 on the L4: 0.018 s.
Octave 16/16 files (246 blocks); pytest 465 + 1 optional-JAX skip.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Second hardware run (Cloud Run L4, europe-west1): JAX CUDA clean on all 44
outputs (jitted 2e6-point elliptic12: 0.018 s vs 7.30 s numpy); the pytest
backend matrix passed on the device; two further defects only a real
device exposes:
- weierstrassP/PPrime/Zeta/Sigma (python): _broadcast4 turned the
Python-float roots into 0-d CPU tensors, and torch's broadcast_arrays
refuses to mix them with CUDA tensors ("Expected all tensors to be on the
same device"). Scalars are now materialised with full_like on the
device of the array input. Located without a GPU via torch 'meta'
tensors, which carry a device but no data and still enforce the
same-device rule -- that probe now passes for every public function.
_agm.py's bare xp.ones(N) (currently unused) had the same hazard.
- gpu_ellipticBDJ (MATLAB): handed ocl arrays to the Carlson duplication,
which needs logical indexing ocl does not support, and carlsonRF's
isreal() rejects ocl arrays outright. There is no OpenCL kernel for it;
the GPU dispatch now evaluates the serial core on host arrays (identical
results), which is stated in the code.
Octave 16/16 files; pytest 465 + 1 optional-JAX skip; torch CPU smoke of
the full public API 44/44.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… core Third L4 run: torch CUDA 44/44, JAX CUDA 44/44, pytest backend matrix 10/10, every Octave OpenCL kernel within 4.8e-16 of the CPU, testGpu 9/9 on the device. The forced-GPU pass over testEdgeCases then hit the last sibling of the ellipticBDJ defect: gpu_ellipticBD wrapped m in gpuArray before calling the Carlson-based core, which carlsonRF rejects (isreal(ocl array) is false). Host arrays now, identical results. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Fifth L4 run's stack trace: ocl's gather() rejects non-ocl arrays (ocl_to_octave: invalid argument type), and the previous fix had left gather(B/D/S) after switching gpu_ellipticBD to the host core. Every other gather() in the tree applies to a real gpuArray result (ellipj, elliptic12, elliptic3, jacobiThetaEta), all verified on the device. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…s, NaN handling (both ports)
Dense random fuzz (m log-uniform to within 1e-14 of both endpoints, phases
over +/-13 periods and at odd multiples of pi/2 with 1e-9 jitter) against
mpmath, plus a logic/API abuse probe. Findings:
- Phase reduction u - k*pi rounded k*pi, costing eps*|u| in the reduced
phase; near pi/2 at m -> 1 that is amplified ~1e5x into Z and Pi (1e-10).
Cody-Waite split (u - k*PI_HI) - k*PI_LO in elliptic12, elliptic3 and
ellipticBDJ, both ports. What remains is the double-precision floor of
the reduced phase itself.
- python ellipj returned sn(u | 0.5) -- the interior AGM placeholder -- for
m = 1.5 or m = NaN, with no error. New check_range() raises on eager
numpy and NaN-masks on device backends; applied to ellipj, elliptic12 and
nomeq; NaN inputs now propagate.
- Every `xp is np` eager check in the package (elliptic12i, elliptic3,
inverse, carlson, nome, applications) was silently skipped for ndarray
and numpy-scalar inputs, because array_namespace() returns the
array_api_compat.numpy module for those. is_numpy() recognises both.
- MATLAB elliptic12 crashed inside unique() for a NaN m ("subscripts must
be ..."); NaN in, NaN out now, without leaking into neighbours.
- elliptic12i(-0.0) and (0) returned eps from the cot(phi) nudge; exact
zero, sign preserved.
- python R_J: 60 -> 100 duplications (argument-ratio limit 3e32 -> 4e56);
the fuzz hit ratio 1.9e44 at 9e-12.
Verified clean in the same round: no input mutation, vector/scalar
agreement bit-for-bit over 300 random points incl. m in {0, 1, 1e-17,
nextafter(1,0)}, matrix shapes preserved by every function, empty in ->
empty out, F(conj u) = conj F(u), F(-u) = -F(u), sn(conj u) = conj sn(u).
Tests: testEdgeCases.m block S, TestAdversarialRound3 (anchors at the
exact double inputs). Octave 16/16 files; pytest 468 + 1 optional-JAX skip.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…v (both ports) Adversarial round 4 (loop iteration 2): complex-plane periods of the Jacobi functions (2K, 2iK', 4iK'), Pi with n > 1 below its pole and n = -1e4, degenerate/near-degenerate Weierstrass lattices and extreme-scale roots, degenerate ellipses -- all clean. Two real findings: - Theta series lost accuracy at large arguments: forming (2n+1)*v and 2n*v as double products rounds by eps*|k v| (python 1.4e-12 at v ~ 1.2e8, 9e-9 at v ~ 1e11; MATLAB 2.2e-10 because theta() also round-tripped v -> u = 2Kv/pi -> v through jacobiThetaEta). All theta series in both ports (theta, theta_prime, jacobiThetaEta, the Weierstrass theta1 loops) now take sin/cos of the multiples from the angle-addition recurrence started at sin v, cos v; MATLAB gains a shared theta_series.m and theta() evaluates on v directly. Verified exact (1e-16) to v = 1e11 against mpmath. - python Weierstrass accepted unsorted roots and returned NaN silently (MATLAB errors); now raises on numpy. Equal neighbouring roots stay accepted as the legitimate m = 0 / m = 1 degenerate lattices. Tests: testEdgeCases.m block T, TestAdversarialRound4. Octave 16/16 files; pytest 471 + 1 optional-JAX skip. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… phase tail
Adversarial rounds 4-5 (see docs/specs/post-0d09740-regression-audit.md):
* carlsonRF/RD/RJ: per-element convergence. The duplication loop stopped
on a whole-vector test, so the number of steps applied to an element
depended on its batch mates: chunked vs serial elliptic3 differed by an
ulp on 6 of 1000 points. An `active` mask now freezes each element at
its own converged step.
* elliptic12 / ellipj (CPU and GPU): use the first converged AGM row per
element, a(n+1), instead of the batch-wide last row a(max(n)); K and the
Landen back-substitution scale no longer depend on batch composition.
* elliptic12 GPU: Cody-Waite tail term in the phase reduction, as in the
CPU path. Without it Z was off by 3.6e-11 and E by 1.2e-10 at u = 1e6,
m = 1 - eps/2 (k * 1.22e-16 with k = 318310).
* par_worker: recursion guard. A worker that sees parallel mode still
enabled re-enters the parallel dispatcher when N is an exact multiple of
chunk_size and recurses until SIGILL; the worker now forces
elliptic_config('parallel', false) for the duration of its call.
* testParallel: new block evaluates every parallel-capable function
serially and chunked (N in {cs-1, cs, 2cs, 3cs, 3cs+7}) under temporary
get_nworkers/parcellfun stubs and requires bit-identical results.
Octave: 16 test files, 0 failures. Python: unchanged, pytest green.
GPU identity-stub probe: CPU vs GPU path max diff 0.0 on 4005 points.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…me m Both ports were evaluated on identical random doubles (m down to 1e-16 and up to 1-1e-16, |u| up to 1e5) and every disagreement above 1e-13 was adjudicated against mpmath at the exact inputs. See round 6 in docs/specs/post-0d09740-regression-audit.md. MATLAB * elliptic12 (CPU+GPU): m in [eps^2, ~5e-16] returned F = E = Inf -- the AGM converges in one step, no Landen step ran and the scale e stayed 0. Scale 2^(n-2) now in closed form. * elliptic12 (CPU+GPU): the E/K sum stopped one AGM term early; E was off by 1.4e-13 near m -> 1. * theta, theta_prime, jacobiThetaEta: nome from ellipke(1-m) rounded 1-m first; theta1 was off by 1e-5 at m ~ 1e-16. K'(m) = R_F(0, m, 1) now. * elliptic3: the reflection used elliptic3(double(pi/2)) as the complete integral; cos(double(pi/2)) = 6e-17 leaves a 2e-7 sliver at m = 1-eps/2 (3e-10 relative). Exact Carlson complete form now. * elliptic3: c < 0 accepted (the Python port already did), via Carlson. * carlsonRJ: series term E3 = XYZ + 2 E2 P + 4 P^3 (DLMF 19.36.2), was 3 P^3 -- O(eps^4) residual of 1e-13 relative at the 0.0015 tolerance. * ellipticBDJ: n > 1 beyond the pole returned complex J silently; n = 1 returned NaN (0*Inf in the period term). Error / masked now. Both ports * k*pi reduction: the "Cody-Waite" split used double(pi) as the head, so k*pi already rounded by eps*|u| (2.3e-10 at u = 1e6). New sub_kpi / _xputils.sub_kpi: 25-bit PI_A, PI_B (k*PI exact for k < 2^28) + PI_C, verified to 4e-16 over 20000 random k < 2^27. Jacobi Zeta at u = 8e4 went from 1.5e-11 to 3e-17. * carlson.py: same E3 coefficient (masked there by 100 fixed duplications). * ellipticBDJ.py: n > 1 beyond the pole raised nothing and returned 1.147 (principal value 0.859) because the private _rj_xp skips the p > 0 check; n = 1 gave NaN. ValueError (NumPy) / NaN (traced) now. Tests * testEdgeCases block U / TestAdversarialRound6: mpmath anchors at the exact doubles for every defect above. * testDocExamples.m runs every docstring Example block; pytest now runs --doctest-modules (one example printed a 0-d array). After the fixes the sweep agrees with mpmath to < 2e-15 relative on F, E, Z, Pi, theta1 in both ports. Octave: 17 files, 0 failures; pytest 478 passed.
cel evaluated through m = 1 - kc^2, which loses kc entirely below ~1e-8: cel1(1e-9) was Inf (MATLAB) / 2e6 (Python) against ln(4/kc) = 22.1. MATLAB also rejected kc > 1 (m < 0) and both ports returned Inf for p < 0, where Bulirsch defines the Cauchy principal value. Both cores are now Bulirsch's algorithm (Numer. Math. 13 (1969) 305): any real kc (the integral depends on kc^2 only), p < 0 gives the principal value (= Re Pi(1-p | m), checked against mpmath), p = 0 gives Inf, kc = 0 gives +-Inf unless b = 0 (finite limit). MATLAB iterates with a per-element active mask; Python runs a fixed 40-step ascent with converged elements frozen (backend-native). Cross-port results are bit-identical on 400 random points spanning kc = 1e-300 .. 100 and p = +-1e-3 .. 1e6, within 1e-16 of mpmath / the identity cel(1,0) + p cel(0,1) = K. Tests: testEdgeCases block V, TestAdversarialRound6.test_bulirsch_cel_is_kc_native. Docs: round 6.9 in docs/specs/post-0d09740-regression-audit.md. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…rclength) Sweep of the remaining outputs (complex E/Z, jacobiEDJ, ellipticBD, agm, arclength_ellipse) on 1000 random points against mpmath; rounds 6.10-6.13 in docs/specs/post-0d09740-regression-audit.md. * elliptic12i.py (Jacobi Zeta output): the complete K, E were taken as F(double(pi/2)|m); cos(double(pi/2)) = 6e-17 leaves K 5.8e-9 relative short at m = 1-eps/2 and Z inherited 1.8e-11. Exact Carlson complete forms now (3.5e-16). * ellipticBDJ (both ports): Delta^2 = 1 - m sin^2 cancels near phi = pi/2 as m -> 1 (relative 2.5e-9 at m = 1-1e-8), which R_D turned into 3e-10 in D(phi|m). Delta^2 = (1-m) + m cos^2 now. * jacobiEDJ (both ports): took am(u) at |u| ~ 1e3 before reducing, where the map phi -> D is steep (1/sqrt(1-m)); D_u(1520|1-1e-8) was off by 3e-10. Reduce u by 2K first, amplitude of the reduced argument, add 2k complete integrals; now at the eps*|u| floor in both ports. * arclength_ellipse.m: if(a<b)/elseif(a>b) on arrays has all-elements semantics, so any mixed array fell through to the circle formula for every element. Elementwise masks after broadcasting scalars. Tests: testEdgeCases block W, TestAdversarialRound6.test_round6c_*. Octave: 17 files, 0 failures; pytest green. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…q near q=1, nome2m Rounds 6.14-6.16 in docs/specs/post-0d09740-regression-audit.md. * elliptic12i (both ports): the period term pi*ceil(phi/pi - 0.5 + eps) (Python: + 1e-14) was counted from a separately rounded quantity, so for phi a few ulps (Python: 3e-14) below pi/2 it added a period that the sign term (-1)^floor(2phi/pi) had not crossed: Re F came out 3K instead of K. The period is now pi*ceil(k/2) from the same k = floor(2phi/pi). asin(sqrt(3)) is 2 ulps below pi/2, so elliptic123 returned K(3) = 3.003 (mpmath: 1.001) -- how the defect was found. * elliptic123.m, complete m > 1: evaluated elliptic12i exactly on the branch point of F(.|1/m), where the decomposition is sqrt(eps)-conditioned (1e-8 even after the fix). DLMF 19.7.3 closed forms now, also for elliptic123(pi/2, m); K(3), E(3), K(5), E(5) match mpmath to 2e-16. * inversenomeq (both ports): above q_max = 0.7789534 the 30-term series is not converged; MATLAB returned m > 1 (1.034 at q = 0.999), Python raised. The true 1-m is below eps/2 there, so both return exactly 1 (clamped <= 1 below q_max). * nome2m.m: captured its whole input array in the fzero objective (errored on any array) and the bracket covered only q < 0.62. Now an alias of inversenomeq. Tests: testEdgeCases blocks X, Y; TestAdversarialRound6 additions; the Python q_max test now asserts the rounded value. Octave 17 files, 0 failures; pytest 481 passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…ce stub Rounds 6.17-6.18 in docs/specs/post-0d09740-regression-audit.md. Shape sweep (matrix / column / row / mixed scalar-array calls of every function against scalar loops; the Python port passed, pinned by TestInputShapes): * ellipj: re-read cn(I) from its column-shaped output against the row m(I): 6x6 broadcast error for any column u (and hence jacobiEDJ). * jacobiThetaEta: returned a row for a matrix input. * inverselliptic2: the vector-wide Newton stop made values depend on the batch by an ulp; per-element mask now. * elliptic123: failed on any matrix; inputs flattened to rows (its legacy kernels preallocate rows), shape restored. GPU: * elliptic12 GPU path: K_per = 2 .* gpuArray(k_per) .* K_vals mixed a device array with a host matrix (K_vals became host in the per-element AGM change); ocl refuses that and the Octave section of the L4 run failed. Host product, then gpuArray. * tests/gpu_stub/gpuArray.m: strict stand-in for an ocl device array that errors on device-by-host-matrix operators and on logical indexing, plus testGpuStrict.m running every GPU path under it against the CPU path (agreement <= 2.8e-16). The identity stubs used before let three host/device defects through to hardware. L4 hardware run (image from d388368): torch and JAX 88/88 OK on device (max rel diff 3.6e-14, x50 over NumPy), backend pytest 10/10; the Octave OpenCL section failed on the K_per mixing fixed here. Octave: 18 files, 0 failures; pytest 483 passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…TLAB
Round 6.19 in docs/specs/post-0d09740-regression-audit.md. A sweep of
[] / NaN / +-Inf / [0.3 NaN 0.7] through every function:
MATLAB
* empty input -> empty output of the same shape in ellipticBDJ,
theta_prime, cel, weierstrassP/Zeta/Sigma/PPrime, carlsonRF/RD/RJ/RC,
arclength_ellipse, elliptic123 (they rejected [] against a scalar).
* nomeq aborted inside ellipke ("algorithm did not converge") on a single
NaN element; inversenomeq rejected NaN as out of [0,1); both now isolate
NaN and keep the domain error for finite out-of-range values.
* elliptic12i raised "Input arguments must be real" for NaN or Inf phases
because (-1)^NaN is complex NaN in Octave.
* Carlson wrappers: NaN, Inf and negative arguments give real NaN (were
complex NaN); carlsonRJ errors for p <= 0 like the Python port instead of
returning complex garbage from sqrt of a negative.
Python
* cel(NaN, ...) returned pi/2: a NaN kc never became active in the ascent.
Tests: testEdgeCases block AA; TestEmptyNaNInf. Octave 18 files, 0
failures; pytest 484 passed.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…portable par_worker Round 6.20 in docs/specs/post-0d09740-regression-audit.md. * L4 run (image from 5b3c561): torch/JAX 88/88 OK on device, backend pytest 10/10, Octave OpenCL parity <= 4.8e-16 (x2.7) -- and the NaN block of testEdgeCases failed under the GPU flag: a NaN m fell through find(m ~= 1 & m ~= 0) in gpu_elliptic12 / gpu_ellipj, the AGM loop exited at once and F came back equal to u. Reproduced locally with the strict device stub (testGpuStrict now carries NaN cases); NaN masks added. * Octave's ellipke aborts ("algorithm did not converge") when any element is NaN, taking theta, theta_prime, jacobiThetaEta, elliptic12i and inverselliptic2 down with it. New ellipke_safe.m (NaN-propagating) at the six call sites; NaN m gives NaN after the q = 0 stand-in. * par_worker.m: try/catch instead of the Octave-only unwind_protect, so the file also parses in MATLAB. The whole Octave suite also passes under --traditional (MATLAB-compatibility mode). * complex_elliptic.py: silence a masked 0/0 in the complex E assembly. * Audit doc: "Automated verification" refreshed for the state after round 6. Octave: 18 files, 0 failures; pytest 484 passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…, elliptic123 Round 6.21: elliptic3 expanded c from the still-scalar u before u was expanded from m and rejected elliptic3(0.3, [0.2 0.5 0.9], 0.3) as "must be the same size"; theta preallocated its output before broadcasting and hit a size error for theta(1, 0.3, [..]); elliptic123 restored the shape of the scalar phase. Broadcast to the largest input first, preallocate after. Empty inputs give empty outputs in elliptic3 and theta as elsewhere. Tests: testEdgeCases block AB. Octave 18 files, 0 failures. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…he CPU core Round 6.22: the OpenCL branch of elliptic3 was the 20-node Gauss-Legendre rule only. The serial core switches to the Carlson form (DLMF 19.25.14) where the rule loses digits -- endpoint denominators below 0.25, and c < 0 since round 6.8 -- but the GPU path had no such fallback: on the L4, Pi(1|0.5,-100) came back 3.8e-9 off and Pi(4|0.9,-100) 5.4e-10, failing the new block U anchors under the GPU flag. The strict device stub reproduces it locally. The Carlson block is now the shared subfunction elliptic3_carlson; the GPU path evaluates the "danger" subset (and NaN elements) on the host with it and sends only the regular elements to the kernel. GPU and CPU paths are bit-identical on the c < 0 / near-pole cases (testGpuStrict extended). Octave: 18 files, 0 failures; pytest green. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…ance in the shape test
* testGpuStrict.m / testDocExamples.m: mfilename('fullpath') is empty inside
test blocks when test() is called with a full path from the repo root
(CircleCI), so the relative '../src' and 'gpu_stub' lookups failed
(1/2 and a 0-example run). Locate the directories from
which('elliptic12') instead. The suite now passes with CircleCI's exact
invocation from the repo root.
* TestInputShapes: batch-vs-scalar tolerance 4e-15 instead of 4e-16 --
NumPy's SIMD sin/cos on Linux differ from the scalar libm path by an ulp
(ubuntu 3.13 and the torch job failed on 1.8e-15).
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Follow-up audit on baseline
0d09740(the issue #35 fix round), plus an independent adversarial review of the audit itself. The initial patch did not fully close the regression class: the #35 phase-reduction fix had not reached the MATLAB GPU kernels, and stale private copies insideelliptic123.mcould still bypass the repaired publicelliptic12i. A subsequent adversarial review (Codex, cross-checked against mpmath at 40–100 digits) of these very fixes produced five further counterexamples, all now fixed and regression-tested.Full findings tables with severity and data impact:
docs/specs/post-0d09740-regression-audit.mdAudit round (commits 1–4)
elliptic12/ellipj: apply the same quasi-period reduction as the serial path (the weird branching with elliptic12i (Matlab/Octave) #35 fix was CPU-only).elliptic123: retire privateelliptic12i/elliptic12iccopies that preserved pre-fix behavior.ellipticBDsmall-m cancellation:(K-E)/m→ CarlsonRD/RFforms (atm=1e-20,Swas wrong by ~1e20).elliptic12atm=1: pole crossings detected from the original phase (F(pi,1)returned 0).Adversarial-review round (commit 5)
elliptic3(-1, .5, 1)→ NaN (0·Inf against a complete-integral pole)F(0.2i|1e-20)→ 0; complex F/E lose √(eps/m) digits as m→0 (both ports)weierstrassP(1e-16)→ Inf though the true value is the finite 1e32 (both ports)inversenomeq(1e-30)off by 9 orders (Python); MATLAB tables unreliable outside [1e-5, 0.76]m = (θ₂/θ₃)⁴, exact at every scalecarlsonRC(1e-20, 2e-20)27% off — absolute branch tolerance broke DLMF 19.20 homogeneity (Python)Review claims not adopted, with documented reasons: the m→1 complex values are the A&S 17.4.11 branch convention (now documented in both ports), and the
u=1e16Jacobi phase loss is the double-precision reduction wall shared by every double implementation (recorded as a deliberate limit).Verification
Adversarial-review round 2 (commit 6)
A second, independent adversarial pass fuzzed every public function of both ports against mpmath (dps=40) over parameter endpoints, extreme scales, near-pole/near-lattice arguments, exact period multiples ± ulps and complex arguments across the branch point, classifying each candidate by whether scipy reaches machine precision at the same double inputs. Nine implementation defects — all cancellation or convergence, all in both ports unless noted — fixed and pinned (
testEdgeCases.mblock R,TestAdversarialRound2):1−m·sin²φ,1−n·sin²φby subtraction (elliptic12py,elliptic3)(1−m)+m·cos²φF(φ|1)=log(tan(π/4+φ/2))F(0|1)=−1.1e-16atanh(sin φ)(ratio−1)/min A&S 17.4.11Im F(π/2+1e-9i)=0; √(eps/m) losstan²μ, verified at 120 digitsasinnear ±1 as m→1 (CPU+GPU)cn(9.4|1−eps/2)off 5e-10atan2form viaa²−c²=b²R_Carctanh branch, then plain logRC(3,1e-10)off 3e-9 → contaminatedR_Jlog1pformR_Jduplication cap 30RF(0,0,1)=2e6Infz=−1e-9·E11−mfrom roots /ellipke(1−m)/ reduction insideellipjq(1e-16)11% off;P(2ω₁+1e-9)40% off on near-m=1 lattice1−m=(e1−e2)/(e1−e3),K'=R_F(0,m,1), reduce by 2ω₁ firstEverything remaining in the fuzz sits at the input's conditioning floor (documented as a deliberate limit). Anchors are evaluated at the exact double the library receives — the decimal input and its rounding differ at 1e-9 near singularities, which is now noted in both test files.
Boards after round 2: Octave 16/16 files (246 blocks), pytest 465 + 1 optional-JAX skip.
Hardware verification on an NVIDIA L4 (commits 7–10)
Six Cloud Run GPU job executions (europe-west1, 1× L4, Octave 7.3 +
ocl1.2.4 over OpenCL 3.0; PyTorch 2.6+cu124; JAX 0.11 CUDA 12) — total cost ≈ $0.37. Four defects only real hardware could expose, all fixed:dn(I) = … m(I).*cn(I).^2mixed a column with a row-indexed vectorgpu_ellipj(MATLAB)gpu_jacobiThetaEta(MATLAB).astypeon torch tensors; Python-float roots became 0-d CPU tensorscarlsonRF/RJ,weierstrass._broadcast4(Python)R_F/R_Jraised; Weierstrass "cuda:0 and cpu" device mixoclarrays handed to the Carlson duplication /gather()on host arraysgpu_ellipticBDJ,gpu_ellipticBD(MATLAB)carlsonRF: must be real,ocl_to_octave: invalid argument typeFinal run (all green): torch CUDA 44/44 public outputs on device and within 1e-12 of numpy; JAX CUDA 44/44; pytest backend matrix 10/10 on the device; every Octave OpenCL kernel (
elliptic12,ellipj,elliptic3,ellipticBDJ) within 4.8e-16 of the CPU path;testGpu9/9 and the fulltestEdgeCases39/39 executed withelliptic_config('gpu', true)on real hardware. Throughput on 2·10⁶ points: JAX jittedelliptic120.018 s vs 7.7 s numpy (≈×400); torch CUDAelliptic12+ellipj0.33 s vs 12.7 s (≈×39); Octave OpenCL 3.6 s vs 15.0 s (≈×4.2).Two local pre-flight tools came out of this and are worth keeping in mind: running every MATLAB GPU kernel through
gpuArray/gatheridentity stubs (catches shape bugs), and torchmetatensors, which carry a device with no data and still enforce same-device rules (catches CPU/CUDA mixing without a GPU).Deliberate limits
GPU kernels now hardware-verified on an L4 (see above); MATLAB parallel needs a real multi-worker run before release claims; Weierstrass functions reject complex input explicitly;
elliptic3does not implement principal-value continuation;ellipjphase reduction is double-precision (full accuracy to |u|~1e12).🤖 Generated with Claude Code
Bug-hunt rounds 4–6 (commits 43884a6 … dfb10b0)
Self-paced adversarial loop after the hardware round: parallel-chunking path, batch independence, cross-port random sweeps at extreme parameters (m down to 1e-16 and up to 1−1e-16, |u| to 1e5) adjudicated by mpmath at the exact doubles, input-shape and empty/NaN/Inf sweeps, docstring examples as tests, and a strict device stub for the GPU kernels. Full tables: rounds 4–6 in
docs/specs/post-0d09740-regression-audit.md.Highest-impact fixes (both ports unless noted):
k·πreduction useddouble(pi)as the "Cody-Waite" headsub_kpi: 25-bit split,k·πexact for k < 2²⁸ (Z now 3e-17)celevaluated throughm = 1 − kc²cel1(1e-9)= Inf / 2e6 (true 22.1); kc > 1 rejected; p < 0 gave InfcarlsonRJseries termE3with3P³(DLMF:4P³)R_Jand everything built on itelliptic12iperiod termceil(φ/π − 0.5 + eps)Re F = 3Kinstead ofK;elliptic123(3)returnedK(3) = 3.003(true 1.001)floor(2φ/π)that picks the sign;elliptic123complete m > 1 via DLMF 19.7.3 closed formselliptic12.mfor m ∈ [eps², 5e-16]F = E = Inf(AGM converges in one step, Landen scale never set)2^(n−2)ellipke(1−m)(MATLAB)K′ = R_F(0, m, 1)elliptic3.mreflection usedΠ(double(π/2))as the complete integralellipticBDJΔ² = 1 − m sin²φ;jacobiEDJtookam(u)before reducingellipticBDJn > 1 beyond the poleelliptic3; period term maskedarclength_ellipse.mon arraysinverselliptic2Newton stop)testParallelrequires bit-identical resultselliptic3kernel without the Carlson fallback (3.8e-9 for c = −100)tests/gpu_stubstrict stub +testGpuStrict.mreproduce all three classes on a laptopinversenomeqnear q → 1[]; a NaN abortedellipke-based functionsellipke_safeDeliberate limits documented (conditioning, not defects): complex sn/cn/dn near the poles u = iK′ and Weierstrass functions near lattice points (
eps·K′/|u − iK′|),elliptic12iexactly at its branch point (√eps),theta/jacobiThetaEtaat |v| ~ 1e3–1e4 (eps·|v|),ellipjat |u| ~ 1e5 (eps·|u|, the period 4K is not a constant).Verification after round 6
--traditional; newtestDocExamples.m,testGpuStrict.m, chunking block intestParallel.m.testGpu9/9 andtestEdgeCases50/50 under the GPU flag, OpenCL parity ≤ 4.8e-16 (×3.1).