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

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

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

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.