I'm writing a code with the aim of automating operator mixing calculations in Euclidean QCD. I have alreaady written a module for generating all graphs and another for generating and assigning Feynman rules. Eventually I would like to be able to obtain exact solutions for generic parameters. This would involve making no assumptions about the relevant kinematic variables. To that end, I am compiling a list of so-called "master integrals" of the form
$$I(a,b;m_1,m_2)=\int\frac{d^dk}{(2\pi)^d}\frac{1}{(k^2+m_1^2)^a((p+k)^2+m_2^2)^b}$$
for the two-point functions. (I haven't even gotten to the three-point or higher integrals yet). While these integrals should in principle be expressible in terms of Appell functions, I want to start small with the simple cases $(m_1,m_2) = (0,0)$, $(m,0)$, and $(m,m)$. The first two cases are very easy. I started encountering problems with the third case, $I(a,b;m,m)$. The most direct method I can think of for this integral is the Mellin-Barnes representation,
$$\frac{1}{(k^2+m^2)^a}=\frac{1}{2\pi i}\int_\mathcal{L}dz\frac{(m^2)^z}{(k^2)^{z+a}}B(z+a,-z)$$
(where the countor is along the entire imaginary axis, deformed to separate the poles of $\Gamma(a+z)$ and $\Gamma(-z)$), but there are of course two solutions corresponding to the choice of closing the contour on the left or the right half-plane. In the semifinal step, I am left with an integral of the form
$$\int_\mathcal{L}dz\ B(a+z,b+z)\Gamma(a+b-d/2+z)\Gamma(-z)R^{-z},$$
where $R=m^2/p^2$. I can easily calculate this integral; closing the contour on the right gives a single ${}_3F_2$, while closing it on the left give a triple sum of ${}_3F_2$s, both as expected according to the pole structure. Here is, however, where I feel that I've been playing fast and loose in defining $\mathcal{L}$. The two results I name are only valid (I think) for $a\neq b$, $a+b\neq d/2$, $\textit{etc.}$; that is to say, they are only valid for those cases where $a$ and $b$ have values such that all poles are simple.
As a sanity check, I compute the one-loop self energy for gauge fields and confirm that the left and right semicircles are inequivalent for the diagram containing a fermion loop. Only the left contour works in this case, but of course an automated calculation should not require a manual selection of contour.
Here is my question. Is there a simple(-ish) way to cover all of the cases generically? I'd imagine that I will have to define $I(a,b;m,m)$ in a piecewise fashion over the space of $(a,b)$ values, but I'm having trouble determining whether or not I'm covering all my bases.
Thank you in advance for the help!