`expint(nu,z)` has bad precision for `nu != 1`
Nobody has claimed this yet.
Assessment
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Newbie friendliness
- 48/100
Research direction
Start by running the Julia comparison code in the issue with SpecialFunctions and QuadGK, reproducing the relative-error plot for nu = 2 over z from 0.01 to 60. Then locate the expint(nu,z) implementation and compare its results with the quad-precision benchmark. Done means materially reducing the systematic precision errors for nu values other than 1 while preserving the existing nu == 1 behavior.
Written by the indexing model from the issue text.
Description
The function expint(nu,z) presents systematic large errors for nu != 1 (for nu == 1 the library calls an optimized routine that in my experiments has good precision). I computed the relative error of expint(nu,z) for $z \in [0.01, 60]$ using as benchmark numerical quadrature in quad precision. Code for nu == 2:
import Plots as plt
import QuadGK as GK
using SpecialFunctions
function compare_plot()
setprecision(BigFloat, 128)
int_tol = BigFloat(1e-24)
nu = Float64(2)
f0 = s -> GK.quadgk(u -> exp(-s*u) / u^nu, BigFloat(1), BigFloat(Inf);
rtol=int_tol, order=15)[1]
f = s -> expint(nu, Float64(s))
srange = 0.01:0.02:60
plt.gr(size=(160*3,120*3), legend=:none)
p1 = plt.plot(srange, BigFloat.(f.(srange)) ./ f0.(BigFloat.(srange)) .- 1)
plt.display(p1)
end
compare_plot();
The code outputs this plot:
Similar plots are obtained for other values of nu (except for 1, as mentioned before), and in general the errors have the same form, with peaks of the order of 1e-13 in magnitude.
- 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
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
More from JuliaMath/SpecialFunctions.jl
-
Difficulty 2/5 1-3 hours Newbie friendliness 68/100
JuliaMath/SpecialFunctions.jl#570 ·
-
Difficulty 5/5 Over a week Newbie friendliness 40/100
JuliaMath/SpecialFunctions.jl#569 · 2 comments ·
-
Difficulty 4/5 3-5 days Newbie friendliness 25/100
JuliaMath/SpecialFunctions.jl#567 ·
-
Difficulty 4/5 1-2 days Newbie friendliness 56/100
JuliaMath/SpecialFunctions.jl#563 ·
-
Difficulty 3/5 1-2 days Newbie friendliness 72/100
JuliaMath/SpecialFunctions.jl#560 ·
All issues in JuliaMath/SpecialFunctions.jl
Similar issues
-
Difficulty 2/5 1-3 hours Newbie friendliness 74/100
Maintainers usually reply within 1 day
-
Difficulty 2/5 1-3 hours Newbie friendliness 74/100
-
Difficulty 2/5 1-3 hours Newbie friendliness 67/100
-
Difficulty 2/5 1-3 hours Newbie friendliness 78/100
-
backport 1.14 GC
Difficulty 2/5 1-3 hours Newbie friendliness 66/100
Maintainers usually reply within 1 day