How did wolframalpha evaluate integrate $\int_0^\infty \sin(x)/(x^2+1) dx =(\operatorname{Ei}(1)-e^2\operatorname{Ei}(-1))/(2e)$? I have been working on th integral
$$\int_0^\infty \frac{\sin x}{1+x^2} dx$$
for a short while now trying substitutions and even Laplace transforms and other stuff but gave up. I looked to Wolfram Alpha but how did it get the series expansion of the constant that it did (see here)? I tried doing a little algebra with the expression from Wolfram Alpha  it looks like an integral and there are more things that could be done but I don't know. Thanks.
 A: $$\newcommand{\Re}{\operatorname{Re}}\newcommand{\Im}{\operatorname{Im}}\newcommand{\Ei}{\operatorname{Ei}}\newcommand{\PV}{\operatorname{PV}}
\begin{align}
\int_0^\infty\frac{\sin(x)}{1+x^2}\,\mathrm{d}x
&=\Im\left(\int_0^\infty\frac{e^{ix}}{1+x^2}\,\mathrm{d}x\right)\tag1\\
&=\frac12\Im\left(\int_0^\infty e^{ix}\left(\frac1{1-ix}+\frac1{1+ix}\right)\mathrm{d}x\right)\tag2\\
&=\frac12\Im\left(\int_0^{i\infty}e^{ix}\left(\frac1{1-ix}+\frac1{1+ix}\right)\mathrm{d}x\right)\tag3\\
&=\frac12\Re\left(\int_0^\infty e^{-x}\left(\frac1{x+1}-\frac1{x-1}\right)\mathrm{d}x\right)\tag4\\
&=\frac e2\int_1^\infty\frac{e^{-x}}{x}\,\mathrm{d}x-\frac1{2e}\PV\!\!\int_{-1}^\infty\frac{e^{-x}}{x}\,\mathrm{d}x\tag5\\
&=\frac1{2e}\left(\Ei(1)-e^2\Ei(-1)\right)\tag6
\end{align}
$$
Explanation:
$(1)$: $\sin(x)=\Im\left(e^{ix}\right)$
$(2)$: partial fractions
$(3)$: move contour to the positive imaginary axis
$\phantom{\text{(3):}}$ with a small counterclockwise bump around $i$
$(4)$: substitute $x\mapsto ix$ using $\Im(iz)=\Re(z)$
$(5)$: the small counterclockwise bump around $1$
$\phantom{\text{(5):}}$ gets turned into $\PV$ by $\Re$
$(6)$: apply $\Ei(z)=-\PV\!\int_{-z}^\infty\frac{e^{-x}}{x}\,\mathrm{d}x$
A: In the hope that it becomes a bit more visible, what happens behind the scenes, I tried to extract all steps from the linked notebook:
$\frac{\text{Ei}(1)-e^2 \text{Ei}(-1)}{2 e}=\sum _{k=1}^{\infty } -\frac{e^2 (-1)^k-1}{2 e k k!}-\frac{e \gamma }{2}+\frac{\gamma }{2 e}$
$\frac{\text{Ei}(1)-e^2 \text{Ei}(-1)}{2 e}=-\frac{e \sum _{k=1}^{\infty } \sum _{j=1}^k \frac{(-1)^k}{j k!}-e \sum _{k=1}^{\infty } \sum _{j=1}^k \frac{(-1)^{2 k}}{j k!}+e^2 \gamma -\gamma }{2 e}$
$\frac{\text{Ei}(1)-e^2 \text{Ei}(-1)}{2 e}=\sum _{k=1}^{\infty } \frac{(-1)^k \left(e^2 \left(-z_0-1\right){}^k-\left(1-z_0\right){}^k\right) \Gamma \left(k,-z_0\right)}{2 e k! z_0^k}+\frac{\text{Ei}\left(z_0\right)}{2 e}-\frac{e \text{Ei}\left(z_0\right)}{2}$
$\frac{\text{Ei}(1)-e^2 \text{Ei}(-1)}{2 e}=\frac{e^{x+2} \sum _{k=1}^{\infty } \sum _{j=0}^{k-1} \frac{(-1)^{k-j} (-x-1)^k x^{j-k}}{k j!}-e^x \sum _{k=1}^{\infty } \sum _{j=0}^{k-1} \frac{(-1)^{k-j} (1-x)^k x^{j-k}}{k j!}-(2 i) e^2 \pi  \left\lfloor \frac{\arg (-x-1)}{2 \pi }\right\rfloor +(2 i) \pi  \left\lfloor \frac{\arg (1-x)}{2 \pi }\right\rfloor +\text{Ei}(x)-e^2 \text{Ei}(x)+i \pi }{2 e}$
$\frac{\text{Ei}(1)-e^2 \text{Ei}(-1)}{2 e}=\sum _{k=1}^{\infty } -\frac{\left((1-x)^k-e^2 (-x-1)^k\right) \left((-1)^k k!-k x \, _2\tilde{F}_2(1,1;2,2-k;x)\right)}{2 e k k! x^k}-i e \pi  \left\lfloor \frac{\arg (-x-1)}{2 \pi }\right\rfloor +\frac{i \pi  \left\lfloor \frac{\arg (1-x)}{2 \pi }\right\rfloor }{e}+\frac{\text{Ei}(x)}{2 e}-\frac{e \text{Ei}(x)}{2}+\frac{i \pi }{2 e}$
Note that $\Gamma(x)$ is the Gamma function, $\text{Ei}(x)$ denotes the exponential integral, $\gamma$ is the Euler-Mascheroni Constant, $\arg(z)$ is the complex argument and $\tilde{F}$ denotes the Generalized Hypergeometric Function.
