I just found the following proof of Euler’s reflection formula. The thing I like about it is that allĀ  it requires is dumbly marching forwards, there is no trickery.

When s and 1-s are both positive, that is, when 0<s<1,

    \[\Gamma(s)\Gamma(1-s)=\int_0^\infty t^{s-1}e^{-t}dt\int_0^\infty t^{-s}e^{-t}dt.\]

I’m going to join those integrals, so giving different names to each copy of t,

    \[=\int_0^\infty\int_0^\infty x^{s-1}e^{-x}y^{-s}e^{-y}dxdy\]

    \[=\int_0^\infty\int_0^\infty \frac 1x \left(\frac xy\right)^s e^{-(x+y)}dxdy.\]

Now change variables u=\frac xy and v=x+y. The inverse transformation is x=\frac{uv}{u+1} and y=\frac v{u+1} and the Jacobian is \frac v{(u+1)^2}, then

    \[=\int_0^\infty\int_0^\infty\frac{u+1}{uv} u^s e^{-v}\frac v{(u+1)^2}dudv\]


To evaluate this integral, split into the two parts


    \[=\int_0^1\frac{u^{s-1}}{1+u}du+\int_1^\infty\frac{u^{s-2}}{1+\frac 1u}du.\]

Now use the series for the denominators

    \[=\int_0^1u^{s-1}\sum_{h=0}^\infty (-1)^hu^hdu+\int_1^\infty u^{s-2}\sum_{h=0}^\infty (-1)^hu^{-h}du\]

and now change the order of the sums and integrals

    \[=\sum_{h=0}^\infty\int_0^1(-1)^hu^{s+h-1}du+\sum_{h=0}^\infty\int_1^\infty (-1)^hu^{s-h-2}du\]

    \[=\sum_{h=0}^\infty\frac {(-1)^h}{s+h}+\frac {(-1)^h}{s-h-1}\]

    \[=\frac 1s+\sum_{h=1}^\infty \frac {(-1)^h 2h}{s^2-h^2}.\]

Now define the function

    \[f(x)=x\prod_{n=1}^\infty \frac{4n^2-x^2}{4n^2}\frac{(2n-1)^2}{(2n-1)^2-x^2}\]

and observe that

    \[f(x)=\left(x\prod_{n=1}^\infty 1-\frac{x^2}{4n^2}\right)\left(\prod_{n=0}^\infty 1-\frac{x^2}{(2n+1)^2}\right)^{-1}\]

and use Euler’s products for the sine and cosine

    \[\sin(\pi x)=\pi x\prod_{n=1}^\infty 1-\frac{x^2}{n^2}\quad \cos\left(\frac\pi 2 x\right)=\prod_{n=0}^\infty 1-\frac{x^2}{(2n+1)^2}\]

to notice

    \[f(x)=\frac 2\pi \sin\left(\frac\pi 2 x\right)\frac 1{\cos\left(\frac\pi 2 x\right)}=\frac 2\pi\tan\left(\frac\pi 2 x\right)\]

now take the logarithmic derivative

    \[\log(f(x))'=\frac{\pi}{\sin(\pi x)}=\frac 1x+\sum_{n=1}^\infty \frac {4n}{4n^2-x^2}-\frac{2(2n-1)}{(2n-1)^2-x^2}\]

the right hand side is what we had previously found with the variable s instead of x, therefore

    \[\Gamma(s)\Gamma(1-s)=\frac{\pi}{\sin(\pi s)}.\]

We have the equality when 0<s<1, and since these functions are holomorphic, and they coincide in an open set they must be equal.