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?

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$.

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.$$

Let,

$$I=\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$$

As noted by Jon Haussmann,

$$2I=\int_0^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx + \int_0^{\infty}\frac{(x+1)^2}{x^6 +4x^5 + 3x^4 - 4x^3 - 2x^2 + 2x + 1}dx$$

Perform the change of variable $x=\dfrac{1}{u-1}$ in the second integral,

$\begin{align}2I&=\int_0^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\int_1^{\infty} \frac{x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\\ &=\left(\int_0^1 \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\\ \int_1^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\right)+\\&\int_1^{\infty} \frac{x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\end{align}$

Perform the change of variable $x=\dfrac{u-1}{u}$ in the first integral of the latter equality,

$\begin{align}2I&=\int_1^{\infty} \frac{x^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\int_1^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\int_1^{\infty} \frac{x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\\ &=\int_1^{\infty} \frac{x^2+(x-1)^2+x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\\ &=\Big[\arctan\left(\dfrac{x^3-2x^2-x+1}{x(x-1)}\right)\Big]_0^{+\infty}\\ &=\pi \end{align}$

(Problem found in American Mathematical Monthly, vol. 112, april 2005.

Solution found in vol. 114)