Add a type stability sweep over the exported functions (supersedes #43)

[andreasnoack]

Reworks the sweep proposed by @giordano in #43, which has been waiting for a rebase since 2017. The idea is his; the implementation is new, because the file it patched (test/runtests.jl) has since been split into per-topic files.

Based on #548 — 12 of the assertions here depend on it. Please merge that one first; this PR targets its branch and GitHub will retarget it to master afterwards.

What it checks

A new test/type_stability.jl asserting

@test @inferred(f(x)) isa T

over a grid of Float16, Float32, Float64, ComplexF32 and ComplexF64 arguments, for the error functions, the gamma family (including zeta(s, z), polygamma, beta/logbeta, gamma(a, x)), the incomplete gamma and beta functions, the Airy functions plain and scaled, the elliptic integrals, the sine and cosine integrals, the exponential integrals and the Bessel functions.

Asserting the type rather than only that inference succeeds is the part that differs from #43, and it matters: of the 43 methods that get this wrong, only 14 fail a bare @inferred. The other 29 are inferred perfectly and simply return a Float64 — for example besseli(2, Float16(1)) === 0.13574766976703834.

The rule being encoded is that the result carries the precision that promoting the arguments gives: an Integer order does not by itself force Float64, and a Float16 or Float32 argument is never widened.

Result

243 passing, 43 broken, 0 failing, identical on 1.10, 1.11, 1.12 and 1.13-rc3, so no @test_broken unexpectedly passes on any supported version.

Every broken entry is in the Bessel family and they share one cause, documented at the top of the file: the fall-through in

function besselj(nu::Real, x::AbstractFloat)
    if isinteger(nu)
        if typemin(Cint) <= nu <= typemax(Cint)
            return besselj(Cint(nu), x)
        end
    ...
    real(besselj(float(nu), complex(x)))
end

promotes the order with float(nu), which is Float64 for an Integer order, and drags the argument up with it. Hence besselj(2, 1.0f0) is Union{Float32,Float64} — the two branches disagree — while besselj(2.0f0, 1.0f0) is fine. The breakdown is 12 Float32, 12 Float16 and 19 ComplexF32 entries.

besselj(2, Float16(1)) and bessely(2, Float16(1)) are deliberately not called: the same fall-through recurses until the stack overflows, which would leave the test process in a corrupted state. That is #547, and fixing it should clear most of this list at the same time. See also #234.

The intent is that the broken entries be turned into passing ones one at a time; @test_broken fails loudly once a case starts working, so nothing can be fixed silently without the list shrinking.

Not covered yet: Complex{Float16} arguments, and BigFloat.

---

Prepared by Claude Code on behalf of @andreasnoack; the investigation and the text above are Claude's.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.67%. Comparing base (adbeb4b) to head (c26740d).
⚠️ Report is 8 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #549      +/-   ##
==========================================
+ Coverage   94.49%   94.67%   +0.17%     
==========================================
  Files          14       14              
  Lines        3016     3023       +7     
==========================================
+ Hits         2850     2862      +12     
+ Misses        166      161       -5     

Flag	Coverage Δ
unittests	94.67% <ø> (+0.17%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ 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.