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: shift x upward by $K = 6$ 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. Accuracy is roughly seven decimal digits across the positive-real regime.

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.