pub fn expint_double_primal_n<P: Policy, E, V, const N: usize>(x: V) -> (V, V)where
E: FloatElementWithBits + ExpIntConsts<N>,
V: FloatVectorWithBits<Element = E> + SpecializedSpecialMath<E>,Expand description
$E_N(x)$ together with the adjacent lower order $E_{N-1}(x)$, which is
$-E_N'(x)$ by differentiation under the integral sign.
The lower order comes from whichever direction is stable in the regime the value
itself was computed in: below ExpIntConsts::RECURRENCE_THRESHOLD the forward
recurrence is running anyway, so E_{N-1} is just its previous iterate. Above it,
where the asymptotic series takes over, the recurrence is inverted instead –
$E_{N-1}(x) = (e^{-x} - (N-1) E_N(x)) / x$. Inverting is the stable direction
(it damps by 1/x where the forward one amplifies by x) and its only weakness,
the cancellation as x -> 0, is unreachable here because that branch only runs for
very large x.