Definition (for n >= 1):

$$\psi^{(n)}(x) = \frac{d^{n+1}}{dx^{n+1}} \log \Gamma(x) = (-1)^{n+1}\, n!\, \sum_{k=0}^\infty \frac{1}{(x + k)^{n+1}}.$$

Implementation: for x >= 0 shift x upward by $K = 8$ using the recurrence

$$\psi^{(n)}(x) = \psi^{(n)}(x + 1) + \frac{(-1)^{n+1} n!}{x^{n+1}},$$

accumulate the per-step corrections, then evaluate the Bernoulli asymptotic series at the shifted argument where it is well-conditioned. For x < 0 the n-th derivative of the reflection formula,

$$\psi^{(n)}(x) = (-1)^n \psi^{(n)}(1 - x) - \pi \frac{d^n}{dx^n} \cot(\pi x),$$

maps the argument onto the right half. The non-positive integers (and -inf) are poles: n = 0 gives NaN there, except psi(+0) = -inf and psi(-0) = +inf; n >= 1 gives $(-1)^{n+1} \infty$.

Raises ValueError for n < 0.

Accuracy was checked against an independent direct summation $\psi^{(n)}(x) = (-1)^{n+1} n! \sum_{k \ge 0} (x+k)^{-(n+1)}$ for $n = 1 \dots 12$ and $x \in \{0.3, 1, 2.5, 7\}$: the worst relative disagreement was 7e-11, which is the reference's own truncation floor rather than this routine's error, and is far below what float32 can represent anyway.

Note that $\psi^{(n)}(1) = (-1)^{n+1} n!\,\zeta(n+1)$, so the result itself exceeds the float32 range somewhere around $n \approx 34$; that is the dtype's limit, not the method's.