[tex]\Gamma(z):=\displaystyle\int_0^\infty s^{z-1}e^{-s}\,\mathrm ds[/tex]
[tex]\displaystyle\Gamma(z)\Gamma(1-z)=\left(\int_0^\infty s^{z-1}e^{-s}\,\mathrm ds\right)\left(\int_0^\infty t^{(1-z)-1}e^{-t}\,\mathrm dt\right)[/tex]
[tex]=\displaystyle\int_0^\infty\int_0^\infty\frac1s\left(\frac st\right)^ze^{-(s+t)}\,\mathrm ds\,\mathrm dt[/tex]
Consider the change of coordinates
[tex]\begin{cases}\mathbf s(u,v)=\dfrac v{u+1}\\\\\mathbf t(u,v)=\dfrac{uv}{u+1}\end{cases}[/tex]
which yields a Jabobian determinant of
[tex]\det\mathbf J=\begin{vmatrix}\mathbf s_u&\mathbf s_v\\\mathbf t_u&\mathbf t_v\end{vmatrix}=-\dfrac v{(u+1)^2}[/tex]
Then the double integral is equivalent to
[tex]\displaystyle\int_0^\infty\int_0^\infty\frac{u^{z-1}}{u+1}e^{-v}\,\mathrm du\,\mathrm dv=\int_0^\infty\frac{u^{z-1}}{u+1}\,\mathrm du[/tex]
As a function over the positive reals, this will only converge if [tex]0<z<1[/tex]. We can find the value of the integral by considering the complex-valued function,
[tex]f(\zeta)=\dfrac{\zeta^{z-1}}{\zeta+1}[/tex]
and integrate it over the contour [tex]\mathcal C[/tex] consisting of a circle [tex]\mathcal C_R[/tex] of radius [tex]R>0[/tex] connected to a smaller circle [tex]\mathcal C_\varepsilon[/tex] of radius [tex]\varepsilon[/tex] (both centered at the origin) by two line segments parallel and close to the positive real axis (but not touching it), oriented in the opposite direction relative to one another and denoted [tex]\gamma_{\varepsilon\to R}[/tex] and [tex]\gamma_{R\to\varepsilon}[/tex], respectively. (See attachment)
By the residue theorem, the value of the contour integral will be the sum of the residues at the poles of [tex]f(\zeta)[/tex] multiplied by [tex]2\pi i[/tex]. We have only one simple pole at [tex]\zeta=-1[/tex], which has residue
[tex]\mathrm{Res}\left(\dfrac{\zeta^{z-1}}{\zeta+1},-1\right)=\displaystyle\lim_{\zeta\to-1}(\zeta+1)\frac{\zeta^{z-1}}{\zeta+1}=\lim_{\zeta\to-1}\zeta^{z-1}=-e^{i\pi z}[/tex]
So we have
[tex]-e^{i\pi z}=\displaystyle\left\{\int_{\mathcal C_R}+\int_{\gamma_{R\to\varepsilon}}+\int_{\mathcal C_\varepsilon}+\int_{\gamma_{\varepsilon\to R}}}\right\}f(\zeta)\,\mathrm d\zeta[/tex]
By the ML lemma and the restriction of [tex]0<z<1[/tex], we have as [tex]R\to\infty[/tex] and [tex]\varepsilon\to0[/tex] that
[tex]\displaystyle\left|\int_{\mathcal C_R}f(\zeta)\,\mathrm d\zeta\right|=\left|\int_0^{2\pi}\frac{(Re^{i\theta})^{z-1}}{Re^{i\theta}+1}iRe^{i\theta}\,\mathrm d\theta\right|\le\dfrac{2\pi R^z}{|Re^{i\theta}+1|}\to0[/tex]
[tex]\displaystyle\left|\int_{\mathcal C_\varepsilon}f(\zeta)\,\mathrm d\zeta\right|=\left|\int_0^{2\pi}\frac{(\varepsilon e^{i\theta})^{z-1}}{\varepsilon e^{i\theta}+1}i\varepsilon e^{i\theta}\,\mathrm d\theta\right|\le\dfrac{2\pi \varepsilon^z}{|\varepsilon e^{i\theta}+1|}\to0[/tex]
We're left with
[tex]\displaystyle-2\pi ie^{i\pi z}=\left\{\int_{\gamma_{\varepsilon\to R}}+\int_{\gamma_{R\to\varepsilon}}\right\}f(\zeta)\,\mathrm d\zeta[/tex]
Note that the integral along [tex]\gamma_{\varepsilon\to R}[/tex] corresponds to the integral we wanted to compute in the first place, so we can replace [tex]\zeta=u[/tex]. For the other, we write the numerator of [tex]f(\zeta)[/tex] as [tex]\zeta^{z-1}=r^{z-1}e^{i(z-1)\theta}e^{2i\pi(z-1)}[/tex], to account for the fact that we're considering a particular branch of [tex]f(\zeta)[/tex]. As [tex]R\to\infty[/tex] and [tex]\varepsilon\to0[/tex], we're left with
[tex]-2i\pi e^{i\pi z}=\left(1-e^{2i\pi(z-1)}\right)\displaystyle\int_0^\infty\frac{u^{z-1}}{u+1}\,\mathrm du[/tex]
[tex]\implies\displaystyle\int_0^\infty\frac{u^{z-1}}{u+1}\,\mathrm du=\frac\pi{\sin\pi z}=\Gamma(z)\Gamma(1-z)[/tex]
as required.
