I encountered an astonishing integral (numerically verified):
$$\int_{-\infty}^\infty \frac{x\ \mathrm dx}{\sin ^2(\sqrt{x}) \sinh ^2(2 \sqrt{2 x})+\pi ^2 \cos ^2(\sqrt{x}) \cosh ^2(2 \sqrt{2 x})}\\ \small =-\frac{1}{262144 \pi ^3}\left(18 \sqrt{2}+4 \sqrt{34 \sqrt{2}+46}+13\right) \Gamma \left(\frac{1}{8}\right)^4 \Gamma \left(\frac{3}{8}\right)^4$$
What technique should be used to establish it? The gamma factors suggest that Elliptic Singular Values may be involved. Any help will be appreciated.
Let $K,K'$ be positive real numbers, $\mathcal{H} = \{z\in \mathbb{C} | \Im z > 0 \}$, $$f(z) = \cos (\sqrt{z} \frac{K-iK'}{2}) = \underbrace{\cos(\sqrt{z}\frac{K}{2})\cosh(\sqrt{z}\frac{K'}{2})}_{b(z)}+i\times\underbrace{\sin(\sqrt{z}\frac K2)\sinh(\sqrt{z}\frac{K'}{2})}_{d(z)}$$ $f,b,d$ are entire functions and all zeroes of $f$ are in $\mathcal{H}$. Write $\overline{f}(z) = \overline{f(\overline{z})}$, one sees $|f(z)/\overline{f}(z)| \leq 1$ for $z\in \mathcal{H}$.
For $a>0$, let $$\DeclareMathOperator{\sn}{sn} \DeclareMathOperator{\cd}{cd} \mu(x;a) = \frac{a}{d(x)^2 + a^2 b(x)^2} \qquad x\in \mathbb{R}$$
Fix a non-negative integer $n$, a non-trivial observation is that
Proof: Let $w = (1+ai)/(1-ai)$, we have $|w|<1$. One easily checks (which is true for generic $f$ and $w$) $$\mu(x;a) = \mu(x;1) + \frac{w}{\overline{f} (\overline{f}-wf)} + \frac{\overline{w}}{f (f-\overline{w}\overline{f})}$$ so for $R$ large $$\int_{|x|<R} x^n \mu(x;a) dx = \int_{|x|<R} x^n \mu(x;1) dx + 2 \Re \underbrace{\int_{|x|<R} \frac{x^n w}{\overline{f}(x) (\overline{f}(x) - w f(x))} dx}_{I(R)}$$ We argue $\lim_{R\to\infty} I(R) = 0$. Note that $\overline{f}$ has no zero in $\mathcal{H}$ and also $|1-w f/\overline{f}| > 1-|w|$ in $\mathcal{H}$. Thus we can deform the contour to semicircle of radius $R$ in $\mathcal{H}$, $$|I(R)| \leq R^{n+1} \frac{1}{1-|w|} \int_0^\pi \frac{d\theta}{|\overline{f}(Re^{i\theta})|^2}$$ now $\overline{f}$ increases fast enough on $\mathcal{H}$, such that RHS tends to zero as $R\to \infty$, as desired. QED.
So we denote $$I_n = \int_\mathbb{R} x^n \mu(x;1) dx = \int_\mathbb{R} \frac{1}{2} \frac{x^n dx}{\cos(\sqrt{x}K)+\cosh(\sqrt{x}K')}$$ We need to evaluate the last integral.
So far we treat $K,K'$ as generic positive reals, now we demand $K=K(m), K'=K(1-m)$ for some $0<m<1$ with $K(x) = \frac{\pi}{2} {_2F_1}(\frac{1}{2},\frac{1}{2};1;x)$ the complete elliptic integral. Under this assumption:
Proof sketch: Write LHS as $J$, making $x\mapsto x^2$, the path of integration becomes a L-shape from $i\infty\to 0 \to \infty$, summing residue in first quadrant: $$J=\frac{2\pi i}{K-iK'} \sum_{n\geq 0} (-1)^{n+1} \frac{\sin(\pi(n+1/2) \frac{2u}{K-iK'})}{\cos(\pi(n+1/2) \frac{K+iK'}{K-iK'})}$$ This is the Fourier series of Jacobi's $\text{sd}$, so $$J = \frac{2i}{K-iK'} (-Kkk') \times \text{sd}(K\frac{2u}{K-iK'} | \tau_1)$$ where $\tau_1 = (\tau+1)/(-\tau+1)=\left(\begin{smallmatrix}1& 1 \\ -1& 1\end{smallmatrix}\right)\tau$ with $\tau$ corresponds to our original $K,K'$. One wants to transform above $\text{sd}(\cdot | \tau_1)$ to Jacobi elliptic function with respect to $\tau$ instead of $\tau_1$.
To do this, expand all terms of $J$ into Jacobi's theta function, from which $\tau \mapsto 1+\tau, \tau\mapsto -1/\tau$ (which generates $\text{SL}_2(\mathbb{Z})$) are easily done. The above matrix has determinant $2$, so a $\tau \mapsto 2\tau$ has to be performed in between, this is done via the classical Landen transformation. After all these steps, one puts the result back to Jacobi's elliptic function, to get RHS. See here for more details. QED.
By expanding around $u=0$: for $n=0,1,2,3$ respectively
$$\tag{2}\frac{(-1)^n}{(2n+1)!} I_n = \left\{1,\frac{1-2m}{3},\frac{2(m^2-m+1)}{15}, \frac{(1-2m)(2 m^2-2 m+17)}{315},\cdots\right\} $$
We have now proved, for any $a>0$ and $I_n$ in $(2)$: $$\int_{-\infty}^\infty \frac{ a x^n \mathrm dx}{\sin^2(\sqrt{x})\sinh^2(\frac{K'}{K}\sqrt{x})+a ^2 \cos ^2(\sqrt{x})\cosh^2(\frac{K'}{K}\sqrt{x})} = (\frac{K}{2})^{2n+2} I_n$$
All remain is a straightforward calculation of elliptic singular value. Here we want to find $m$ and $K$ such that $K(1-m)/K(m) = 2\sqrt{2}$. For small cases like this, their values are well-tabulated:
looking up here $$m = (\sqrt{2}+1-\sqrt{2\sqrt{2}+2})^2$$
looking up here $$K = \frac{1}{16} \sqrt[4]{\frac{1}{2} \left(\sqrt{2}+1\right)} \sqrt{\frac{2 \sqrt{2}+\sqrt{5 \sqrt{2}+1}}{\pi }} \Gamma \left(\frac{1}{8}\right) \Gamma \left(\frac{3}{8}\right)$$
OP's integral is the case $n=1$, a comment mentions the case $n=3$.
Another comment gives an article for $m=1/2 \iff K=K'$, which is a special case of our general considerations.