Function polygamma_impl 

pub fn polygamma_impl<P, E, V>(x_in: V, n: u32) -> V
where
    P: Policy,
    E: FloatElementWithBits + BernoulliNumbers + CotPiDerivatives + Factorials,
    V: FloatVectorWithBits<Element = E> + SpecializedSpecialMath<E>,

Shared polygamma (psi_n) implementation for all real element types.

$\psi_n(x) = \frac{\mathrm{d}^n}{\mathrm{d}x^n}\psi(x)$, the (n+1)-th derivative of $\ln\Gamma$. n = 0 and n = 1 delegate to the tuned digamma and trigamma kernels. n >= 2 runs the two-region scheme Boost.Math’s polygamma_imp uses on the positive axis, restructured for vectors:

• a masked forward recurrence $\psi_n(x) = \psi_n(x+1) + (-1)^{n-1} n!\,x^{-(n+1)}$ walks every lane up to the transition point $N = 0.4\,d_{10} + 4n$ (with $d_{10}$ the format’s decimal digits), then
• the asymptotic expansion at large x,

\psi_n(x) = (-1)^{n-1}\left[\frac{(n-1)!}{x^n} + \frac{n!}{2x^{n+1}}
    + \sum_{k\ge1} B_{2k}\,\frac{(2k+n-1)!}{(2k)!\,x^{2k+n}}\right]

evaluated by the term-ratio recurrence, so the scalar order-dependent coefficients $(n+2k-2)(n+2k-1)/((2k-1)\,2k)$ are splatted and the vector work per term is one multiply by $1/x^2$. The $B_{2k}$ come from BernoulliNumbers::B2N, and the series converges well before that table ends for every x past the transition point.

Unlike Boost there is no separate near-zero zeta series: the leading $n!/x^{n+1}$ term dominates so completely below the recurrence range that the walk loses nothing.

Negative arguments reflect through $\psi_n(x) = (-1)^n\left[\psi_n(1-x) + \pi\,\frac{\mathrm{d}^n}{\mathrm{d}z^n}\cot(\pi z)\big|_{z=1-x}\right]$, with the cot derivative’s cosine polynomial from CotPiDerivatives::COT_PI_ROWS and both sin_pi/cos_pi evaluated at x itself, the smaller-magnitude representative (they agree with the 1 - x values exactly, by periodicity, but carry less argument error). At the poles (zero and the negative integers) odd n yields +inf, the correct two-sided limit. Even n has one-sided limits of opposite sign and yields NaN when overflow checking is enabled.

Current limits (deliberate, documented rather than patched)

• Reflection stops at n = 20, the cot-pi table’s reach (Boost tabulates the same range, with its runtime coefficient recurrence past it queued work). Negative arguments at n > 20 return NaN. The positive axis is unaffected.
• Direct powers bound the domain. x^(n+1) is formed directly at arguments up to max(x, N), so lanes where (n+1) log10(max(x, N)) exceeds the format’s decimal exponent range (~300 for f64, ~36 for f32) flush to zero even where psi_n itself is representable (e.g. psi_100(1e4) ~ -9.4e-245), and n! likewise overflows at n >= 171 (f64) / n >= 35 (f32). Boost rescues both with log-domain arithmetic. That is queued work, and the tests pin the boundary.