Residue of $\Gamma^{2}$ and $\Gamma^{3}$ Based on wiki, the residues of $\Gamma$ at non positive integers are given by:
$$\text{Res}\left ( \Gamma(z),z=-n \right )=\frac{(-1)^{n}}{n!}.$$
I have been trying to find residue for $\Gamma^{2}$ and $\Gamma^{3}$ and could not find any. I tried deriving but not successful. How to find their residue? Or are they given?
 A: Fix $n$ and define
$$
\begin{align}
f(z)
&=\frac{\Gamma(z)}{\Gamma(z+n)}
=(-1)^n\frac{\Gamma(1-(z+n))}{\Gamma(1-z)}\tag{1}
\end{align}
$$
then
$$
\Gamma(z)=\frac{f(z)\Gamma(z+n+1)}{z+n}\tag{2}
$$

Using
$$
\frac{\Gamma'(z)}{\Gamma(z)}
=\psi(z)
=-\gamma+\sum_{k=0}^\infty\left(\frac1{k+1}-\frac1{k+z}\right)\tag{3}
$$
and
$$
\frac{f'(z)}{f(z)}
=\psi(1-z)-\psi(1-(z+n))\tag{4}
$$
we can get the series
$$
\begin{align}
\Gamma(z+n+1)
&=\Gamma(1)+\Gamma'(1)(z+n)+\frac{\Gamma''(1)}2(z+n)^2+O(z+n)^3\\
&=1+\psi(1)(z+n)+\frac{\psi'(1)+\psi(1)^2}2(z+n)^2+O(z+n)^3\tag{5}
\end{align}
$$
and
$$
\begin{align}
f(z)
&=f(-n)+f'(-n)(z+n)+\frac{f''(-n)}2(z+n)^2+O(z+n)^3\\
&=\frac{(-1)^n}{n!}\left(\vphantom{\frac12}\right.1+(\psi(n+1)-\psi(1))(z+n)\\
&\hphantom{=\frac{(-1)^n}{n!}\left(\vphantom{\frac12}\right.}+\frac{(\psi(n+1)-\psi(1))^2-(\psi'(n+1)-\psi'(1))}2(z+n)^2\\
&\hphantom{=\frac{(-1)^n}{n!}\left(\vphantom{\frac12}\right.}+O(z+n)^3\left.\vphantom{\frac12}\right)\tag{6}
\end{align}
$$

Multiplying $(5)$ and $(6)$ and dividing by $z+n$ gives $\Gamma(z)$ near $z=-n$ to be
$$
\frac{(-1)^n}{n!}\left(\frac1{z+n}+\psi(n+1)+{\small\frac{\pi^2+3\psi(n+1)^2-3\psi'(n+1)}6}(z+n)+O(z+n)^2\right)\tag{7}
$$
Looking at the coefficient of $\frac1{z+n}$ in $(7)$, we get
$$
\bbox[5px,border:2px solid #C0A000]{\operatorname*{Res}_{z=-n}\left(\Gamma(z)\right)=\frac{(-1)^n}{n!}}\tag{8}
$$
Looking at the coefficient of $\frac1{z+n}$ in the square of $(7)$, we get
$$
\bbox[5px,border:2px solid #C0A000]{\operatorname*{Res}_{z=-n}\left(\Gamma(z)^2\right)=\frac2{(n!)^2}\psi(n+1)}\tag{9}
$$
Looking at the coefficient of $\frac1{z+n}$ in the cube of $(7)$, we get
$$
\bbox[5px,border:2px solid #C0A000]{\operatorname*{Res}_{z=-n}\left(\Gamma(z)^3\right)=\frac{(-1)^n}{2(n!)^3}\left(\pi^2+9\psi(n+1)^2-3\psi'(n+1)\right)}\tag{10}
$$

To put the answers above in the same terms that Jack D'Aurizio uses, we have
$$
\psi(n+1)=H_n-\gamma\tag{11}
$$
and
$$
\psi'(n+1)=\frac{\pi^2}6-H_n^{(2)}\tag{12}
$$
A: Since $\Gamma(z)=\frac{1}{z}\Gamma(z+1)=\frac{1}{z(z+1)}\Gamma(z+2)=\frac{1}{z(z+1)(z+2)}\Gamma(z+3)=\ldots$ it follows that $z=-n$ is a double pole for $\Gamma(z)^2$, hence:
$$\text{Res}\left(\Gamma(z)^2,z=-n\right) = \left. \frac{d}{dz}(z+n)^2\Gamma(z)^2\right|_{z=-n}.$$
Since, from: 
$$\psi(z)=\frac{\Gamma'(z)}{\Gamma(z)}=-\gamma+\sum_{k=0}^\infty\left(\frac1{k+1}-\frac1{k+z}\right)$$
we have:
$$\Gamma(x+n+1)=1-\gamma(x+n)+O((x+n)^2)\tag{1}$$
while:
$$\frac{1}{x(x+1)\cdot\ldots\cdot(x+n)}=\frac{(-1)^n}{n!(x+n)}+A_n+O((x+n))\tag{2}$$
where:
$$ A_n = \left.\frac{d}{dx}\frac{1}{x(x+1)\cdot\ldots\cdot(x+n-1)}\right|_{x=-n}=\frac{(-1)^n H_n}{n!}.$$
By multiplying $(1)$ and $(2)$ it follows that:
$$ \Gamma(x) = \frac{(-1)^n}{n!(x+n)}+\frac{(-1)^n}{n!}\left(H_n-\gamma\right)+O((x+n))\tag{3}$$
and by squaring such Laurent series it follows that:
$$ \text{Res}\left(\Gamma(z)^2,z=-n\right)=\color{red}{\frac{2(H_n-\gamma)}{n!^2}}.\tag{4}$$
To compute the residues of $\Gamma(z)^3$, we have to compute one extra term in $(1)$ and $(2)$ then follow the same lines. That is tedious but not difficult:
$$ \text{Res}\left(\Gamma(z)^3,z=-n\right)=\color{red}{\frac{3(-1)^n\left(3\gamma^2+\zeta(2)-6\gamma H_n+3H_n^2+H_{n}^{(2)}\right)}{2\, n!^3}}.\tag{5}$$
A: For the case of $\Gamma^{4}(x)$ the residue is
\begin{align}
Res_{x \rightarrow -n}( \Gamma^{4}(x)) = \frac{1}{3!} \left( \frac{(-1)^{n}}{n!}\right)^{4} \left[ 8 H_{n,3} + 12 (H_{n} - \gamma) H_{n,2} + 36 (H_{n} - \gamma) \zeta(2) + 46 (H_{n} - \gamma)^{3} - 4 \zeta(3) \right]
\end{align}
where $\gamma$ is the Euler-Mascheroni constant, $H_{n}$ are the Harmonic numbers, and 
\begin{align}
H_{n,r} = \sum_{k=1}^{n} \frac{1}{k^{r}}.
\end{align}
