As noted, this integral is reducible to an elliptic integral and the form of the proposed solution may suggest that its modulus could be a singular value. Indeed, from the form
\begin{align}
I&=2\int_0^\infty \frac{\mathrm{d}u}{\sqrt {u^4+21u^2+112}} \\
&=2\int_{0}^\infty \frac{\mathrm{d}u}{\sqrt {u^2+\frac{21}{2}+\sqrt{\frac{7}{4}}i}{\sqrt {u^2+\frac{21}{2}-\sqrt{\frac{7}{4}}i}}} \\
&=2\int_{0}^\infty \frac{\mathrm{d}u}{\sqrt {\left( u^2+u_+^2\right)\left( u^2+u_-^2 \right)}}
\end{align}
The opposite of the roots of the polynomial under the square root ($u_{\pm}^2=(21\pm\sqrt{7})/2$) are complex conjugate, as the coefficients of the polynomial are real. In the following, we first apply a Landen's transformation where these roots are replaced by their geometric and arithmetic means which are real.
Defining $p=\sqrt{u_+u_-}$ and $q=(u_++u_-)/2$, we find from the coefficients of the polynomial
\begin{align}
p^2&=\sqrt{112}=4\sqrt{7}\\
q^2&=\frac{1}{4}\left( u_+^2+u_-^2+2u_+u_- \right)\\
&=\frac{1}{4}\left(21+2\sqrt{112} \right)=\frac{1}{4}\left(21+8\sqrt{7} \right)
\end{align}
With $u=x+\sqrt{x^2+p^2}$, remarking that
\begin{equation}
\sqrt {\left( u^2+u_+^2\right)\left(u^2+u_-^2 \right)}=2u\sqrt{x^2+q^2}
\end{equation}
and that
\begin{equation}
du=\frac{u}{\sqrt{x^2+p^2}}
\end{equation}
The integral becomes
\begin{equation}
I=\int_{-\infty}^\infty\frac{dx}{\sqrt{\left( x^2+p^2 \right)\left( x^2+q^2 \right)}}
\end{equation}
Now, changing $x=q\tan\theta$ and taking into account the parity of the function, one obtains an integral representation of the complete elliptic integral of the first kind:
\begin{align}
I&=2\int_0^{\pi/2} \frac{d\theta}{\sqrt{p^2\cos^2\theta+q^2\sin^2\theta}}\\
&=\frac{2}{p}\int_0^{\pi/2} \frac{d\theta}{\sqrt{1-\frac{p^2-q^2}{p^2}\sin^2\theta}}\\
&=\frac{2}{p}\mathbf{K}\left( \sqrt{1-\frac{q^2}{p^2}} \right)%\\
%&=\frac{2}{p}K'\left( \frac{q}{p} \right)
\end{align}
Now, $\sqrt{1-\frac{q^2}{p^2}} =\frac{\sqrt{2}}{8}\left( 3-\sqrt{7} \right)=\lambda^*(7)$, where $\lambda^*(r)$ gives a singular value of the elliptic modulus $k_r$ for which $\mathbf{K}\left( \sqrt{1-k_r^2} \right)=\sqrt{r}\mathbf{K}(k_r)$ see here, here and here. Then,
\begin{align}
I&=\frac{2}{p}\mathbf{K}\left( k_7 \right)\\
&=\frac{2}{2.7^{1/4}}\frac{\Gamma\left( \frac{1}{7} \right)\Gamma\left( \frac{2}{7} \right)\Gamma\left( \frac{4}{7} \right)}{4\pi 7^{1/4}}\\
&=\frac{\Gamma\left( \frac{1}{7} \right)\Gamma\left( \frac{2}{7} \right)\Gamma\left( \frac{4}{7} \right)}{4\pi \sqrt{7}}
\end{align}
as expected.
