Function bd0_series 

pub fn bd0_series<P, E, V>(k: V, diff: V, v: V) -> V
where
    P: Policy,
    E: FloatElement,
    V: FloatVector<Element = E> + SpecializedSpecialMath<E>,

The binomial/Poisson deviance bd0(k, lambda) = k ln(k/lambda) + lambda - k >= 0 in its peak form: (k - lambda) v + 2k sum_{j>=1} v^{2j+1}/(2j+1) for v = (k - lambda)/(k + lambda) (from ln((1+v)/(1-v)) = 2 atanh v), bd0_terms terms, full precision inside |v| < 1/5. Away from the peak the direct form is fine and pmf_parts uses it (folded with the rest of the exponent); this is only the part that needs care. diff and v are passed in because callers have them.