How does one evaluate $I(\alpha,\beta,\gamma)= \int_{0}^{\infty} \cos\left(\frac{x(x^2-\alpha^2)}{x^2-\beta^2}\right)\frac{1}{x^2+\gamma^2} \, dx.$ I have been trying to evaluate this integral for a while now,
$$I(\alpha,\beta,\gamma)= \int_{0}^{\infty} \cos\left(\frac{x(x^2-\alpha^2)}{x^2-\beta^2}\right)\frac{1}{x^2+\gamma^2} \, dx.$$
I first came across this integral on Quora (https://www.quora.com/What-improper-integrals-are-hard-to-solve).
Someone stated that this integral was submitted as a problem to the Gazette of the Royal Mathematics Society of Spain and was (at the time) still open, so the complete solution would not be published there. (have not been able to find the publication). 
Although no derivation, a closed form expression was given:
$$I(\alpha,\beta,\gamma)=\frac{\pi}{2\gamma}e^{-\frac{\gamma(\alpha^2+\gamma^2)}{\beta^2+\gamma^2}}.$$
One might be wondering, why not just respond to the post on Quora? That I did. However, the thread seems to inactive.
I have tried several methods including the one employed by Mark Viola in this post (Computing $\int_{-\infty}^{\infty} \frac{\cos x}{x^{2} + a^{2}}dx$ using residue calculus). Trying to work backwards from the closed form expression, by treating the term $ \frac{\gamma(\alpha^2+\gamma^2)}{\beta^2+\gamma^2}$ in exponent as one variable. However, I have not been able to succeed. 
Hopefully, you can help me solve this integral/ give suggestions/ different methods.
edit: I found that when $\alpha=\beta$, the integral becomes
$I= \int_{0}^{\infty} \cos(x)\frac{1}{x^2+\gamma^2}dx$ , which evaluates to $I=\frac{\pi}{2\gamma}e^{-\gamma}$. This is comparable to the desired result. The argument in the exponent is different however.   
 A: For $\alpha, \beta \geq 0$ and $\gamma > 0$ let
$$ f_{\alpha,\beta,\gamma} (z) = \frac{\exp \left(\mathrm{i} z \frac{z^2 - \alpha^2}{z^2 - \beta^2}\right)}{z^2 + \gamma^2} \, , \, z \in \mathbb{C} \setminus \{\pm\beta, \pm \mathrm{i} \gamma\} \, . $$
By symmetry we have
$$ I (\alpha,\beta,\gamma) = \frac{1}{2} \int \limits_{-\infty}^\infty f_{\alpha,\beta,\gamma} (x) \, \mathrm{d} x \, .$$
Clearly, the integral exists for every combination of parameters and is bounded by $\frac{\pi}{2 \gamma}$ .

We can compute its value using the residue theorem. First note that we have
$$ \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) = \frac{1}{2  \mathrm{i} \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) $$
and (by Jordan's lemma)
$$ \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z = 0 \, ,$$
where $\Gamma_R$ is a semi-circle of radius $R > 0$ in the upper half-plane. A naive application of the residue theorem therefore yields
$$ I (\alpha,\beta,\gamma) = \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z \right]= \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) \, .$$
While the result is indeed correct, this calculation is only valid for $\alpha = \beta$ . In this case we obtain (as already mentioned in the question) $I(\alpha,\alpha,\gamma) = \frac{\pi}{2 \gamma}\mathrm{e}^{- \gamma}$ .
If $\alpha \neq \beta$ , $f_{\alpha,\beta,\gamma}$ has essential singularities on the real axis, namely at $\pm \beta$ . Thus we need to deform our contour using small semi-circles $\gamma_\varepsilon (\pm \beta)$ of radius $\varepsilon > 0$ centred at these points and show that their contribution to the integral vanishes as $\varepsilon \to 0$ .

For simplicity, I will only discuss the case $\beta = 0$ in detail. We want to find
$$ \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z = \varepsilon \int \limits_0^\pi \frac{\mathrm{e}^{\mathrm{i} \left[\phi + \varepsilon \mathrm{e}^{\mathrm{i}\phi} - \alpha^2 \varepsilon^{-1} \mathrm{e}^{-\mathrm{i} \phi}\right]}}{\gamma^2+\varepsilon^2 \mathrm{e}^{2 \mathrm{i} \phi}} \, \mathrm{d} \phi = \frac{\varepsilon}{\gamma^2} \int \limits_0^\pi \mathrm{e}^{\mathrm{i} \left[\phi - \alpha^2 \varepsilon^{-1} \mathrm{e}^{-\mathrm{i} \phi}\right]} \left[1 + \mathcal{O} (\varepsilon) \right] \, \mathrm{d} \phi \, .$$
The leading-order term can actually be calculated analytically: for $z \in \mathbb{C}$ we have
$$ \int \limits_0^\pi \mathrm{e}^{\mathrm{i} \left[\phi - z \mathrm{e}^{-\mathrm{i} \phi}\right]} \, \mathrm{d} \phi = \mathrm{i} \left[(2 \operatorname{Si}(z) - \pi) z + 2 \cos(z)\right] \, . $$
The sine integral $\operatorname{Si}$ satisfies $\lim_{x \to \infty} \operatorname{Si}(x) = \frac{\pi}{2}$, so we obtain
$$ \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z \sim \frac{\mathrm{i} \varepsilon}{\gamma^2} \left[(2 \operatorname{Si}(\alpha^2 \varepsilon^{-1}) - \pi) \alpha^2 \varepsilon^{-1} + 2 \cos(\alpha^2 \varepsilon^{-1})\right] \stackrel{\varepsilon \to 0}{\longrightarrow} 0 $$
as desired. Now we are allowed to apply the residue theorem , which yields
\begin{align}
I(\alpha,0,\gamma) &= \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,0,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (0)} f_{\alpha,0,\gamma} (z) \, \mathrm{d} z \right] \\
&= \frac{\pi}{2 \gamma} \exp \left(- \frac{\alpha^2 + \gamma^2}{\gamma}\right) \, .
\end{align}

Almost the same calculation yields 
\begin{align}
I(\alpha,\beta,\gamma) &= \frac{1}{2} \left[ 2 \pi \mathrm{i} \operatorname{Res} (f_{\alpha,\beta,\gamma}, \mathrm{i} \gamma) - \lim_{R \to \infty} \int \limits_{\Gamma_R} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (\beta)} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z - \lim_{\varepsilon \to 0} \int \limits_{\gamma_\varepsilon (-\beta)} f_{\alpha,\beta,\gamma} (z) \, \mathrm{d} z \right] \\ 
&= \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right)
\end{align}
for $\beta > 0$ . 

Combining these results we conclude that the integral is given by
$$ I(\alpha,\beta,\gamma) = \frac{\pi}{2 \gamma} \exp \left(- \gamma \frac{\alpha^2 + \gamma^2}{\beta^2 + \gamma^2}\right) $$
for every $\alpha,\beta \geq 0$ and $\gamma > 0$ .
