Residues of dual-pole Γ functions

Question

In working with a certain Meijer G function, I encounter an integral that requires I compute the following residue: $$\mathrm{Res}\left(f(s); s = a+n\right)$$ for $$ f(s)= Γ(a-s) Γ(a - 1 - s) Γ(b - s) Γ(c - s) z^s,$$ for $a$, $b$, and $c$ rational; $z>0$; and $n \geq 0$ is a non-negative integer. The idea is that there's always exactly one dual-pole ($Γ(a-s) Γ(a - 1 - s)$ will both blow up/down as $s → a+n$ because they get an argument of 0 or a negative integer) and a few simple poles—in this example, two simple poles, $b$, and $c$, though in my applications there can be more. I can deal with the residue at the simple poles, my difficulty is with the residue at the double-poles, at $a+n$.

A lot of the discussion about residues online involves rational functions with polynomial denominators, which I obviously don't have. Wikipedia suggests a limit formula for higher-order poles according to which the above residue is $$ \lim_{s \rightarrow a+n} \frac{\partial}{\partial s} (s-(a+n))^2 f(s) $$ but this can't work because the derivative is $$(s-(a+n)) \cdot f(s) \cdot \Big( 2 + (a+n-s) \cdot \big( \psi(a-s) + \psi(a-1-s) + \psi(b-s) + \psi(c-s) - \log(z) \big) \Big),$$ and because each of the $\psi$ terms (the polygamma function) that are a function of $a$ involve

1. $(s-(a+n))^2$ both of which go to zero in the limit as $s \rightarrow a+n$,
2. $\Gamma(a-s) \Gamma(a-1-s)$ in $f(s)$, both of which go to $\pm \infty$, and finally
3. $\psi(a-s)$ or $\psi(a-1-s)$ as a "tie-breaker" also goes to $\pm \infty$,

so the limit seeeeeems to not converge?

However. I can use Sympy to derive the residues for various $n=0,1,...$ and the residues it gives me seem correct (I can cross-check them by (1) evaluating the Meijer G function via hypergeometric functions and also (2) numerical integration of a probabilistic expectation that yields that Meijer G expression and (3) Monte Carlo). I provide Sympy's residue for a few $n$ below, where $\gamma$ is the Euler-Mascheroni constant:

n=0 $$z^{a} \big(- \log{\left(z \right)} + \psi{\left(- a + b \right)} + \psi{\left(- a + c \right)} - 2 \gamma + 1\big) \Gamma\left(- a + b\right) \Gamma\left(- a + c\right)$$

n=1 $$\frac{1}{4}z^{a + 1} \big(- 2 \log{\left(z \right)} + 2 \psi{\left(- a + b - 1 \right)} + 2 \psi{\left(- a + c - 1 \right)} - 4 \gamma + 5\big) \Gamma\left(- a + b - 1\right) \Gamma\left(- a + c - 1\right)$$

n=2 $$\frac{1}{36} z^{a + 2} \big(- 3 \log{\left(z \right)} + 3 \psi{\left(- a + b - 2 \right)} + 3 \psi{\left(- a + c - 2 \right)} - 6 \gamma + 10\big) \Gamma\left(- a + b - 2\right) \Gamma\left(- a + c - 2\right)$$

n=3 $$\frac{1}{1728}z^{a + 3} \big(- 12 \log{\left(z \right)} + 12 \psi{\left(- a + b - 3 \right)} + 12 \psi{\left(- a + c - 3 \right)} - 24 \gamma + 47\big) \Gamma\left(- a + b - 3\right) \Gamma\left(- a + c - 3\right)$$

I can soooort of relate these Sympy residues to the limit formula above, but I totally can't seem to find the pattern in the denominators above, e.g., 1, 4, 36, 1728, …

Where am I going wrong?

---

(Python code follows)

import sympy as sp
s, z, a, b, c, d, e = sp.symbols('s z a b c d e')
f = sp.gamma(a - s) * sp.gamma(a - 1 - s) * sp.gamma(b - s) * sp.gamma(c - s) * z**s
for added in range(4):
  print(f'\nn={added}')
  res = sp.residue(f, s, a + added)
printable = sp.simplify(res.subs({sp.exp(a * sp.log(z)): z**a}))
  sp.pprint(printable)
for Math Exchange:
print('$$' + sp.latex(printable).replace("\operatorname{polygamma}", "\psi") + '$$')

Answer 1

Mess with the orders of zeros/poles. Let's redo the computations: \begin{align*} \operatorname*{Res}_{s=a+n}f(s)&=\lim_{s\to a+n}\big((s-a-n)^2f(s)\big)'=L_1\cdot L_2, \\L_1&=\lim_{s\to a+n}(s-a-n)^2f(s), \\L_2&=\lim_{s\to a+n}\left(\frac2{s-a-n}+\frac{f'(s)}{f(s)}\right). \end{align*} Now $L_1$ equals $\Gamma(b-a-n)\Gamma(c-a-n)z^{a+n}\cdot L_{11}\cdot L_{12}$, where \begin{align*} L_{11}&=\lim_{s\to a+n}(s-a-n)\Gamma(a-s)=\lim_{t\to 0}t\,\Gamma(-n-t)=\frac{(-1)^{n+1}}{n!}, \\L_{12}&=\lim_{s\to a+n}(s-a-n)\Gamma(a-1-s)=\dots=\frac{(-1)^n}{(n+1)!}. \end{align*} Similarly, $L_2$ equals $\log z-\psi(b-a-n)-\psi(c-a-n)+L_{21}+L_{22}$, where \begin{align*} L_{21}&=\lim_{s\to a+n}\left(\frac1{s-a-n}-\psi(a-s)\right)=\dots=\gamma-H_n, \\L_{22}&=\lim_{s\to a+n}\left(\frac1{s-a-n}-\psi(a-1-s)\right)=\gamma-H_{n+1}. \end{align*}