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.

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.

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.

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.