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?
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} $$
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}$$
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}
