Let $a_1,a_2>0$ be real numbers and $n_1,n_2\ge 1$ be integers. We consider a following family of definite integrals:
\begin{equation} {\mathfrak I}_{n_1,n_2}:= \int\limits_{-\infty}^\infty \frac{1}{(x^2+a_1^2)^{n_1}(x^4+a_2^4)^{n_2}} dx \end{equation}
By decomposing the integrand into partial fractions, meaning by using A partial fraction decomposition of an inverse of a generic polynomial of an arbitrary order. we have found the result for $n_1=1$. It reads:
\begin{equation} {\mathfrak I}_{1,n_2}= \pi \left(\frac{\sum\limits_{j=0}^{4n_2-1} {\mathfrak P}_{n_2,j} a_1^{4 n_2-j} a_2^j}{a_1 a_2^{4 n_2-1} \left(a_1^4+2 \sqrt{2} a_1^3 a_2+4 a_1^2 a_2^2+2 \sqrt{2} a_1 a_2^3+a_2^4\right)^{n_2}}\right) \end{equation}
Here the matrix ${\mathfrak P} := \left({\mathfrak P}_{n_2,j}\right)_{n_2=1,j=0}^{3,4n_2-1} $ reads:
\begin{equation} {\mathfrak P}:= \left( \begin{array}{c} \left\{\frac{1}{\sqrt{2}},2,\frac{3}{\sqrt{2}},1\right\} \\ \left\{\frac{3}{4 \sqrt{2}},3,\frac{47}{4 \sqrt{2}},14,\frac{87}{4 \sqrt{2}},11,\frac{27}{4 \sqrt{2}},1\right\} \\ \left\{\frac{21}{32 \sqrt{2}},\frac{63}{16},\frac{751}{32 \sqrt{2}},45,\frac{1959}{16 \sqrt{2}},\frac{247}{2},\frac{3009}{16 \sqrt{2}},108,\frac{2937}{32 \sqrt{2}},\frac{441}{16},\frac{339}{32 \sqrt{2}},1\right\} \\ \end{array} \right) \end{equation}
In[544]:= {a1, a2} = RandomReal[{0, 5}, 2, WorkingPrecision -> 50];
n1 = 1;
l1 = Table[
NIntegrate[
1/(x^2 + a1^2)^n1 1/(x^4 + a2^4)^n2, {x, -Infinity, Infinity},
WorkingPrecision -> 30], {n2, 1, 3}];
PP = {{1/Sqrt[2], 2, 3/Sqrt[2], 1}, {3/(4 Sqrt[2]), 3, 47/(4 Sqrt[2]),
14, 87/(4 Sqrt[2]), 11, 27/(4 Sqrt[2]), 1}, {21/(32 Sqrt[2]), 63/
16, 751/(32 Sqrt[2]), 45, 1959/(16 Sqrt[2]), 247/2, 3009/(
16 Sqrt[2]), 108, 2937/(32 Sqrt[2]), 441/16, 339/(32 Sqrt[2]), 1}};
l2 = Table[
Pi Sum[PP[[n2, 1 + j]] a1^(4 n2 - 1 - j) a2^j, {j, 0, 4 n2 - 1}]/(
a1 a2^(-1 + 4 n2) (a1^4 + 2 Sqrt[2] a1^3 a2 + 4 a1^2 a2^2 +
2 Sqrt[2] a1 a2^3 + a2^4)^n2), {n2, 1, 3}];
l1 - l2
Out[549]= {0.*10^-32, 0.*10^-32, 0.*10^-33}
Apart from the obvious question as to how does the result looks like for generic values of $n_1 \ge 1 $ another question would be how does this method of computing integrals compare the one based on calculating residues? More generally, which method, the partial fraction decomposition or the residue theorem is better suited -- more efficient-- in computing such integrals?
Let us take $d=6$ and define $\left( a_\xi\right)_{\xi=1}^d = \left( \imath a_1, -\imath a_1, e^{\imath \pi/4}a_2, e^{\imath 3 \pi/4}a_2, e^{\imath 5\pi/4}a_2, e^{\imath 7\pi/4}a_2\right) $ then $\left(n_\xi \right)_{\xi=1}^d = \left(n_1,n_1,n_2,n_2,n_2,n_2 \right) $. Then define $\chi(i):= 1_{i=1}-1_{i=2} +1_{i=3}+1_{i=4}-1_{i=5}-1_{i=6}$ Then by applying the partial fraction decomposition formula from A partial fraction decomposition of an inverse of a generic polynomial of an arbitrary order. and by the virtue of the fact that $\int\limits_{-\infty}^\infty 1/(z+a) dz = -\imath \pi \cdot \mbox{sign[Im[a]]} $ we have:
\begin{eqnarray} &&{\mathcal I}_{n_1,n_2} = \\ &&\sum\limits_{i=1}^d \sum\limits_{0 \le j_1 \le j_2 \le \cdots \le j_{d-1} \le n_i-1} \left[ \prod\limits_{\begin{array}{lll} \xi=1,\xi\neq i \\ \lambda:= \xi - 1_{\xi \ge i+1}\end{array}}^d \binom{n_\xi+j_\lambda-j_{\lambda-1}-1}{n_\xi-1} \frac{1}{\left(a_\xi - a_i \right)^{n_\xi+j_\lambda-j_{\lambda-1}}} \right] \cdot (-1)^{j_{d-1}} 1_{n_i-j_{d-1}==1} \cdot \chi(i) \cdot (-\imath \pi) \end{eqnarray}