Post-0d09740 regression audit: GPU period fix, elliptic123 routing, scipy-free runtime, adversarial-review round

.circleci/config.yml

@@ -33,7 +33,7 @@ jobs:
       - run:
           name: Run Octave tests
           command: |
-            octave --eval "addpath(fullfile(pwd, 'matlab', 'src')); cd matlab/tests; test testElliptic12; test testElliptic3; test testEllipj; test testThetaPrime; test testAgm; test testJacobiThetaEta;"
+            octave --eval "addpath(fullfile(pwd, 'matlab', 'src')); files=dir(fullfile(pwd, 'matlab', 'tests', 'test*.m')); failed=false; for k=1:numel(files); [n,nmax]=test(fullfile(files(k).folder,files(k).name),'quiet'); fprintf('%s: %d/%d\n',files(k).name,n,nmax); failed=failed || n~=nmax; end; if failed; exit(1); end"
 workflows:
   build:
     jobs:

.github/workflows/python.yml

@@ -2,9 +2,9 @@ name: tests
 
 on:
   push:
-    paths: ["python/**"]
+    paths: ["python/**", ".github/workflows/python.yml"]
   pull_request:
-    paths: ["python/**"]
+    paths: ["python/**", ".github/workflows/python.yml"]
   release:
     types: [published]
 
@@ -13,6 +13,32 @@ defaults:
     working-directory: python
 
 jobs:
+  test-standalone:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+      - name: Install runtime dependencies only
+        run: pip install -e .
+      - name: Exercise public functions without SciPy
+        run: |
+          python - <<'PY'
+          import importlib.util
+          import elliptic
+
+          assert importlib.util.find_spec("scipy") is None
+          elliptic.theta(1, 0.2, 0.5)
+          elliptic.theta_prime(1, 0.2, 0.5)
+          elliptic.jacobiThetaEta(0.2, 0.5)
+          elliptic.nomeq(0.5)
+          elliptic.inversenomeq(0.04)
+          elliptic.inverselliptic2(0.5, 0.5)
+          elliptic.elliptic12i(0.7 + 0.2j, 0.5)
+          elliptic.arclength_ellipse([2.0, 3.0], [3.0, 2.0])
+          PY
+
   test:
     runs-on: ${{ matrix.os }}
     strategy:
@@ -62,7 +88,7 @@ jobs:
   publish:
     name: Build and publish to PyPI
     runs-on: ubuntu-latest
-    needs: [test, test-torch, test-jax]
+    needs: [test-standalone, test, test-torch, test-jax]
     if: github.event_name == 'release' && github.event.action == 'published'
     environment: pypi
     permissions:

docs/GPU.md

@@ -6,6 +6,10 @@ configuration flag.  The same source code paths work in both **MATLAB**
 (Parallel Computing Toolbox / CUDA) and **Octave** (ocl Forge package /
 OpenCL).
 
+For `elliptic3`, regular inputs run on the GPU. Inputs close to an endpoint
+pole are intentionally gathered and evaluated by the Carlson CPU path because
+fixed GPU quadrature is not accurate enough there.
+
 ---
 
 ## Benchmark results
@@ -20,9 +24,9 @@ Octave 6.4.0 · ocl 1.2.4 · CUDA driver 535, OpenCL 3.0.
 | `elliptic3` | 4 M | 2.668 s | 1.167 s (2.3×) | 0.199 s | **13.4×** |
 | `jacobiThetaEta` | 1 M | 1.460 s | 0.583 s (2.5×) | 0.494 s | **3.1×** |
 
-`elliptic3` achieves 13× because it is a pure Gauss-Legendre quadrature
-(no AGM, no sequential dependencies) — every element is completely
-independent and maps trivially to GPU threads.
+These historical `elliptic3` measurements use regular inputs, for which the
+Gauss-Legendre GPU path remains active. Every element is independent and maps
+directly to GPU threads.
 
 ### Hardware utilisation at N = 1 M
 

docs/specs/codebase-weakness-audit.md

@@ -1,5 +1,10 @@
 # Codebase Weakness Audit — Python Elliptic Library
 
+> Historical snapshot. The high-severity items and the numerical defects
+> confirmed from this list were addressed in the
+> [post-0d09740 follow-up audit](post-0d09740-regression-audit.md). Keep this
+> document as the pre-fix evidence, not as the current defect list.
+
 **Date:** 2026-04-21
 **Scope:** `/home/igor/Work/elliptic/python/elliptic/`
 **Purpose:** Identify weak or breaking points (special values, complex numbers, inf/NaN, memory, exception handling, backend tracing) to guide future implementation work.

docs/wiki/Elliptic-Integrals.md

@@ -160,7 +160,7 @@ The library provides the AGM function directly:
 
 ## Carlson's Method
 
-The conventional methods for computing elliptic integrals are Gauss and Landen transformations, which converge quadratically and work well for elliptic integrals of the first and second kinds. Unfortunately they suffer from loss of significant digits for the third kind. Carlson's algorithm provides a unified method for all three kinds with satisfactory precision. The third kind integral in this library uses a Gauss-Legendre 10-point quadrature instead.
+The conventional methods for computing elliptic integrals are Gauss and Landen transformations, which converge quadratically and work well for elliptic integrals of the first and second kinds. Unfortunately they suffer from loss of significant digits for the third kind. Carlson's algorithm provides a unified method for all three kinds with satisfactory precision. Python evaluates the third kind directly with Carlson RF/RJ forms. MATLAB/Octave uses a hybrid: vectorised 20-node Gauss-Legendre quadrature on regular inputs and Carlson RF/RJ near endpoint poles, where fixed quadrature loses precision.
 
 ---
 

docs/wiki/elliptic.md

@@ -148,7 +148,7 @@ _See also_ `ELLIPKE`, `ELLIPJ`, `ELLIPTIC3`, `THETA`.
 
 [ELLIPTIC3](https://github.com/moiseevigor/elliptic/blob/master/src/elliptic3.m) evaluates incomplete elliptic integral of the third kind `Pi = ELLIPTIC3(U,M,C)` where `U` is a phase in radians, `0 < M < 1` is the module and `0 < C < 1` is a parameter.
 
-`ELLIPTIC3` uses Gauss-Legendre 10 points quadrature template described in [3] to determine the value of the Incomplete Elliptic Integral of the Third Kind (see [1, 2]).
+`ELLIPTIC3` uses vectorised 20-node Gauss-Legendre quadrature on regular inputs and switches to Carlson RF/RJ symmetric forms near endpoint poles. This preserves the fast path while avoiding fixed-quadrature precision loss as `M` or `C` approaches one.
 
 **General definition:**
 ```

matlab/src/arclength_ellipse.m

@@ -86,6 +86,17 @@
 %     Moiseev Igor
 
 %arguments
+% Empty input -> empty output of the same shape (elementwise semantics; the
+% size checks below would otherwise reject [] against a scalar).
+if nargin >= 2 && (isempty(a) || isempty(b) || (nargin == 4 && (isempty(theta0) || isempty(theta1))))
+    sz = size(a);
+    if isempty(b), sz = size(b); end
+    if nargin == 4 && isempty(theta0), sz = size(theta0); end
+    if nargin == 4 && isempty(theta1), sz = size(theta1); end
+    arclength = zeros(sz);
+    return;
+end
+
 if nargin ~= 2 && nargin ~= 4,
  error('ARCLENGTH_ELLIPSE: Requires two or four inputs.')
  return
@@ -96,23 +107,33 @@
  theta1 = 2*pi;
 end
 
+% Broadcast scalars to a common size, then branch ELEMENTWISE.  The previous
+% if(a<b) / elseif(a>b) on arrays used all-elements semantics, so any mixed
+% array fell through to the circle formula for every element.
+sz = size(a);
+for x = {b, theta0, theta1}
+    if numel(x{1}) > 1, sz = size(x{1}); end
+end
+a = a + zeros(sz);  b = b + zeros(sz);  theta0 = theta0 + zeros(sz);  theta1 = theta1 + zeros(sz);
+
 % Default solution for a==b (circles)
 arclength = a.*(theta1-theta0);
 
-% Ellipses (a<b or a>b)
-if(a<b)
-    % Theta measured from a axis = semi-MINOR axis
-    % Use standard formulation for E(phi,m)
-    [F1, E1] = elliptic12( theta1, 1 - (a./b).^2 );
-    [F0, E0] = elliptic12( theta0, 1 - (a./b).^2 );
-    arclength = b.*(E1 - E0);
-elseif(a>b)
-    % Theta measured from a axis = semi-MAJOR axis
-    % Standard formulation will not work ((1-(a/b)^2) < 0); instead use PI/2 - phi and b/a instead of a/b
-    [F1, E1] = elliptic12( pi/2 - theta1, 1 - (b./a).^2 );
-    [F0, E0] = elliptic12( pi/2 - theta0, 1 - (b./a).^2 );
-    % d(PI/2 - phi)/dphi = -1, so reverse operands in this difference to flip sign:
-    arclength = a.*(E0 - E1);
+% Ellipses: theta measured from the a axis
+lt = a < b;                   % a is the semi-MINOR axis: standard E(phi|m)
+if any(lt(:))
+    m = 1 - (a(lt)./b(lt)).^2;
+    [~, E1] = elliptic12(theta1(lt), m);
+    [~, E0] = elliptic12(theta0(lt), m);
+    arclength(lt) = b(lt).*(E1 - E0);
+end
+gt = a > b;                   % a is the semi-MAJOR axis: (1-(a/b)^2) < 0, use pi/2 - phi and b/a
+if any(gt(:))
+    m = 1 - (b(gt)./a(gt)).^2;
+    [~, E1] = elliptic12(pi/2 - theta1(gt), m);
+    [~, E0] = elliptic12(pi/2 - theta0(gt), m);
+    % d(pi/2 - phi)/dphi = -1, so reverse the operands to flip the sign:
+    arclength(gt) = a(gt).*(E0 - E1);
 end
 
 return;

matlab/src/carlsonRC.m

@@ -29,6 +29,15 @@
 %   [2] B.C. Carlson, "Computing Elliptic Integrals by Duplication,"
 %       Numer. Math. 33 (1979), 1–16.
 
+% Empty input -> empty output of the same shape (elementwise semantics; the
+% size checks below would otherwise reject [] against a scalar).
+if nargin >= 2 && (isempty(x) || isempty(y))
+    sz = size(x);
+    if isempty(y), sz = size(y); end
+    RC = zeros(sz);
+    return;
+end
+
 if nargin < 2, error('carlsonRC: requires two arguments (x, y).'); end
 if ~isreal(x) || ~isreal(y)
     error('carlsonRC: all input arguments must be real.');
@@ -38,7 +47,11 @@
 origSize = size(x);
 x = x(:).';  y = y(:).';
 
+% NaN, Inf and negative x give NaN (y < 0 is the principal value, handled by the core).
+bad = ~(x >= 0) | isinf(x) | isnan(y) | isinf(y);
+x(bad) = 1;  y(bad) = 1;
 RC = carlsonRC_core(x, y);
+RC(bad) = NaN;
 RC = reshape(RC, origSize);
 
 
@@ -73,8 +86,12 @@
 end
 
 if any(lt)
-    d      = sqrt((x(lt) - y(lt)) ./ x(lt));   % DLMF 19.2.18: (x-y)/x, not (x-y)/y
-    RC(lt) = atanh(d) ./ sqrt(x(lt) - y(lt));
+    % log((sqrt(x)+sqrt(x-y))/sqrt(y))/sqrt(x-y) == atanh(sqrt(1-y/x))/sqrt(x-y)
+    % without the 1 - sqrt(1-eps) cancellation (RC(3,1e-10) lost 8 digits).
+    % ... as log1p: log((sx+sxy)/sy) = log1p(((x-y)/(sx+sy) + sxy)/sy); the
+    % plain log lost 9 digits again for tiny x - y (RC(1+1e-13, 1)).
+    xl = x(lt);  yl = y(lt);  sx = sqrt(xl);  sy = sqrt(yl);  sxy = sqrt(xl - yl);
+    RC(lt) = log1p(((xl - yl)./(sx + sy) + sxy) ./ sy) ./ sxy;
 end
 
 

matlab/src/carlsonRD.m

@@ -23,6 +23,16 @@
 %   [2] B.C. Carlson, "Numerical Computation of Real or Complex Elliptic
 %       Integrals," Numer. Algorithms 10 (1995), 13–26.
 
+% Empty input -> empty output of the same shape (elementwise semantics; the
+% size checks below would otherwise reject [] against a scalar).
+if nargin >= 3 && (isempty(x) || isempty(y) || isempty(z))
+    sz = size(x);
+    if isempty(y), sz = size(y); end
+    if isempty(z), sz = size(z); end
+    RD = zeros(sz);
+    return;
+end
+
 if nargin < 3, error('carlsonRD: requires three arguments (x, y, z).'); end
 if ~isreal(x) || ~isreal(y) || ~isreal(z)
     error('carlsonRD: all input arguments must be real.');
@@ -32,7 +42,12 @@
 origSize = size(x);
 x = x(:).';  y = y(:).';  z = z(:).';
 
+% NaN, Inf and negative arguments give NaN (see carlsonRF).
+bad = ~(x >= 0 & y >= 0 & z > 0) | isinf(x) | isinf(y) | isinf(z);
+x(bad) = 1;  y(bad) = 1;  z(bad) = 1;
 RD = carlsonRD_core(x, y, z);
+RD(bad) = NaN;
+RD((x == 0) & (y == 0)) = Inf;          % diverges (DLMF 19.16.5)
 RD = reshape(RD, origSize);
 
 
@@ -44,19 +59,21 @@
 S     = zeros(size(x));
 fac   = ones(size(x));    % 4^{-n}
 
-for iter = 1:30
+% Per-element convergence: each element stops when IT has converged, so a
+% value never depends on what else is in the batch (a vector-wide break
+% made chunked and serial evaluations differ by an ulp).
+active = true(size(x));
+for iter = 1:200   % per-element break decides; cap guards pathological input
     lam  = sqrt(x.*y) + sqrt(y.*z) + sqrt(z.*x);
     sz   = sqrt(z);
-    S    = S + fac ./ (sz .* (z + lam));
-    fac  = fac ./ 4;
-    x    = (x + lam) ./ 4;
-    y    = (y + lam) ./ 4;
-    z    = (z + lam) ./ 4;
+    S(active)   = S(active) + fac(active) ./ (sz(active) .* (z(active) + lam(active)));
+    fac(active) = fac(active) ./ 4;
+    x(active) = (x(active) + lam(active)) ./ 4;
+    y(active) = (y(active) + lam(active)) ./ 4;
+    z(active) = (z(active) + lam(active)) ./ 4;
     A    = (x + y + 3.*z) ./ 5;
-    rng  = max([abs(x-A); abs(y-A); abs(z-A)]);
-    if rng < cr * min(A)
-        break;
-    end
+    active = active & (max([abs(x-A); abs(y-A); abs(z-A)], [], 1) >= cr * A);
+    if ~any(active), break; end
 end
 
 A  = (x + y + 3.*z) ./ 5;

matlab/src/carlsonRF.m

@@ -27,6 +27,16 @@
 %   [2] B.C. Carlson, "Numerical Computation of Real or Complex Elliptic
 %       Integrals," Numer. Algorithms 10 (1995), 13–26.
 
+% Empty input -> empty output of the same shape (elementwise semantics; the
+% size checks below would otherwise reject [] against a scalar).
+if nargin >= 3 && (isempty(x) || isempty(y) || isempty(z))
+    sz = size(x);
+    if isempty(y), sz = size(y); end
+    if isempty(z), sz = size(z); end
+    RF = zeros(sz);
+    return;
+end
+
 if nargin < 3, error('carlsonRF: requires three arguments (x, y, z).'); end
 if ~isreal(x) || ~isreal(y) || ~isreal(z)
     error('carlsonRF: all input arguments must be real.');
@@ -36,7 +46,15 @@
 origSize = size(x);
 x = x(:).';  y = y(:).';  z = z(:).';
 
+% NaN, Inf and negative arguments give NaN (R_F is defined for x, y, z >= 0);
+% they used to reach the duplication and come back as complex NaN.
+bad = ~(x >= 0 & y >= 0 & z >= 0) | isinf(x) | isinf(y) | isinf(z);
+x(bad) = 1;  y(bad) = 1;  z(bad) = 1;
 RF = carlsonRF_core(x, y, z);
+RF(bad) = NaN;
+% Two zero arguments: the integral diverges (DLMF 19.16.1); the duplication
+% loop just stalls and returned a finite 2e6 for R_F(0, 0, 1).
+RF((x == 0) + (y == 0) + (z == 0) >= 2) = Inf;
 RF = reshape(RF, origSize);
 
 
@@ -49,15 +67,20 @@
 
 x0 = x;  y0 = y;  z0 = z;
 
-for iter = 1:20
+% The adaptive break decides; the cap only guards pathological input
+% (20 was too few for R_F(0, 1e-16, 1) and every K(m) at tiny m).
+% Per-element convergence: each element stops when IT has converged, so a
+% value never depends on what else is in the batch (a vector-wide break
+% made chunked and serial evaluations differ by an ulp).
+active = true(size(x));
+for iter = 1:200
     lam = sqrt(x.*y) + sqrt(y.*z) + sqrt(z.*x);
-    x   = (x + lam) ./ 4;
-    y   = (y + lam) ./ 4;
-    z   = (z + lam) ./ 4;
+    x(active) = (x(active) + lam(active)) ./ 4;
+    y(active) = (y(active) + lam(active)) ./ 4;
+    z(active) = (z(active) + lam(active)) ./ 4;
     A   = (x + y + z) ./ 3;
-    if max(max(abs(x - A)), max(max(abs(y - A)), abs(z - A))) < cr * min(A)
-        break;
-    end
+    active = active & (max([abs(x - A); abs(y - A); abs(z - A)], [], 1) >= cr * A);
+    if ~any(active), break; end
 end
 
 A  = (x + y + z) ./ 3;