[NAThompson]

Eigenvalues, angular functions (DLMF normalization), radial functions of both kinds and derivatives,
tested against Arb references; see #303 and #1630.

Requested by @dschmitz89 in #303 for SciPy, whose specfun versions are inaccurate (pro_rad2, obl_rad2 have no correct digits at many points).

Conventions. The angular functions are DLMF's Ps_n^m(η, ±c²) (Meixner–Schäfke: the norm and Condon–Shortley phase of the Ferrers functions, sign fixed at η = 0), as discussed in #1630. Flammer's normalization (SciPy's pro_ang1/obl_ang1) is a per-(n, m, c) rescaling, documented in the docs; it was avoided because for oblate functions with large c the value at 0 is exponentially small (norms near 1e24 at c = 60), so it cannot even be computed reliably in double. The eigenvalue is DLMF's λ (SciPy's pro_cv = λ + c², obl_cv = λ − c²). The radial functions agree in all conventions.

API. {prolate,oblate}_spheroidal_eigenvalue(n, m, c), {prolate,oblate}_spheroidal_angular(n, m, c, eta, T* derivative = nullptr) and {prolate,oblate}_spheroidal_radial_{1,2}(n, m, c, xi, T* derivative = nullptr). c may be negative: the eigenvalues and angular functions are even in c, and R1, R2 change by (−1)ⁿ, (−1)ⁿ⁺¹, as in Mathematica. The derivative shares the expensive part of the computation (as SciPy's functions return both); it uses forms that avoid the cancellation autodiff of the value inherits near ξ = 1 (e.g. 2079 ε vs ~2 ε for R1′ at ξ = 1.0001). prolate_spheroidal<Real>(n, m, c) and oblate_spheroidal<Real>(n, m, c) compute the eigenvalue and coefficients once and have members eigenvalue(), angular(eta, &d), radial_1(xi, &d), radial_2(xi, &d), bit-identical to the free functions, which are implemented through them.

Methods. Sturm bisection + Newton on Bouwkamp's continued fraction for λ; Legendre coefficients from both continued fractions, matched at the largest coefficient. The classic Bessel expansion of R1 cancels catastrophically for prolate functions (16 digits at c = 50), so prolate R1 uses Van Buren & Boisvert's expansion in j(c√(ξ²−1)); oblate is the mirror image. R2 comes from the Neumann series where it is well-conditioned and otherwise from an asymptotic expansion with phase c√(ξ²−τ), τ = ±1 − λ/c², then a Taylor-series integration inward. Where the Legendre series of the angular function cancels (oblate near η = 0, prolate near η = ±1), the angular equation is integrated from η = 0 or from the power series at η = 1 and matched.

Accuracy (worst error on a 1500-point Arb grid, c ≤ 50, scaled by the local envelope max(|f|, |x f'|, |f'|/k)):

	prolate c ≤ 8	prolate c ≤ 50	oblate c ≤ 8	oblate c ≤ 50
eigenvalue	1 ε	1 ε	1 ε	1 ε
angular	6 ε	16 ε	10 ε	15 ε
angular derivative	26 ε	36 ε	51 ε	49 ε
radial, first kind	3 ε	4 ε	6 ε	10 ε
radial, second kind	2 ε	3 ε	6 ε	25 ε

The public eigenvalues are correctly rounded except in a few cases (worst 0.84 ulp). At large c, λ ≈ ∓c² carries an absolute rounding error near εc², which biases solutions where they grow and cancels in the low-order Legendre coefficients. So, as Arb does with guard bits, a double-length low part of λ is found from one Newton step on a double-double Bouwkamp residual, and the coefficient ratio chain is rerun in double-double. Near ξ = 1, prolate R2 at large c is a small multiple of R1, so any multiple of R1 the inward integration picks up is amplified, and the Wronskian rescale cannot remove it. For prolate R2 with c > 8 and ξ − 1 < 1/8 the integration therefore runs in double length from the asymptotic expansion at ξ = min(1 + c/8, 4): 26 ε at c = 50, ξ − 1 = 10⁻³ (Arb), and at most 6.4× the ξ − 1 envelope on a grid c ∈ [12, 64], ξ − 1 ∈ [10⁻⁴, 10⁻¹]; what remains comes from the asymptotic start value.

	prolate	oblate
angular
R1
R2
R2 error / envelope over (c, ξ)

Performance (google/benchmark, Apple M1 Pro, -O2, reporting/performance/spheroidal_performance.cpp): one call costs 16–32 µs for the eigenvalue at c = 1–64 and 20–100 µs for the angular and radial functions; prolate R2 near ξ = 1 costs 0.22 ms at c = 16 and 0.49 ms at c = 64 (double length; 1.8× and 5.4×), and oblate R2 near ξ = 0 0.46 ms and 1.2 ms (eigenvalue splitting). Over a 64-point sweep (angular and R1) the class is 21× faster than the free functions at c = 1–4, 10× at c = 16 and 4.8× at c = 64. The reference generator is reporting/accuracy/spheroidal_reference.jl.

🤖 Generated with Claude Code

[dschmitz89]

Thanks @NAThompson .

@steppi Could you take a look here? Boost will add spheroidal wave functions. In its current form, is the code compatible with your new gufunc paradigm? From what I understand, you prefer a clear distinction between coefficient computation kernels and the Eigenvalue computation.

[NAThompson]

@dschmitz89 : Check the results here -the new special functions are in the dropdown.

[dschmitz89]

@dschmitz89 : Check the results here -the new special functions are in the dropdown.

That site is incredible. Is its code/frontend purely LLM based?

[NAThompson]

Is its code/frontend purely LLM based?

Yup-all I did was invoke the frontend design skill and then added some nitpicks about styling and color tables. The C++ is evaluated via the WASM module.

[NAThompson]

@dschmitz89 : Just a heads up: We'll follow the Mathematica argument ordering convention:

SpheroidalEigenvalue[2,1, 1]
5.4246991437751347200272962780281450832808406774641549564802859050...

>>> import boost.math as bm
>>> n = 2, m =1, c = 1
>>> bm.prolate_spheroidal_eigenvalue(n,m,c)
5.424699143775134

scipy currently switches the ordering of the first two arguments and differs from Mathematica/Boost by c^2:

>>> import scipy
>>> scipy.special.pro_cv(m,n,c)
np.float64(6.424699143775136)

. . . which I imagine your LLM will have no problem resolving as the convention ambiguity is also stated in DLMF.

[NAThompson]

Build failure is:

--2026-11-05 04:58:39--  https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB
Resolving apt.repos.intel.com (apt.repos.intel.com)... 23.198.69.169, 2600:1407:3c00:148f::a87, 2600:1407:3c00:1491::a87
Connecting to apt.repos.intel.com (apt.repos.intel.com)|23.198.69.169|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 4738 (4.6K) [application/vnd.exstream-package]
Saving to: ‘GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB’

     0K ....                                                  100%  131M=0s

2026-11-05 04:58:40 (131 MB/s) - ‘GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB’ saved [4738/4738]

sudo: 'apt-key': command not found
Error: Process completed with exit code 1.

merging over it.