expint truncates its series after a number of terms fixed for Float64, silently capping BigFloat accuracy

`expint` (and everything that routes through it) truncates its series and continued fractions after a number of terms that is fixed for `Float64`, while the methods accept any number type. For `BigFloat` the result therefore silently stops improving well before the requested precision is reached. There are two independent limits, plus a third place where a fixed-order expansion is (correctly) restricted but leaves `BigFloat` without a substitute.

All numbers below are SpecialFunctions v2.8.3 on Julia 1.12.7. The reference values come from MPFR's `expinti` via `E₁(x) = -Ei(-x)`, `E₁(-x + 0im) = -Ei(x) - πi` and `E₂(z) = exp(-z) - z E₁(z)`, so they are independent of the code under test.

## 1. `niter = 1000` caps the continued fraction at ~92 digits

```julia
julia> using SpecialFunctions

julia> setprecision(BigFloat, 512) do # eps(BigFloat) = 1.5e-154
x = BigFloat(3)
ref = exp(-x) + x*expinti(-x) # = E₂(3)
(default = Float64(abs(expint(2, x) - ref)/abs(ref)),
niter_1e5 = Float64(abs(expint(2, x, 100_000) - ref)/abs(ref)))
end
(default = 1.1031886851560754e-92, niter_1e5 = 1.3250352104469716e-152)
```

`expint`/`expintx`/`_expint` take `niter::Int = 1000`. The continued fraction is slowest just above the cutoff `abs(z) == 3`, where it gains only about `10*sqrt(n)` bits from `n` terms:

| `niter` | 1000 | 2000 | 4096 | 8000 | 16000 |
|---|---|---|---|---|---|
| bits of `expint(2, big(3))` correct | 305 | 436 | 629 | 883 | ≥1024 |

so the number of terms needed grows *quadratically* with the precision: `En_cf(big(2), big(3), 400_000)` reports convergence after 10562 iterations at 1024 bits and 168194 at 4096 bits. 1000 terms are ample for `Float64` but cap every higher-precision result at roughly 92 digits.

## 2. `En_taylor` is a hard 101-term series, capping the negative real axis at ~97 digits

```julia
julia> setprecision(BigFloat, 512) do
y = BigFloat(4)
ref = exp(y) + y*(-expinti(y) - im*BigFloat(π)) # = E₂(-4 + 0im)
z = Complex{BigFloat}(-y, 0)
(default = Float64(abs(expint(2, z) - ref)/abs(ref)),
niter_1e5 = Float64(abs(expint(2, z, 100_000) - ref)/abs(ref)))
end
(default = 5.2022497324938007e-98, niter_1e5 = 8.768948249862689e-101)
```

Note that raising `niter` does *not* help here: the limit is the

```julia
for k = 0:100
```

loop in `En_taylor`, which is used by the `real(z) < 0` procedure to step towards the axis. It has a proper `ϵ = 10*eps(real(asum))` convergence check but a term count fixed for `Float64`. Raising just that bound to 1000 takes the same computation from 323 to 1017 correct bits at 1024-bit precision, at no extra cost in time.

## 3. The near-pole δ-expansion has no high-precision counterpart

In `En_expand_origin_general` the correction used when `ν` is close to a positive integer is a fixed 5-term expansion in `δ = round(ν) - ν` (orders `δ⁰` … `δ⁴`), and it is guarded by

```julia
if real(ν+z) isa Union{Float64, Float32} && abs(gammaterm - blowup) < 1e-3 * abs(blowup)
```

That guard is the right call — with only five terms the truncation error `O(δ⁵)` is not below `eps(BigFloat)` — but it means `BigFloat` falls back to the cancelling `gammaterm - (blowup + sumterm)` form and loses about `log2(1/δ)` bits:

```
expint(2 + 1e-20, big(0.5)) at 128-bit precision: 61.8 of 128 bits correct
expint(2 + 1e-40, big(0.5)) at 512-bit precision: 382.6 of 512 bits correct
```

Extending the series is unattractive (the coefficients are polynomials in `polygamma` values). For `BigFloat` the natural fix is guard bits instead: evaluate the plain expression with the working precision raised by `ceil(log2(1/abs(δ)))`. There is precedent for that in the package — `loggamma(z::Complex{BigFloat})` already guardrails its precision.

## Suggested fixes for 1 and 2

Since every loop stops as soon as its `eps`-based criterion is met, a generous limit costs nothing for arguments that converge quickly. Making the two limits scale with the precision of the arguments is enough:

```julia
# quadratic, because the continued fraction gains ~10*sqrt(n) bits from n terms
_expint_niter(ν, z) = max(1000, div(precision(float(real(promote_type(typeof(ν), typeof(z))))), 8)^2)

# linear is ample here, the series converges factorially
for k = 0:max(100, precision(real(asum))) # En_taylor
```

`precision` is 53, 24 and 11 for `Float64`, `Float32` and `Float16`, so all three defaults evaluate to the current literals and hardware-float results are unchanged (I checked 144 values of `expint`/`expintx`/`gamma(a,x)` across real/complex `Float16/32/64` orders and arguments: bit-identical). With both limits scaled, the three probes above reach ~`100*eps(BigFloat)` at 512 bits, and the negative-axis case goes from 323 to 2045 correct bits at 2048-bit precision. It is also *faster* at 1024 bits (13 s versus 20 s), because the old code compensated for the small `quick_niter = niter >> 4` by starting much further off the axis and taking more Taylor steps. That path remains expensive at high precision, though: about 2 s at 512 bits and 190 s at 2048 bits for `expint(2, Complex{BigFloat}(-4))`.

## Related: non-convergence is silent

`En_cf` returns the iteration count it used, and `_expint` discards it (`g, cf, _ = En_cf(ν, z, niter)`). Once the limit tracks the requested precision, reaching it means genuine non-convergence rather than "`Float64` was already enough", so it would be a good place to warn or throw instead of returning a quietly inaccurate result.

For context, I went looking for this pattern across the package and `expint` is the only place where it produces silently inaccurate results. `_ellipk`/`_ellipe` (`Float64`-only), `sinint`/`cosint` (explicit `error` for other `AbstractFloat`s), `_jinc` (Taylor shortcut restricted to `Union{Float32,Float64}` with per-type thresholds), `@E₁_cf64`/`@E₁_taylor64` (gated on `Float64`/`ComplexF64`), and the `beta_inc*`, `_gamma_inc`, `_gamma_inc_inv`, `erfi`/`dawson`/`faddeeva` internals all either restrict their signatures or route `BigFloat` to MPFR.

---

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

Start by tracing _expint and En_cf for the continued-fraction limit, then inspect En_taylor and En_expand_origin_general for the other precision-dependent paths. Compare the existing BigFloat guardrails in loggamma and run the reported high-precision expint probes. Done means BigFloat accuracy is no longer capped by Float64-sized limits or cancellation, with non-convergence handled explicitly.

Written by the indexing model from the issue text.