One-line proof of the Euler's reflection formula 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?
 A: 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.
A: It can be shown (from the Beta function) that
\begin{eqnarray}
  \Gamma(1-x) \Gamma(x) = \mathrm{B}(x, 1 -x) 
  = \int_0^{\infty} \frac{s^{x-1} d s}{s+1} \quad \quad (1)
\end{eqnarray}
Now we show that
\begin{eqnarray*}
  \int_0^{\infty} \frac{s^{x-1} d s}{s+1} =  \frac{\pi}{\sin \pi x}
\end{eqnarray*}
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
\begin{eqnarray*}
  \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)}.
\end{eqnarray*}
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,
\begin{eqnarray*}
  |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 } 
\end{eqnarray*}
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
\begin{eqnarray*}
  |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 } 
\end{eqnarray*}
We use the L'H^{o}pital rule to find that
\begin{eqnarray*}
  \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
\end{eqnarray*}
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 (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,
\begin{eqnarray*}
  \lim_{\epsilon \to 0, R \to \infty} \int_{C_3} \frac{s^{x-1} ds}{(s+1)^{x+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)^{x+1}} = -\mathrm{e}^{2 \pi \mathrm{i} ( x-1)} \int_0^{\infty} 
    \frac{s^{x-1} ds }{(s+1)^{x+1}}
\end{eqnarray*}
Putting all integrals together we find that
\begin{eqnarray*}
  (1 - \mathrm{e}^{2 \mathrm{i} \pi ( x-1) }) 
  \int_0^{\infty} \frac{s^{x-1} ds}{(s+1)^{x+1}} = 
  2 \pi \; \mathrm{i} \; \mathrm{e}^{ \mathrm{i} \pi ( x-1)}.
\end{eqnarray*}
Hence
\begin{eqnarray*}
 \int_0^{\infty} \frac{s^{x-1} ds}{ (s+1)^{x+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}^{- \pi x} + \mathrm{e}^{\mathrm{i} \pi x})/2 i}
   = \frac{\pi}{\sin \pi x}
\end{eqnarray*}
We then showed that Euler's reflection formula
A: 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.
A: 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
http://warwickmaths.org/files/gamma.pdf
A: 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.
A: 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. 
