expint(ν, z) loses 8-17 bits to cancellation; one- and two-argument forms disagree by up to 439 ulps in Float32

The two-argument `expint(ν, z)` loses roughly 8 bits, and up to 17, to cancellation in its series about the origin and in its continued fraction. At `Float64` this is invisible in most uses; at `Float32` it is not, and it shows up as the one- and two-argument forms disagreeing about the same value.

SpecialFunctions v2.8.3, Julia 1.12.7.

```julia
julia> mism = 0; worst = 0.0;
julia> for x in range(0.01f0, 20f0, length=2000)
a, b = expint(x), expint(1, x) # E₁(x), same value
a == b || (mism += 1)
end
julia> mism
1550 # of 2000 points, worst 439 ulps at x = 2.94
```

The one-argument entry evaluates `E₁` in `Float64` and rounds (`expint(x::Float32) = Float32(expint(Float64(x)))`), so it is accurate; the two-argument entry keeps `Float32` throughout and is not.

## Accuracy of the two-argument path

Correct bits against a 256-bit reference:

| ν | z | `Float32` | `Float64` |
|---|---|---|---|
| `1` | `1.5` | 19.8 / 24 | 49.9 / 53 |
| `1` | `2.5` | 16.4 / 24 | 53.1 / 53 |
| `1` | `2.94` | 15.1 / 24 | 52.5 / 53 |
| `2` | `2.94` | 16.3 / 24 | 47.9 / 53 |
| `5` | `2.94` | 14.8 / 24 | 44.2 / 53 |
| `1.001` | `2.94` | **7.1 / 24** | **36.2 / 53** |
| `5` | `5.0` | 16.5 / 24 | 45.5 / 53 |
| `2` | `8.0` | 21.4 / 24 | 46.1 / 53 |

The `Float64` column for `ν = 1` looks healthy only because `_expint` routes `ν == 1, real(z) > 0, z isa Union{Float64,Complex{Float64}}` to the tuned `expint_opt`; the series is what the other rows exercise.

## Mechanism

`En_expand_origin_posint` returns `gammaterm - sumterm`, and near the top of the small-|z| region those two nearly cancel. At `ν = 1`, `x = 2.94`:

```
|gammaterm| = 1.656
|sumterm| = 1.670
result = 0.01409
```

so the operands are ~118× the result and about `log2(118) ≈ 7` bits are lost structurally. The terms of `Σ (-z)^k/k!` alternate and peak near `k ≈ |z|`, so the loss grows with `|z|` and is largest just below the `abs2(z) < 9` switch to the continued fraction — which then has its own 5–7 bit loss. For `ν` just above a positive integer the near-pole cancellation stacks on top, which is the `ν = 1.001` row: 17 bits gone at `Float64`.

## Directions

* **Integer ν**: the upward recurrence `E_{n+1}(z) = (e^{-z} - z·E_n(z))/n` from the tuned `E₁` amplifies error by `|z|/n` per step, so it is stable exactly where the problem is (`n ≳ |z|`, `|z| < 3`).
* **Non-integer ν** needs an accurate base value, i.e. the incomplete gamma connection `E_ν(z) = z^(ν-1) Γ(1-ν, z)`. `gamma_inc`'s NSWC kernel is cancellation-free for `a > 0`, but `a = 1-ν < 0` for `ν > 1`, so it would need the reflection `Γ(a,x) = (Γ(a+1,x) - x^a e^{-x})/a` — which reintroduces a subtraction and needs its own error analysis.
* Widening the narrow types at the boundary (as the one-argument entry already does) would close the `Float32` discrepancy and hide the loss below `eps`, but does nothing for `Float64`.

## Note

#550 restructures these functions for precision handling but deliberately does not touch this; the numbers above are unchanged by it, except that `Float16` near a positive integer order no longer returns `0.0`.

---

*Filed by Claude Code on behalf of @andreasnoack; the investigation and the text above are Claude's.*

Start with `_expint`, `En_expand_origin_posint`, and `expint_opt`; trace the Float32 two-argument path and the `abs2(z) < 9` switch. Compare the integer-order recurrence and incomplete-gamma direction against the 256-bit reference data, then verify that one- and two-argument forms agree without the reported cancellation losses.

Written by the indexing model from the issue text.