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

Open
  • Is something wrong?
#545 3 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
55/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Active
Tech stack
julia
Domain
backend, testing

Research direction

Start by reading the real(z) < 0 branch of _expint in src/expint.jl, especially the imstart, nsteps, and Taylor-walk logic. Then inspect test/expint.jl and reproduce the reported Float32 and Float16 cases against Float64 or BigFloat references. Done means low-precision results near the negative real axis are accurate while the existing Float64 behavior remains correct.

Written by the indexing model from the issue text.

Description

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.

Dominant language
Julia
Stars
384
Forks
115
Avg merge
2d 9h
Merged PRs (30d)
2

Getting set up

This project ships no dev container, Dockerfile or contributing guide, so setting up is up to you: start from its README, and see our first-contribution guide for the general steps.

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from JuliaMath/SpecialFunctions.jl

All issues in JuliaMath/SpecialFunctions.jl

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.