expint truncates its series after a number of terms fixed for Float64, silently capping BigFloat accuracy

[andreasnoack]

expint (and everything that routes through it) truncates its series and continued fractions after a number of terms that is fixed for Float64, while the methods accept any number type. For BigFloat the result therefore silently stops improving well before the requested precision is reached. There are two independent limits, plus a third place where a fixed-order expansion is (correctly) restricted but leaves BigFloat without a substitute.

All numbers below are SpecialFunctions v2.8.3 on Julia 1.12.7. The reference values come from MPFR's expinti via E₁(x) = -Ei(-x), E₁(-x + 0im) = -Ei(x) - πi and E₂(z) = exp(-z) - z E₁(z), so they are independent of the code under test.

1. niter = 1000 caps the continued fraction at ~92 digits

julia> using SpecialFunctions

julia> setprecision(BigFloat, 512) do            # eps(BigFloat) = 1.5e-154
           x = BigFloat(3)
           ref = exp(-x) + x*expinti(-x)         # = E₂(3)
           (default   = Float64(abs(expint(2, x)          - ref)/abs(ref)),
            niter_1e5 = Float64(abs(expint(2, x, 100_000) - ref)/abs(ref)))
       end
(default = 1.1031886851560754e-92, niter_1e5 = 1.3250352104469716e-152)

expint/expintx/_expint take niter::Int = 1000. The continued fraction is slowest just above the cutoff abs(z) == 3, where it gains only about 10*sqrt(n) bits from n terms:

niter	1000	2000	4096	8000	16000
bits of expint(2, big(3)) correct	305	436	629	883	≥1024

so the number of terms needed grows quadratically with the precision: En_cf(big(2), big(3), 400_000) reports convergence after 10562 iterations at 1024 bits and 168194 at 4096 bits. 1000 terms are ample for Float64 but cap every higher-precision result at roughly 92 digits.

2. En_taylor is a hard 101-term series, capping the negative real axis at ~97 digits

julia> setprecision(BigFloat, 512) do
           y = BigFloat(4)
           ref = exp(y) + y*(-expinti(y) - im*BigFloat(π))   # = E₂(-4 + 0im)
           z = Complex{BigFloat}(-y, 0)
           (default   = Float64(abs(expint(2, z)          - ref)/abs(ref)),
            niter_1e5 = Float64(abs(expint(2, z, 100_000) - ref)/abs(ref)))
       end
(default = 5.2022497324938007e-98, niter_1e5 = 8.768948249862689e-101)

Note that raising niter does not help here: the limit is the

for k = 0:100

loop in En_taylor, which is used by the real(z) < 0 procedure to step towards the axis. It has a proper ϵ = 10*eps(real(asum)) convergence check but a term count fixed for Float64. Raising just that bound to 1000 takes the same computation from 323 to 1017 correct bits at 1024-bit precision, at no extra cost in time.

3. The near-pole δ-expansion has no high-precision counterpart

In En_expand_origin_general the correction used when ν is close to a positive integer is a fixed 5-term expansion in δ = round(ν) - ν (orders δ⁰ … δ⁴), and it is guarded by

if real(ν+z) isa Union{Float64, Float32} && abs(gammaterm - blowup) < 1e-3 * abs(blowup)

That guard is the right call — with only five terms the truncation error O(δ⁵) is not below eps(BigFloat) — but it means BigFloat falls back to the cancelling gammaterm - (blowup + sumterm) form and loses about log2(1/δ) bits:

expint(2 + 1e-20, big(0.5)) at 128-bit precision:  61.8 of 128 bits correct
expint(2 + 1e-40, big(0.5)) at 512-bit precision: 382.6 of 512 bits correct

Extending the series is unattractive (the coefficients are polynomials in polygamma values). For BigFloat the natural fix is guard bits instead: evaluate the plain expression with the working precision raised by ceil(log2(1/abs(δ))). There is precedent for that in the package — loggamma(z::Complex{BigFloat}) already guardrails its precision.

Suggested fixes for 1 and 2

Since every loop stops as soon as its eps-based criterion is met, a generous limit costs nothing for arguments that converge quickly. Making the two limits scale with the precision of the arguments is enough:

# quadratic, because the continued fraction gains ~10*sqrt(n) bits from n terms
_expint_niter(ν, z) = max(1000, div(precision(float(real(promote_type(typeof(ν), typeof(z))))), 8)^2)

# linear is ample here, the series converges factorially
for k = 0:max(100, precision(real(asum)))     # En_taylor

precision is 53, 24 and 11 for Float64, Float32 and Float16, so all three defaults evaluate to the current literals and hardware-float results are unchanged (I checked 144 values of expint/expintx/gamma(a,x) across real/complex Float16/32/64 orders and arguments: bit-identical). With both limits scaled, the three probes above reach ~100*eps(BigFloat) at 512 bits, and the negative-axis case goes from 323 to 2045 correct bits at 2048-bit precision. It is also faster at 1024 bits (13 s versus 20 s), because the old code compensated for the small quick_niter = niter >> 4 by starting much further off the axis and taking more Taylor steps. That path remains expensive at high precision, though: about 2 s at 512 bits and 190 s at 2048 bits for expint(2, Complex{BigFloat}(-4)).

Related: non-convergence is silent

En_cf returns the iteration count it used, and _expint discards it (g, cf, _ = En_cf(ν, z, niter)). Once the limit tracks the requested precision, reaching it means genuine non-convergence rather than "Float64 was already enough", so it would be a good place to warn or throw instead of returning a quietly inaccurate result.

For context, I went looking for this pattern across the package and expint is the only place where it produces silently inaccurate results. _ellipk/_ellipe (Float64-only), sinint/cosint (explicit error for other AbstractFloats), _jinc (Taylor shortcut restricted to Union{Float32,Float64} with per-type thresholds), @E₁_cf64/@E₁_taylor64 (gated on Float64/ComplexF64), and the beta_inc*, _gamma_inc, _gamma_inc_inv, erfi/dawson/faddeeva internals all either restrict their signatures or route BigFloat to MPFR.

---

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

[stevengj]

I'm confused, the Float64 implementation shouldn't be used for BigFloat, since we have a specialized BigFloat method that calls mpfr_eint

SpecialFunctions.jl/src/expint.jl

Lines 144 to 148 in af329fa

	function expint(x::BigFloat)
	iszero(x) && return Inf
	signbit(x) && throw(DomainError(x, "negative argument, convert to complex first"))
	return -expinti(-x)
	end

SpecialFunctions.jl/src/expint.jl

Lines 608 to 612 in af329fa

	function expinti(x::BigFloat)
	z = BigFloat()
	ccall((:mpfr_eint, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode), z, x, Base.MPFR.ROUNDING_MODE[])
	return z
	end

[stevengj]

Oh, I see, we are talking about the 2-argument version, expint(n, x)