A nasty integral of a rational function I'm having a hard time proving the following $$\int_0^{\infty} \frac{x^8 - 4x^6 + 9x^4 - 5x^2 + 1}{x^{12} - 10 x^{10} + 37x^8 - 42x^6 + 26x^4 - 8x^2 + 1} \, dx = \frac{\pi}{2}.$$
Mathematica has no problem evaluating it while I haven't the slightest idea how to approach it. Of course, I would like to prove it without the use of a computer. Is this an explicit form of a special function I fail to recognize?
 A: The integral over $\mathbb{R}$ of a meromorphic function $f(z)$, $O(|z|^{-2})$ at infinity, non-vanishing over $\mathbb{R}$, is equal to $2\pi i$ times the sum of residues in the poles located in the complex upper-half plane. Since:
$$p(y) = y^6-2y^5-2y^4+4y^3+3y^2-4y+1 = p_{+}(y)\cdot p_{-}(y),$$
$$p_{+}(y)= y^3-(i+1)y^2+(i-2)x+1,\qquad p_{-}(y)=y^3+(i-1)y^2-(2+i)y+1,$$
(I got this through a numerical calculation of the roots of $p(y)$, followed by a separation of the roots with positive and negative imaginary part, say $\zeta_1,\zeta_2,\zeta_3$ and $\bar{\zeta_1},\bar{\zeta_2},\bar{\zeta_3}$ - so $p_{+}(z)$ is just $\prod_{j=1}^3 (z-\zeta_j)$) we have:
$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = 2\pi i\sum_{j=1}^{3}\operatorname{Res}_{z=\zeta_j}\left(\frac{z^2}{p_{+}(z)\cdot p_{-}(z)}\right),$$
but $p_{-}(x)-p_{+}(x)=2i(x^2-x)$, so:
$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = 2\pi i\sum_{j=1}^{3}\operatorname{Res}_{z=\zeta_j}\left(\frac{z^2}{p_{+}^2(z)+2i(z^2-z)p_{+}(z)}\right).$$
By De l'Hopital theorem, and since $\zeta_j$ is a double zero of $p_{+}^2(x)$:
$$\lim_{z\to\zeta_j}\frac{z^2(z-\zeta_j)}{p_{+}^2(z)+2i(z^2-z)p_{+}(z)}=\frac{\zeta_j^2}{2i(\zeta_j^2-\zeta_j)p_{+}'(\zeta_j)},$$
so:
$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = \pi\sum_{j=1}^{3}\frac{\zeta_j}{(\zeta_j-1)p_{+}'(\zeta_j)}.$$
Now we compute the remainder between $(z-1)p_{+}'(z)$ and $p_{+}(z)$, in order to have:
$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = \pi\sum_{j=1}^{3}\frac{\zeta_j}{-(1+i)+6\zeta_j-(2-i)\zeta_j^2}.$$
If now we take $\alpha=\frac{3+\sqrt{6-i}}{2-i}$ and $\beta=\frac{3-\sqrt{6-i}}{2-i}$ we can re-write the last line as:
$$ I = \frac{\pi}{(i-2)(\alpha-\beta)}\left(\sum_{j=1}^{3}\frac{\alpha}{\zeta_j-\alpha}-\sum_{j=1}^{3}\frac{\beta}{\zeta_j-\beta}\right)=-\frac{\pi}{2\sqrt{6-i}}\left(\Sigma_1-\Sigma_2\right).$$
Now $\Sigma_1$ is the sum of the reciprocal of the roots of the polynomial $p_{+}(\alpha(z+1))$, and $\Sigma_2$ is the sum of the reciprocal of the roots of the polynomial $p_{+}(\beta(z+1))$. This quantities can be computed through the coefficients of $p_{+}$, since the sum of the reciprocal of the roots of a polynomial $q(z)$ is just $-\frac{q'(0)}{q(0)}$. This gives:
$$\Sigma_1 = -\alpha\frac{p_{+}'(\alpha)}{p_{+}(\alpha)},\qquad \Sigma_2 = -\beta\frac{p_{+}'(\beta)}{p_{+}(\beta)}.$$
Up to a massive amount of long but straightforward computations, we get:
$$\Sigma_1 = (i-2)-\sqrt{6-i},\qquad \Sigma_2 = (i-2)+\sqrt{6-i}, $$
from which $\color{red}{I=\pi}$ finally follows.
I am really grateful to Jon Haussmann for the proof that
$$\int_0^{\infty} \frac{x^8 - 4x^6 + 9x^4 - 5x^2 + 1}{x^{12} - 10 x^{10} + 37x^8 - 42x^6 + 26x^4 - 8x^2 + 1}dx = \frac{1}{2}\int_{\mathbb{R}}\frac{y^2 dy}{p(y)},$$
where only the second integral is treated here.
IMPORTANT UPDATE:
In fact, there is no need to compute the coefficients of $p_{+}(x)$ and $p_{-}(x)$ (we only need the identity $p_{-}(x)-p_{+}(x)=2i(x^2-x)$), or introduce $\alpha$ and $\beta$. Since $p_{+}(x)$ is a third-degree polynomial with roots in the upper half-plane,
$$0=\int_{\mathbb{R}}\frac{dz}{p_{+}(z)}=\sum_{j=1}^{3}\frac{1}{p_{+}'(\zeta_j)}.$$
This gives:
$$ I = \pi\sum_{j=1}^{3}\frac{\zeta_j}{(\zeta_j-1)p_{+}'(z)} = \pi\sum_{j=1}^{3}\frac{1}{(\zeta_j-1)p_{+}'(\zeta_j)}, $$
but if we decompose $\frac{1}{p_{+}(z)}$ in simple fractions, we get:
$$\frac{1}{p_{+}(z)}=\sum_{j=1}^{3}\frac{1}{p_{+}'(\zeta_j)(z-\zeta_j)}, $$
so the magic gives:
$$ I = -\frac{\pi}{p_{+}(1)}.$$
Since $p(x)=p_{+}(x)\cdot p_{-}(x)$, $p(1)=1$, $I\in\mathbb{R}^+$ and $p_{-}(1)$ is the conjugate of $p_{+}(1)$, $p_{+}(1)$ can be only $+1$ or $-1$, so $I=\pi$.
A: Some progess: The integrand actually decomposes as
$$\frac{1}{2} \left( \frac{x^2 + 2x + 1}{x^6 + 4x^5 + 3x^4 - 4x^3 - 2x^2 + 2x + 1} + \frac{x^2 - 2x + 1}{x^6 - 4x^5 + 3x^4 + 4x^3 - 2x^2 - 2x + 1} \right).$$
Note that the second term is the same as the first term, except with $-x$ instead of $x$.  Thus, with some substitutions, the integral becomes
$$\frac{1}{2} \int_{-\infty}^\infty \frac{y^2}{y^6 - 2y^5 - 2y^4 + 4y^3 + 3y^2 - 4y + 1} \; dy.$$
