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

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`:

```julia
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.