Take same-precision arguments in expint's internal helpers, restrict the incomplete gamma

[andreasnoack]

Stacked on #549 — please merge that one first; this PR targets its branch and GitHub will retarget it to master afterwards. It flips the twelve expint/gamma(a, x) entries in the sweep from @test_broken to @test.

The problem

The internal helpers took (n::Real, z::Number) and (ν::Number, z::Number), so an order and an argument of unrelated precisions met inside them and every operation adopted whichever operand it touched: loggamma(n+1), logabsgamma(ν) and polygamma followed the order — Float64 for an Integer one — while the rest followed the argument. _expint promoted z against ν but deliberately left ν alone, so the mismatch was the normal case, not an edge case.

Consequences: expint(::Float32, ::Float32) and gamma(::Float32, ::Float32) inferred Union{Float32,Float64} (what #520 ran into), a BigFloat result on the negative real axis was capped at Float64 accuracy, and Float16 returned 0.0 near a positive integer order.

The change

Every internal helper takes its arguments at one precision — Union{T,Complex{T}} with T one of Float16, Float32, Float64 — and the promotion happens once, at the public entry. No promoting methods on the helpers, so a mismatch is a MethodError at the call site rather than a silent mixed-precision computation.

• An integer order now adopts the argument's precision instead of dragging the computation to Float64.
• En_safe_expfact(n::T, z::T) / (n::T, z::Complex{T}). isodd works on floats, so Int(n), the n <= typemax(Int) guard and the (-1)^n fallback are gone, and with them the typemax(Int) dispatch in En_expand_origin. The switch to the closed form moves from n < 100 to n < 12: measured, the product form is the more accurate of the two everywhere, and above n ≈ 12 the term is negligible next to the sums the callers add it to, in every precision.
• En_expand_origin_general keeps the near-pole correction in a separate En_origin_pole_series, defined only for the precisions its δ⁴ truncation can serve, and triggered by abs(δ) < eps(T)^(1/6)/2 instead of a cancellation ratio: the series error is O(δ⁵) and the cancelling form loses accuracy like eps/δ, so the two meet at δ ~ eps(T)^(1/6). The old ratio saturated at ~0.036 for δ ≥ 0.1 and could not express the Float16 boundary. Its π^2/π^4 constants are evaluated in T.
• En_safeexpmult and En_imagbranchcut are split into real and complex methods rather than branching on isa, and use sincospi.
• The unreachable En_safe_gamma_term(::Integer, ::Real) method and the float/promote/oftype shims inside the helpers are gone.

Incomplete gamma

gamma(a, x) and loggamma(a, x) now throw DomainError for real x < 0, where the value is complex except at integer a. This is a breaking change: gamma(2, -1.3) returned -1.1008 and now throws, and the negative-x assertions in test/gamma_inc.jl — including the hand-computed big"-1.1007890…" and big"-2.0801…e101" references — become @test_throws.

The rationale is that both available kernels reject it anyway: MPFR's mpfr_gamma_inc returns NaN for x < 0, and NSWC's GRATIO documents "IT IS ASSUMED THAT A AND X ARE NONNEGATIVE". The old support came from _gamma_big routing negative x with integer a through x^a * expint(1-a, x), i.e. through the very function this PR is fixing.

So _gamma_big is deleted along with its four _gamma methods and the separate gamma(a::Integer, x); the BigFloat method is now a plain mpfr_gamma_inc ccall. Only x is restricted, not a: expint delegates as gamma(1-ν, z), so the first argument is negative whenever ν > 1, and mpfr_gamma_inc accepts a negative first argument — it is only x < 0 that returns NaN.

Fixes

Float16 near a positive integer order returned zero, because the old real(ν+z) isa Union{Float64,Float32} guard excluded it and the fallback then adds blowup ≈ 1/δ ≈ 500 to sumterm ≈ 0.057 — Float16's spacing at 500 is 0.5, so the sum is annihilated:

julia> expint(Float16(3.002), Float16(1.5))     # true value 0.0567
Float16(0.0)                                    # master
Float16(0.05627)                                # this PR

Float16 was never a numerical problem there: at its spacing near a pole the truncation error is δ⁵ ≈ 3e-14, eleven orders below eps(Float16).

Type stability, all now concrete: expint(::Int,::Float16), expint(::Int,::Float32), expint(::Float32,::Float32), expint(::Int,::ComplexF32), expint(::Float32,::ComplexF32), gamma(::Float32,::Float32).

BigFloat on the negative real axis — from the signature, not a conversion: expint(5, big(-4)+0im) goes from 53 to 244 correct bits of 256, expint(9, …) from 50 to 245.

Real BigFloat near a positive integer order is delegated to MPFR via E_ν(z) = z^(ν-1) Γ(1-ν, z), which has no cancellation there: 258 correct bits of 256, against 119 for the series at δ = 1e-40. Validated against an independent reference (the series at 4096 bits), not against itself.

expintx has no BigFloat method any more. MPFR has no scaled incomplete gamma, and building exp(z)·E_ν(z) from the unscaled value would defeat the point of the scaled function — it forms the underflowing quantity and multiplies by an overflowing one. Complex BigFloat has no route either (polygamma is undefined for it, MPFR is real only, and gamma(a,x) for complex arguments routes back through expint), so it throws rather than returning a badly cancelled value. A follow-up can restore both with guard bits: recomputing at precision + ceil(log2(1/abs(δ))) + 10 and rounding back gives full precision.

Also

JET.report_package loses the _polygamma(::Int64, ::BigFloat) finding, since the correction no longer reaches polygamma with a BigInt-derived order.

Full suite passes on 1.10, 1.12 and 1.13-rc3. New tests cover the near-pole values for all three hardware precisions, the BigFloat delegation, and the complex BigFloat error.

Known limitation, not addressed here

The origin series and the continued fraction lose accuracy to cancellation independently of any of this: at x = 2.94 the two terms are 119–255× the result, so ~8 bits go structurally. It is invisible at Float64 but not at Float32, where it makes expint(1, x) and expint(x) disagree by up to 439 ulps (the one-argument entry evaluates E₁ in Float64 and rounds). Filing separately.

---

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

[codecov]

Codecov Report

❌ Patch coverage is 98.70130% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 94.68%. Comparing base (b2a7190) to head (87a4eb6).

Files with missing lines	Patch %	Lines
src/expint.jl	98.46%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #550      +/-   ##
==========================================
+ Coverage   94.67%   94.68%   +0.01%     
==========================================
  Files          14       14              
  Lines        3023     3010      -13     
==========================================
- Hits         2862     2850      -12     
+ Misses        161      160       -1     

Flag	Coverage Δ
unittests	94.68% <98.70%> (+0.01%)	⬆️

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.

[andreasnoack]

@stevengj Would you be able to take a look? The problem I'm trying to solve here is partly type stability, but also precision. The solution is to restrict methods to have matching floating-point types for both arguments and only allow types for which we can ensure that the results have the expected precision. Hence, I had to drop the claim that the implementation supported arbitrary float inputs since I don't think the implementation actually supported that claim, see #546. Instead, I'm using MPFR when possible and throwing method errors rather than results with large errors.