Function phi_internal_n 

pub fn phi_internal_n<V, E, P, const N: usize, const ADAPTIVE: bool>(
    z: V,
    terms: usize,
) -> V
where
    E: FloatElement,
    V: FloatVector<Element = E> + SpecializedTranscendentalMath<E>,
    P: Policy,

phi_N(z) = sum_{n>=0} z^n/(n+N)!, the exponential-integrator functions.

N = 0 is exp and N = 1 is expm1(z)/z. Beyond that, two arms split at |z| = N:

• Below, the series, summed forward from 1/N! with each term the previous times z/(N+k). Every coefficient is one small ratio, so this stays exact for any element type. terms bounds the loop, and with ADAPTIVE it also stops as soon as the term it just added is under half an ulp of the sum, which is how an element whose precision is not known statically (Compensated) converges to its own epsilon.
• Above, the recurrence phi_{k+1} = (phi_k - 1/k!)/z upward from phi_1 = expm1(z)/z. Each step subtracts a constant from something that is only just larger than it while |z| is small (that is the cancellation the series exists to avoid), but the amplification per step is phi_k/(phi_k - 1/k!), which is bounded once |z| >= k. |z| >= N covers every step, and measured against mpmath the recurrence stays under 6 ulp for N <= 8 in both f32 and f64. The same bound is why the series arm stops at N: its terms are monotone there, so the alternating negative side does not cancel either.

Both arms overflow gracefully: expm1 saturates to +inf and each division by z leaves it there, and -inf gives -1 * -0 and then a run of +0s, the limit. Only +inf itself, inf * (1/inf), and the 0/0 of N = 1 at the origin need patching.