expint is inaccurate for Float32/Float16 arguments near the negative real axis

[andreasnoack]

expint(ν, z) loses most of its accuracy for Float32 and Float16 arguments when z is on or close to the negative real axis. The error is far above the resolution of the type: up to ~85% relative error, i.e. the result is simply wrong rather than merely rounded.

Tested with SpecialFunctions v2.8.3 on Julia 1.12.7.

julia> using SpecialFunctions

julia> expint(2, ComplexF32(-3.5))
-2.9899444580078125 - 10.99557399108051im     # ← wrong

julia> expint(2, ComplexF64(-3.5))
-15.623287024340843 - 10.995574287564276im

julia> setprecision(BigFloat, 128) do; expint(2, Complex{BigFloat}(-7//2, 0)); end
-15.62328702434085976994732524631780778557 - 10.99557428756427633461925184147826009474im

so the Float64 result is correct and the Float32 one is off by a factor of 5 in the real part. Float16 is likewise wrong (-3.19140625 - 10.994015112188281im).

A few more points, Float32 versus Float64 input:

ν	z	Float32 input	Float64 input	rel. error
1	-3.5f0 + 0.0f0im	-11.107379f0 - 3.1415927f0im	-13.92535399515229 - 3.141592653589793im	0.197
2	-3.5f0 + 0.0f0im	-2.9899445f0 - 10.995574f0im	-15.623287024340843 - 10.995574287564276im	0.661
2	-3.5f0 + 0.01f0im	-3.0192275f0 - 0.69151247f0im	-15.654229875324825 - 10.856321873984799im	0.851
3	-3.5f0 + 0.0f0im	-8.832054f0 - 19.242254f0im	-10.783026313250343 - 19.24225500323748im	0.0884
2	-4.0f0 + 0.0f0im	-25.642519f0 - 12.566371f0im	-23.925347847080612 - 12.566370614359172im	0.0635
2	-10.0f0 + 0.0f0im	-2854.4534f0 - 31.415928f0im	-2895.823967611959 - 31.41592653589794im	0.0143
1	-5.0f0 + 0.0f0im	-39.963936f0 - 3.1415927f0im	-40.18527535580315 - 3.141592653589793im	0.00549
1.5	-10.0f0 + 0.0f0im	-2653.0286f0 - 11.209983f0im	-2676.2869535640534 - 11.209982432795858im	0.00869

The order ν can be an Integer or a float of the same type as z; it makes no difference. Points with a positive real part are fine (relative error ~1e-7 for Float32), and Float64 and BigFloat are fine everywhere I looked, so this is specific to the low-precision types on the negative real axis.

Where it seems to come from

The real(z) < 0 branch of _expint (the Amos procedure in src/expint.jl) walks in from a point with a larger imaginary part:

imstart = (imz == 0) ? abs(z)*sqrt(eps(typeof(real(z)))) : imz
z₀ = rez + imstart*im
...
nsteps = ceil(2 * (imstart - imz))
Δ = (imz - imstart)*im / nsteps
for j = 1:nsteps
    E_start = En_taylor(ν, E_start, z₀, Δ)
    z₀ += Δ
end

Several parts of this are tied to the precision of the argument in a way that does not hold up for Float32/Float16:

• imstart is abs(z)*sqrt(eps(T)), i.e. ~3.4e-4 for Float32 and ~3.4e-2 for Float16, versus ~1.5e-8 for Float64. The starting point is much further from the axis.
• nsteps = ceil(2*(imstart - imz)) is then 1 for all these cases, so the whole distance is covered in a single Taylor step of En_taylor.
• the while i == quick_niter loop that doubles imstart uses niter >> 4 iterations irrespective of the precision.

So the number of steps and the truncation of the Taylor step are effectively tuned for Float64 while the starting distance grows as sqrt(eps(T)), and for Float32/Float16 the walk back to the axis is too coarse.

Not covered by the tests

test/expint.jl exercises the negative real axis only in Float64, so nothing catches this. Some Float32/Float16 points with real(z) < 0 compared against a BigFloat or Float64 reference would.

Note that the return type is currently widened to ComplexF64 even for ComplexF32 input, which hides the problem in a typeof check; that part is a separate (type stability) issue.

---

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

[stevengj]

I'm confused, the expint implementation is only designed for Float64, whereas it explicitly defines method for Float16 and Float32 that convert to Float64 and then convert the result back:

SpecialFunctions.jl/src/expint.jl

Lines 152 to 159 in af329fa

	for f in (:expint, :expintx)
	for T in (Float16, Float32)
	@eval $f(x::$T) = $T($f(Float64(x)))
	end
	for CT in (ComplexF16, ComplexF32)
	@eval $f(x::$CT) = $CT($f(ComplexF64(x)))
	end
	end

[stevengj]

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