A popular method of proving the formula is to use the infinite product representation of the gamma function. See ProofWiki for example.
However, I'm interested in down-to-earth proof; e.g. using the change of variables. As the formula being connected to the beta function, there could be one-line proof for it.
Could anyone help me?
I suspect the following is a three line proof to the right reader: Set $e^g = \Gamma(1+z) \Gamma(1-z) z \sin(\pi z)$. Then $\Re(g)$ is an even harmonic function with $g(z) = O(|z| \log |z|)$, so $g$ is constant. Plugging in $z=1/2$ evaluates the constant.
The "right reader" is someone who already knows good estimates for $|\Gamma(z)|$, who is familiar with the lemma that a harmonic function where $|g(z)| = o(|z|^k)$ is a polynomial of degree $\leq k-1$, and who knows how to compute $\Gamma(1/2)$. I'll try to edit in proofs of these later today. This answer is CW in case someone else wants to do it.
Estimates for $\Gamma(z)$: For $\Re(z)>0$, we have $|\Gamma(z)| \leq \int_0^{\infty} e^{-t} |t^z| dt/t = \int_0^{\infty} e^{-t} t^{\Re(z)} dt/t = \Gamma(\Re(z))$. By Stirling's formula, this shows that $\log |\Gamma(z)| = O(|z| \log |z|)$. Now using the recursion for the $\Gamma$ function lets us extend this estimate to all of $\mathbb{C}$ (details omitted).
We also have the easy estimate $\log |\sin (\pi z)| = O(|z|)$, on a contour which stays well away from the zeroes of $\sin$: Say a circle of radius $N+1/2$. So $\Re(g) = O(|z| \log |z|)$ for $z$ on a circle of radius $N+1/2$. As $g$ is entire, the maximum modulus principle gives us the same bound everywhere in $\mathbb{C}$.
Harmonic functions with slow growth rate: Let $f$ be a harmonic function on the unit disc. We have the Poisson integral formula:
$$f(x+iy) = \int \frac{1-x^2-y^2}{1-2 (x \cos \phi - y \sin \phi) + x^2+y^2} f(e^{i \phi}) d \phi.$$ Differentiating inside the integral sign, there is some smooth function $K(x,y,\phi)$ such that $$\frac{\partial^{a+b} f}{(\partial x)^a (\partial y)^b}(x+iy) = \int K_{ab}(x,y,\phi) f(e^{i \phi}) d \phi.$$
Now, let $f$ be defined on a circle of radius $R$. Making the appropriate variable changes, $$\frac{\partial^{a+b} f}{(\partial x)^a (\partial y)^b}(x+iy) = \frac{1}{R^{a+b+1}} \int K_{ab}(x/R, y/R, \phi) f(R e^{i \phi}) d \phi.$$ So, if $g$ is entire and $g = o(R^k)$, then every $k$-fold derivative of $g$ is $o(1/R)$. Sending $R$ to infinity, every $k$-fold derivative of $g$ is zero, so we deduce that $g$ is a polynomial of degree $\leq k$.
The Gamma function at $1/2$:
We have $\Gamma(1/2) = \int_0^{\infty} e^{-t} t^{-1/2} dt = \int_0^{\infty} e^{-u^2} (2 du)$. This can be evaluated by a variety of methods.
This is an argument offered by Paul Monsky in the comments. Now that I have the details right, I definitely think it is better than mine.
For $x$ between $1$ and $2$, and $y$ positive, we have $$|\Gamma(x+iy)| \leq \Gamma(x) = O(1)$$ by the triangle inequality applied to the defining integral. Using the multiplicative recursion for the $\Gamma$ function, we deduce that $\Gamma(x+iy) = O(1)$ for $x \in [0,1]$ and $y \geq 1$. (An earlier version of this proved that, in fact $\Gamma(x+iy) = O(e^{-y})$ in this range, but I don't need it.)
Setting $f(z) = \Gamma(z) \Gamma(1-z) \sin(\pi z)$, for $x \in [0,1]$ we have $$|f(x+iy)| = O(1) O(1) O(e^{\pi y}) = O(e^{(\pi ) y}) \ \mbox{as}\ y \to\infty$$
Since $f(z) = f(z+1)$, we have $f(z) = F(e^{2 \pi i z})$ for some holomorphic function $F$ defined on $\mathbb{C}^*$. The above bounds show that $|F(w)| = O(|w|^{1/2})$ as $w \to \infty$ and similarly $F(w) = O(|w|^{-1/2})$ as $w \to 0$. Since $F$ is an entire function, this shows that $F$ is constant, and any number of arguments can establish the constant.
It can be shown (from the Beta function) that $$ \Gamma(1-x) \Gamma(x) = \mathrm{B}(x, 1 -x) = \int_0^{\infty} \frac{s^{x-1} d s}{s+1} \label{1}\tag{1} $$ Now we show that $$ \int_0^{\infty} \frac{s^{x-1} d s}{s+1} = \frac{\pi}{\sin \pi x} $$ We use the contour shown in Figure below
The singularities are a pole at $s=-1$ and a branch point at $s=0$ where the function is multivalued. We can use the positive $x$ axis as a branch cut. The contour encloses the pole, so the integral along the contour $C=C_1 \cup C_2 \cup C_3 \cup C_4$ is given by $$ \int_C \frac{s^{x-1} ds}{s+1} = 2 \pi \mathrm{i} (-1)^{x-1} = 2 \pi \mathrm{i} \, \mathrm{e}^{ \mathrm{i} \pi (x-1)}. $$ since $(-1)^{x-1}=\mathrm{e}^{- \mathrm{i} \pi (x-1)}$ is the only residue of the integrand. We now evaluate the four individual integrations for the four paths in the figure. Let us call $I_i = \int_{C_i} s^{x-1}/(1+s) ds$. We start with the integrals along the circular paths. For the small circle we can write $s= \epsilon \mathrm{e}^{\mathrm{i} \theta}$ where $\theta \in[ \delta, 2 \pi - \delta]$ with $\epsilon$ the radius of the disk, and $\delta$ the initial angle of integration. We then change variables from $s$ to $\theta$ with $ds = \mathrm{i} \epsilon \mathrm{e}^{\mathrm{i} \theta}$. That is, $$ |I_4| = \left | \int_ {2 \pi - \delta_{\epsilon}}^ {\delta_{\epsilon}} \frac{\mathrm{i} \epsilon^x \mathrm{e}^{\mathrm{i} \theta (x+1) } d \theta } {\epsilon \mathrm{e}^{\mathrm{i} \theta} -1 } \right | \le \int_{\delta_{\epsilon}}^{2 \pi - \delta_{\epsilon}} \frac{|\epsilon^x| d \theta } {1 - \epsilon } = (2 \pi - 2 \delta_{\epsilon} ) \frac{|\epsilon^x| } {1 - \epsilon } $$ which goes to $0$ as $\epsilon, \delta_{\epsilon} \to 0$, since $0 < x < 1$. Likewise along the big circle $s = R \mathrm{e}^{\mathrm{i} \theta}$ and $ds = \mathrm{i} R \mathrm{e}^{\mathrm{i} \theta} d \theta$, and so $$ |I_2| = \left | \int_{\delta_{R}}^{2 \pi - \delta_{R}} \frac{\mathrm{i} R^x \mathrm{e}^{\mathrm{i} \theta (x+1) } d \theta } {R \mathrm{e}^{\mathrm{i} \theta} -1 } \right | \le \int_{\delta_{R}}^{2 \pi - \delta_{R}} \frac{|R^x| d \theta } {R-1 } = (2 \pi - 2 \delta_{R} ) \frac{R^x } {R-1 } $$ We use the L'Hôpital rule to find that $$ \lim_{R \to \infty} |I_2| = \lim_{R \to \infty} 2 (\pi - \delta_R) \frac{x R^{x-1}}{1} = \lim_{R \to \infty} \frac{x}{R^{1-x}} = 0 $$ since $1 > 1-x > 0$. We are left with the integrals along $C_1$ and $C_3$. If in the integral along $C_1$ we take the limit as $\epsilon \to 0 , R \to \infty$, is the original integral \eqref{1}. On the other hand, the integral over $C_3$ has the argument shifted by $2 \pi$ with respect to the original integral. That is, $$ \lim_{\epsilon \to 0, R \to \infty} \int_{C_3} \frac{s^{x-1} ds}{s+1} = \int_{\infty}^0 \frac{ \mathrm{e}^{2 \pi \mathrm{i} (x-1)} s^{x-1} ds}{\mathrm{e}^{2 \pi \mathrm{i}} s + 1} = -\mathrm{e}^{2 \pi \mathrm{i} ( x-1)} \int_0^{\infty} \frac{s^{x-1} ds }{s+1} $$ Putting all integrals together we find that $$ (1 - \mathrm{e}^{2 \mathrm{i} \pi ( x-1) }) \int_0^{\infty} \frac{s^{x-1} ds}{s+1} = 2 \pi \; \mathrm{i} \; \mathrm{e}^{ \mathrm{i} \pi ( x-1)}. $$ Hence $$ \int_0^{\infty} \frac{s^{x-1} ds}{ s+1} = \frac{2 \pi \mathrm{i} \mathrm{e}^{ \mathrm{i} \pi ( x-1)}}{ 1 - \mathrm{e}^{2 \mathrm{i} \pi(x-1)}} = \frac{\pi}{( -\mathrm{e}^{- \mathrm{i}\pi x} + \mathrm{e}^{\mathrm{i} \pi x})/2 i} = \frac{\pi}{\sin \pi x} $$ We then showed that Euler's reflection formula.
This can be done without invoking the Beta function explicitly. Directly from the integral definition of $\Gamma(x)$,the product in question is a double integral: $$ \Gamma(x)\Gamma(1-x) =\int^\infty_0\int^\infty_0 (u/v)^x u^{-1} e^{-(u+v)} du\, dv $$ Switching to polar coordinates via $v=x^2=r^2\cos^2\phi$, $u=y^2=r^2\sin^2\phi$, we have $$ \Gamma(x)\Gamma(1-x)= 4\int_0^{\pi/2} (\tan\phi)^{2x-1} \, d\phi \int_0^\infty r e^{-r^2}\, dr = 2\int_0^{\pi/2} (\tan\phi)^{2x-1} d\phi $$ Finally, $\phi= \tan^{-1}\sqrt{s}$, brings us to $$ \Gamma(x)\Gamma(1-x)=\int^\infty_0 {s^{x-1} ds\over 1+s} $$ which can be evaluated straightforwardly by contour integration as shown elsewhere on this page.
Apparently, there is no short proof of the reflection formula (unless you count the one you are alluding to, which does need a fair bit of background). There is, however, at least one down-to-earth one, see
I have not seen a one line proof of this statement. I believe there is no such a thing. The statement usually takes some work.
Here is my proof: The formula states that \begin{eqnarray} \Gamma(1 -z ) \Gamma(z) = \frac{\pi}{\sin \pi z} \end{eqnarray} We show this formula using contour ingegration. We start with equation for the Beta function in terms of the $\Gamma$ function (second property) with $y=1-x$, and $0 < x < 1$. That is
\begin{eqnarray*} \Gamma(x) \Gamma(1-x) = \mathrm{B}(x, 1-x). \end{eqnarray*} We now show, by using contour integration, that
\begin{eqnarray*} \mathrm{B}(x, 1-x) = \frac{\pi}{\sin \pi x}. \end{eqnarray*} For this, we use a Beta function representation ($\mathrm{B}=\int_0^{\infty} s^{x-1}/(1+s)^{x+y}$) which for $y=1-x$ becomes
\begin{eqnarray*} \mathrm{B}(x, 1-x) = \int_0^{\infty} \frac{s^{x-1} ds}{s+1} \end{eqnarray*} We use contour integration. Let us first make the substitution $s=e^t$, $ ds = e^t dt$, and $t \in (-\infty, \infty)$. So we need to compute
We compute \begin{eqnarray*} \mathrm{B}(x, 1-x) = \int_{-\infty}^{\infty} \frac{\mathrm{e}^{t(x-1)} \mathrm{e}^t dt}{\mathrm{e}^t+1} = \mathrm{B}(x, 1-x) = \int_{-\infty}^{\infty} \frac{\mathrm{e}^{tx} dt}{\mathrm{e}^t+1} \quad , \quad 0 < x < 1. \end{eqnarray*}
Let us consider the contour integral
\begin{eqnarray} I = \int_C f(z) dz, \end{eqnarray} with \begin{eqnarray*} f(z) = \frac{\mathrm{e}^{zx}}{\mathrm{e}^z+1} \quad , \quad 0 < x < 1. \end{eqnarray*} and $C$ is the contour that we need to determine. In the complex plane, the poles of the integrand are the roots of $e^z+1$, that is $\mathrm{e}^z = -1 = \mathrm{e}^{(2k+1) \mathrm{i} \pi}$ so the roots are $z_k= (2k+1) \mathrm{i} \pi$, for $k=0, \pm 1, \pm 2, \cdots$. Then $f(z)$ as an infinite number of poles all lying on the imaginary axis. We will select a contour that has only one pole as shown in the Figure below
The contour $C$ can be seen as the union of $C=C_1 \cup C_2 \cup C_3 \cup C_4$,
where $C_1$ and $C_3$ are horizontal lines from $-R$ to $R$ with opposite orientation.
We want to let $R$ grow to $\infty$. The paths $C_2$ and $C_4$ are vertical lines
between $0$ and $2 \pi \mathrm{i}$ with opposite orientations showed in the figure.
From the Residue Theorem we evaluate the integral over $C$. The residue corresponding to the pole $z_0= \pi \mathrm{i}$, is computed using the expression
\begin{eqnarray*} \lim_{z \to z_0} (z-z_0) f(z) = \lim_{z \to z_0} \frac{ (z-z_0) \mathrm{e}^{z x}} {\mathrm{e}^z + 1} = \lim_{z \to z_0} \frac{\mathrm{e}^{zx} + (z-z_0) \mathrm{e}^{zx}}{e^z} = \mathrm{e}^{z_0 (x-1)}. \end{eqnarray*} where we use L'H^{o}pital's rule.
Hence $I = 2 \pi \mathrm{i} \; \mathrm{e}^{\pi i (x-1)}$ since the only residue inside the contour is at $z=\mathrm{i} \pi$ . That is,
\begin{eqnarray*} 2 \pi \; \mathrm{i} \; \mathrm{e}^{ \mathrm{i} \pi (x-1)} =\int_{C_1} f(z) dz + \int_{C_2} f(z) dz + \int_{C_3} f(z) dz + \int_{C_4} f(z) dz . \end{eqnarray*}
We want to find $I_1=\int_{C_1} f(z) dz$, as $R \to \infty$. Let us first find the integral along the vertical path $C_3$.
\begin{eqnarray*} I_3 &=& \int_R^{-R} \frac{\mathrm{e}^{ (t + 2 \pi \mathrm{i} ) x}}{\mathrm{e}^{t+2 \pi \mathrm{i}} + 1 } d t \\ &=& \mathrm{e}^{2 \pi \mathrm{i} x} \int_R^{-R} \frac{\mathrm{e}^{ t x}}{\mathrm{e}^{t+2 \pi \mathrm{i}} + 1 } d t \\ &=& \mathrm{e}^{2 \pi \mathrm{i} x} \int_R^{-R} \frac{\mathrm{e}^{ t x}}{\mathrm{e}^{t} + 1 } d t \\ &=& -\mathrm{e}^{2 \pi \mathrm{i} x} I_1, \end{eqnarray*} where we reversed the sign since $I_1$ is computed from $-R$ to $R$ instead of going in the opposite direction.
The integral $I_2$ along the path $C_2$ is evaluated as follows \begin{eqnarray*} I_2 = \int_0^{2 \pi} \frac{\mathrm{e}^{ (R + \mathrm{i} t ) x}} {\mathrm{e}^{R+\mathrm{i} t} + 1 } d t = \frac{\mathrm{e}^{R x} }{\mathrm{e}^R} \int_0^{2 \pi} \frac{\mathrm{e}^{ (\mathrm{i} t ) x}} {\mathrm{e}^{\mathrm{i} t} + 1/\mathrm{e}^{R} } d t = \mathrm{e}^{R(x-1)} \int_0^{2 \pi} \frac{\mathrm{e}^{ (\mathrm{i} t ) x}} {\mathrm{e}^{\mathrm{i} t} + 1/\mathrm{e}^{R} } d t. \end{eqnarray*} Now, since $0 < x < 1$ (so $x-1 < 0)$, and the last integral is bounded we have that $\lim_{R \to \infty} I_2 = 0$. The same argument applies for the integral $I_4$ along the path $C_4$. We then have that, from $I=I_1+I_2+I_3+I_4$,
\begin{eqnarray*} 2 \pi \mathrm{i} \mathrm{e}^{\mathrm{i} \pi (x-1)} = (1 - \mathrm{e}^{2 \pi \mathrm{i} x}) I_1 \end{eqnarray*} and
\begin{eqnarray*} I_1 = \int_{0}^{\infty} \frac{s^{x-1} ds}{s+1} = \frac{ 2 \pi \mathrm{i} \mathrm{e}^{ \mathrm{i} \pi (x-1)}}{1 - \mathrm{e}^{2 \pi \mathrm{i} x}} = \frac{ \pi}{ \frac{\mathrm{e}^{\mathrm{i} \pi x} - \mathrm{e}^{-\mathrm{i} \pi x}}{2 \mathrm{i}}} = \frac{\pi}{\sin \pi x}. \end{eqnarray*}
We then showed that Euler's reflection formula is correct.
