Function wright_omega_impl 

pub fn wright_omega_impl<P, E, V>(x: V) -> V
where
    P: Policy,
    E: FloatElement,
    V: FloatVector<Element = E> + SpecializedSpecialMath<E>,

The Wright omega function, the w > 0 with w + ln w = x, which is W_0(e^x) without ever forming e^x.

Seeds by region, after Lawrence, Corless and Jeffrey (2012, the algorithm SciPy uses): q(1 - q(1 - 3q/2)) in q = e^x for x <= -2, the series about x = 1 for -2 < x < 1, and x - ln x + ln x / x above, each one tier down. Newton on w + ln w - x with 1 + 1/w from there, bracketed by a factor of two.

Below x = -7 the residual w + ln w - x is the difference of two numbers near x whose true difference is e^x < 1e-3, so Newton is bounded by eps * |x| relative. The Lagrange series w = sum_{n>=1} (-n)^{n-1}/n! q^n to six terms is 2e-17 relative there and is the whole answer, one exp and five FMAs. Above x = 1e20 the seed is the answer to working precision and Newton’s first residual is already within tolerance.