Function inv_digamma_impl 

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

The x > 0 with digamma(x) = y.

Seed (Minka, “Estimating a Dirichlet distribution”, appendix): e^y + 1/2 for y >= -2.22 and -1/(y - digamma(1)) below, both from the asymptotics at the two ends and within a factor of two of the root everywhere. Newton with trigamma from there, bracketed by that factor of two on each side.

Above y = 3 (x > 20) the answer is the Stirling series solved for x instead. digamma(x) = ln x - t(x) with t = 1/(2x) + 1/(12x^2) - 1/(120x^4) + ... - 691/(32760 x^12) (5e-18 at x = 20), so x = e^y e^{t(x)}. Writing x = e^y u, the unknown u = e^{t} is the root of h(u) = u - e^{t(e^y u)}, which Newton takes quadratically from u_0 = e^{t(e^y + 1/2)} in two steps. Nothing in h subtracts y, so the result is within an ulp or two of e^y’s own rounding, where Newton on digamma is bounded by eps * y relative (the residual digamma(x) - y is ln x - y to a rounding of eps * y). The plain fixed point x <- e^y e^{t(x)} contracts only by 1/(2x) per pass and would need seven passes at x = 20.