expint(nu,z) has bad precision for nu != 1

[djturizo]

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.