Function pmf_parts 

pub fn pmf_parts<P, E, V, const ALL_LARGE: bool>(
    k: V,
    lambda: V,
) -> (V, V, V::Mask, V, V)
where
    P: Policy,
    E: FloatElement,
    V: FloatVector<Element = E> + SpecializedSpecialMath<E>,

The shared core of every $x^k e^{-x}/\Gamma(k+1)$ shape here: the Poisson mass, its log, and the Laguerre-function seed. Returns (rest_hi, rest_lo, large, prod, n) such that

P(k; lambda) = exp(rest_hi + rest_lo - [large ? 0 : lambda]) * prod / sqrt(2 pi n)

rest_lo is the second word of the exponent, nonzero only on the shifted lanes (see the split below) and zero for k >= 9 and the peak series. Callers must pass it to [exp_two_sum]; dropping it cost poisson_pmf 14 ulp.

where n, prod are from [shift_to_stirling] (n = k, prod = 1 for k >= 9) and rest is one of three things, per lane, always with -stirlerr(n):

• k >= 9 near the peak (|k - lambda| < 0.2 (k + lambda)): -bd0(k, lambda) as its series, with nothing large in it.
• k >= 9 off the peak: -bd0 directly as -(k ln(k/lambda) - (k - lambda)). The ratio goes into the ln whole and k - lambda is exact when they are close, so this is a few ulp too, unlike k (ln k - ln lambda) (40x worse when k ~ lambda) or splitting -lambda off (k ln(k/lambda) + k is then large on its own, and both were measured at 80-150 ulp for k = 100, lambda = 150).
• k < 9: k ln lambda - n ln n + n, the shifted Stirling form, with - lambda left out so the caller adds it with a TwoSum, it being the one large term there.

The last two share one ln, of k / lambda or n by lane. ln lambda is only formed if some lane is small. Under ALL_LARGE (a caller who knows k >= 9 everywhere) the shift and ln lambda fold away. Uniform vectors skip whichever branch no lane needs. k = 0 gives 0 * ln 0 = NaN at lambda = 0; callers pin that.