First of all, I know this is a long and intimidating question, but I'm sure anyone with a little bit of complex analysis experience can answer it, so please don't turn away.
In this paper by Flajolet and Salvy, they use the series representation of the Digamma function and the fact that, if $p$ is a nonnegative integer, as $s\to p$, $$\psi(-s)+\gamma\sim \frac{1}{s-p}+H_p+...$$ to state that as $s\to p$, $$(\psi(-s)+\gamma)^2\sim \frac{1}{(s-p)^2}+\frac{2H_p}{s-p}+...$$ and so if $r(s)$ is a function with no poles or zeroes at the positive integers, $$\begin{align} \text{Res}\space\space (\psi(-s)+\gamma)^2 r(s)_{s=p} &= \text{Res}\space\space \frac{r(s)}{(s-p)^2}+\frac{2H_p\space r(s)}{s-p}_{s=p}\\ &= \text{Res}\space\space \frac{r(s)}{(s-p)^2}+\frac{2H_p\space r(s)}{s-p}\\ &= \lim_{s\to p} \frac{d}{ds}r(s)_{s=p}+\lim_{s\to p}2H_pr(s)\\ &= r'(p)+2H_pr(p)\\ \end{align}$$ As long as $r(s)$ is $\mathcal O(s^{-2})$ at infinity, the integral $\oint (\psi(-s)+\gamma)^2 r(s) ds$ around a circle with radius approaching infinity should go to zero, meaning that the sum of all residues of $(\psi(-s)+\gamma)^2 r(s)$ should equal zero. Thus, if $q_i$ are the poles of $r(s)$, we have $$\color{green}{\sum_{p=0}^\infty (r'(p)+2H_p r(p))+\sum_{q_i}\text{Res} (\psi(-s)+\gamma)^2 r(s)_{s=q_i}=0}$$ Which is the result given in the paper. However, I tried to extend this further by saying that near its poles, $$(\psi(-s)+\gamma)^3\sim \frac{1}{(s-p)^3}+\frac{3H_p}{(s-p)^2}+\frac{3H_p^2}{s-p}...$$ Which should yield $$\text{Res}\space\space (\psi(-s)+\gamma)^3 r(s)_{s=p}=\frac{1}{2}r''(p)+3H_p r'(p)+3H_p^2 r(p)$$ and the summatory formula $$\color{red}{\sum_{p=0}^\infty \bigg(\frac{1}{2}r''(p)+3H_p r'(p)+3H_p^2 r(p)\bigg)+\sum_{q_i}\text{Res} (\psi(-s)+\gamma)^3 r(s)_{s=q_i}=0}$$ But this isn't what the paper says. The paper says that $$\color{green}{\sum_{p=0}^\infty \bigg(\frac{1}{2}r''(p)-3r(p)\zeta(2)+3H_p r'(p)+3H_p^2 r(p)-3H_p^{(2)}r(p)\bigg)+\sum_{q_i}\text{Res} (\psi(-s)+\gamma)^3 r(s)_{s=q_i}=0}$$
Why is this so? Where was my mistake?
When $f$ has a simple pole - we assume that the residue is $1$ for simplicity of notation, the general case follows by multiplication with a constant - at $p$, write
$$f(s) = \frac{1}{s-p} + g(s)$$
with $g$ holomorphic on a neighbourhood of $p$. Then by the binomial theorem
$$f(s)^m = \sum_{k = 0}^{m} \binom{m}{k} \frac{g(s)^k}{(s-p)^{m-k}}\,,\tag{$\ast$}$$
and to compute the residue of $f(s)^m r(s)$ for a holomorphic $r$ with a zero of order $0 \leqslant\mu < m$ at $p$, we need the terms for $k \leqslant m - \mu - 1$ on the right hand side of $(\ast)$, since all these may contain a term $\frac{c}{s-p}$ with $c\neq 0$ in their Laurent expansion about $p$. Writing $r(s) = (s-p)^{\mu} u(s)$, we obtain
$$f(s)^m\cdot r(s) = \sum_{k = 0}^{m - \mu - 1} \binom{m}{k} \frac{g(s)^k u(s)}{(s-p)^{m-\mu-k}} + \sum_{k = m-\mu}^{m} \binom{m}{k}(s-p)^{k+\mu-m}g(s)^ku(s)$$
where the second sum is holomorphic at $p$, and hence \begin{align} \operatorname{Res} \bigl(f^m(s) r(s)\bigr)_{s = p} &= \sum_{k = 0}^{m - \mu - 1} \binom{m}{k} \operatorname{Res}\biggl(\frac{g(s)^ku(s)}{(s-p)^{m-\mu-k}}\biggr)_{s = p} \\ &= \sum_{k = 0}^{m - \mu - 1} \binom{m}{k} \frac{1}{(m-\mu-k-1)!} \biggl(\frac{d}{ds}\biggr)^{m-\mu-k-1}\bigl(g(s)^ku(s)\bigr)\biggr\rvert_{s = p}\,. \end{align}
For $m = 3$ and $\mu = 0$ that is
$$\operatorname{Res}\bigl(f(s)^3r(s)\bigr)_{s = p} = \frac{1}{2}r''(p) + 3\bigl(g'(p)r(p) + g(p)r'(p)\bigr) + 3g(p)^2r(p)\,.$$