[NAThompson]

Adds the Struve function H_ν(x) and the modified Struve function L_ν(x) for real order and argument, with optional x-derivatives. "Struve H function" is on the feature wishlist in #303.

API. struve_h(nu, x), struve_h(nu, x, &derivative), struve_h(nu, x, pol) and struve_h(nu, x, &derivative, pol), and the same for struve_l. Negative x is allowed for integer ν (by reflection); at x = 0 the result is 0, 2/π at ν = −1, or an overflow error where it is infinite.

Method (switch points chosen from per-method error maps against MPFR):

• Power series (DLMF 11.2.1/2) for small x, accepted only when it does not cancel.
• H for ν ≥ −1/2 below the transition region: the Neumann series in J_{ν+1/2+k} (11.4.18), whose terms are all positive, with the J's from a normalized Miller recurrence. This covers large ν, where the power series cancels and |Y_ν| ≫ |H_ν|.
• Otherwise H = K + Y_ν, with K from the asymptotic expansion (11.6.1) or the Laplace integral (11.5.2). The integral uses exp_sinh with a new opt-in compensated summation, exp_sinh(max_refinements, compensated = true); with the default false the integrator is bit-identical to before (checked on 44 integrals, including complex, norm and finite-interval ones) and as fast.
• L at large x: L = I_ν + M_ν with M from 11.6.2, used only where the omitted multiple of K_|ν| is below rounding.
• Γ(ν + 3/2) in double-length arithmetic (tgamma alone is up to 11 ulps off), and the rounding of ν + 1 and ν + 3/2 carried to first order.

Accuracy (double, against MPFR at the exact double arguments; the envelope is max(1/2, κ/2) ulps with κ = |x f′/f|). Over ν ∈ [−20, 100] × x ∈ [10⁻², 10⁴], 282 of 38,400 pixels of H and 167 of L are outside the envelope. What remains:

• For small ν and x ≈ 2–60, where |Y_ν| ~ |H_ν|, H inherits the error of cyl_neumann (about 7 ulps for Y₀ and Y₁ and up to 13 for negative orders).
• Within a few hundredths of ν = 1/2, H has near-double zeros and K + Y cancels (ν = 1/2 itself is a closed form).
• cpp_bin_float_50 is within about 55× the envelope; float is computed in double.

struve_h, |ulps| / envelope
struve_l, |ulps| / envelope
struve_h(1, x)
struve_h(40.5, x)
struve_h(-7.3, x)
struve_l(-7.3, x)
struve_l(10, x)

Performance (google/benchmark, Apple M1 Pro, -O2, reporting/performance/struve_performance.cpp):

path	time
H, power series	100 ns
H, asymptotic	98 ns
H, Neumann series (ν = 40)	0.98 µs
H, integral	2.4 µs
H, negative order	2.7 µs
H with derivative (integral)	2.9 µs
L, power series	114 ns
L, asymptotic	55 ns
cyl_neumann, for comparison	16 ns

exp_sinh with compensated = true costs 8–28% more than without on the existing quadrature benchmarks.

Under CUDA, exp_sinh_detail.hpp and sinh_sinh_detail.hpp defined __constant__ tables with the same names in the same namespace, so no translation unit could include both; special_functions.hpp now does (through struve.hpp), so exp_sinh's CUDA section moves into its own nested namespace.

References are generated by reporting/accuracy/struve_reference.jl (BigFloat series, checked at two precisions and against an independent MPFR implementation).

🤖 Generated with Claude Code

[NAThompson]

@matwey: Please take a look at your convenience-see either here or test the wheel or wasm here.