Let ${\bf A} \in M_{n\times n}(\mathbb R)$ be a symmetric positive-definite matrix, ${\bf x}\in \mathbb R^n$. I know how to obtain (for $p>n/2$): $$ \int_{\mathbb R^n} \frac{{\rm d}{\bf x}}{({\bf x}^T {\bf A}{\bf x} + 1)^p} = \frac{1}{\sqrt{\det{\bf A}}} \frac{\pi^{n/2}\Gamma(p-\frac{n}{2})}{\Gamma(p)} $$ by moving to a basis that diagonalises matrix ${\bf A}$. I seek to calculate $$ I = \int_{\mathbb R^n} \frac{{\rm d}{\bf x}}{({\bf x}^T {\bf A}{\bf x} + 1)^{p}({\bf x}^T {\bf B}{\bf x} + 1)^{q}}$$ where both ${\bf A}$ and ${\bf B}$ are symmetric positive-definite matrices. The same method won't work, as in general there's no way to simultaneously diagonalize ${\bf A}$ and ${\bf B}$. It can only transform the integral to the form $$ I =\frac{1}{\sqrt{\det{\bf A}}}\int_{\mathbb R^n} \frac{{\rm d}{\bf x}}{({\bf x}^T {\bf x} + 1)^{p}({\bf x}^T {\bf C}{\bf x} + 1)^{q}}$$ where $ {\bf C} = {\bf A}^{-\frac12}{\bf B}{\bf A}^{-\frac12}$, but I don't know if anything more can be done.
EDIT: After searching around a bit more, I've found formula $$ \frac{1}{a^p b^q} = \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)}\int_0^1 \frac{t^{p-1}(1-t)^{q-1}}{\big(a t + b(1-t)\big)^{p+q}}{\rm d}t$$ which allows to write $I$ in the form $$ I = \frac{1}{\sqrt{\det{\bf A}}}\frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} \int_{\mathbb R^n} \int_0^1 \frac{t^{p-1}(1-t)^{q-1}}{\big({\bf x}^T \big(t{\bf 1}+(1-t){\bf C}\big){\bf x} + 1\big)^{p+q}} {\rm d}t {\rm d}{\bf x} = \\ = \frac{1}{\sqrt{\det{\bf A}}} \frac{\pi^{n/2}\Gamma(p+q-\frac{n}{2})}{\Gamma(p)\Gamma(q)} \int_0^1 \frac{t^{p-1}(1-t)^{q-1}}{\sqrt{\det\big(t{\bf 1}+(1-t){\bf C}\big)}} {\rm d}t$$ Therefore if the eigenvalues of ${\bf C}$ are $\lambda_i$ we get $$ I = \frac{1}{\sqrt{\det{\bf A}}} \frac{\pi^{n/2}\Gamma(p+q-\frac{n}{2})}{\Gamma(p)\Gamma(q)} \int_0^1 t^{p-1}(1-t)^{q-1}\prod_{i=1}^n\big(t+(1-t)\lambda_i\big)^{-\frac12} {\rm d}t$$ The remaining integral looks like some kind of generalized hypergeometric function, but I don't know which one, and whether it has a simpler representation.