Function exp_sum 

pub fn exp_sum<P, E, V>(base: V, rest_hi: V, rest_lo: V) -> V
where
    P: Policy,
    E: FloatElement,
    V: FloatVector<Element = E> + SpecializedSpecialMath<E>,

exp(base + rest_hi + rest_lo) where base is a large exact-ish number (-lambda, x/4) and the rest pair is small: TwoSum recovers the rounding of the sum, and e^{s + lo} = e^s (1 + lo) to first order. Without it the sum rounds to half an ulp of base, which the exponential turns into hundreds of ulp.

rest_lo is the caller’s own second word, added to the residual this function already recovers. It exists because the residual alone is not enough: base’s rounding is only half the problem, and rest arrives from pmf_parts carrying an error of its own that no amount of care in this sum can recover. Measured 2026-09-07 on poisson_pmf(0, lambda), which is exactly e^-lambda: 14.39 ulp with rest_lo dropped, 0.33 with it. The TwoSum itself is SpecializedSpecialMath::exp_two_sum, not spelled here: an error-free transformation written as + and - is only error-free when those are strict, and on the scalar backend under algebraic-scalar they are not. The trait method’s default therefore returns a zero residual and this degrades cleanly to exp(base + rest_hi + rest_lo); ps/pd and Dual override it to recover the real thing. See that method’s docs.