Function digamma_impl 

pub fn digamma_impl<P, E, V, const NR: usize, const NL: usize, const NP: usize, const NQ: usize>(
    x_in: V,
    t: &Digamma<E, NR, NL, NP, NQ>,
) -> V
where
    P: Policy,
    E: FloatElementWithBits,
    V: FloatVectorWithBits<Element = E> + SpecializedSpecialMath<E>,

Shared digamma (psi) implementation for all real element types.

psi(x) = d/dx ln(Gamma(x)). The element-specific rational/asymptotic coefficients are passed in so the f32 and f64 specializations can share this body:

• y / roots / p_12 / q_12: the [1, 2] rational psi(x) = (x - root)(Y + R(x-1)), where root is summed from roots via staged subtraction to preserve bits.
• p_large: the x >= 10 asymptotic expansion in 1/(x-1)^2.