Proof involving gamma function, infinite product and Gauss How can I rigorously and directly prove that $$\Gamma (z)=\lim_{n\rightarrow \infty }\frac{n!n^{z}}{z(z+1)\cdots(z+n)}$$
For simplicity define $z=k$ ,so we prove

$$\lim_{n\to \infty}\frac{n! \,\, n^k}{k(k+1)\cdots (k+n)}=(k-1)!$$

Note that

$$\frac{(k+n)!}{(k-1)!} = k(k+1)\cdots (k+n)$$ 

so we have 

$$ (k-1)! \lim_{n\to \infty}\frac{ \, n! \,\, n^k}{(k+n)!}=(k-1)!\lim_{n\to \infty}\frac{n^k}{(k+n)(k+n-1)\cdots (n+1)}$$

Now this can be written as 

$$(k-1)! \lim_{n\to \infty} \frac{n\cdot n\cdot \cdots n}{(n+k)(n+k-1)\cdots (n+1)}=(k-1)!$$

The gamma function is an extension of the factorial by the relation 

$$\Gamma(k)=(k-1)!$$
We define the function

$$\psi_n(t):=\begin{cases}
\left(1-\dfrac{t}{n}\right)^n t^{x-1} &amp; \text{if}\;0\le x\le n\\
0 &amp; \text{if}\;x&gt;n\\\end{cases}.$$

First we prove: for any $\varepsilon&gt;0$ and any $t_0&gt;0$ there exists a $N_0\in\mathbb N$ such that

$$\left\vert e^{-t} t^{x-1}-\psi_n(t)\right\vert&lt;\varepsilon$$

for $n\ge N_0$ and $t\in[0,t_0]$. With 

$$\chi_n(k):=\begin{cases}
1 &amp; \text{if}\;0\le k\le n\\
0 &amp; \text{if}\;k&gt;n\\\end{cases}$$

we establish

$${\left\vert e^{-t} t^{x-1}-\psi_n(t)\right\vert}=
{\left\vert t^{x-1}\left\{ e^{-t} -\left(1-\dfrac{t}{n}\right)^n\right\}\right\vert}=
{\left\vert t^{x-1}\left\{ e^{-t} -\sum_{k=0}^{n}(-1)^k\binom{n}{k}\left(\dfrac{t}{n}\right)^k\right\}\right\vert}=
{\left\vert t^{x-1}\left\{ e^{-t} -\sum_{k=0}^{n}\dfrac{(-1)^k}{k!}t^k\left(1-\dfrac{1}{n}\right)\left(1-\dfrac{2}{n}\right)\cdots\left(1-\dfrac{k-1}{n}\right)\right\}\right\vert}=
{\left\vert t^{x-1}\left\{ \sum_{k=0}^{\infty}\dfrac{(-1)^k}{k!}t^k -\sum_{k=0}^{\infty}\dfrac{(-1)^k}{k!}t^k\left(1-\dfrac{1}{n}\right)\left(1-\dfrac{2}{n}\right)\cdots\left(1-\dfrac{k-1}{n}\right)\chi_n(k)\right\}\right\vert}=
{\left\vert t^{x}\sum_{k=1}^{\infty}\dfrac{(-1)^k}{k!}t^{k-1}\left\{1-\left(1-\dfrac{1}{n}\right)\left(1-\dfrac{2}{n}\right)\cdots\left(1-\dfrac{k-1}{n}\right)\chi_n(k)\right\}\right\vert}\le
{t^{x}\sum_{k=1}^{\infty}\dfrac{t^{k-1}}{k!}\left\{1-\left(1-\dfrac{1}{n}\right)\left(1-\dfrac{2}{n}\right)\cdots\left(1-\dfrac{k-1}{n}\right)\chi_n(k)\right\}}.$$

With 

$$b_k:=\dfrac{1-\left(1-\dfrac{1}{n}\right)\left(1-\dfrac{2}{n}\right)\cdots\left(1-\dfrac{k-1}{n}\right)\chi_n(k)}{k!}$$

we have $0\le b_k\le \dfrac{1}{k!}$.  Therefore the series has the majorant $\sum_k\dfrac{t^{k-1}}{k!}$. The radius of convergence of this  series is $R=\infty$ with the quotient criterion

$$\dfrac{\left(\dfrac{t^k}{(k+1)!}\right)}{\left(\dfrac{t^{k-1}}{k!}\right)}=\dfrac{t}{k+1}.$$

Hence there exists a $K_1\in\mathbb N$ such that

$$\tag{1}0\le{t^{x}\sum_{k=K_1+1}^{\infty}\dfrac{t^{k-1}}{k!}\left\{1-\left(1-\dfrac{1}{n}\right)\left(1-\dfrac{2}{n}\right)\cdots\left(1-\dfrac{k-1}{n}\right)\chi_n(k)\right\}}&lt;\dfrac{\varepsilon}{2}$$

for $t\in[0,t_0]$. With $\zeta=\dfrac{1}{n}$ and $k\le n$ we write  for the coefficient $b_k$ above

$$b_k=b_k(\zeta)=\zeta\cdot p_k(\zeta)$$

with a polynompial $p_k(\zeta)$. We obtain for $n\ge K_1$

$$t^{x}\sum_{k=1}^{K_1}b_kt^{k-1}=\zeta\cdot t^xp(t,\zeta)$$

with another polynomial $p(t,\zeta)$. On the closed set $[0,t_0]\times[0,1]$ the function $\vert t^xp(t,\zeta)\vert$ has a maximum so that we get

$$\left\vert t^xp(t,\zeta)\right\vert\le B$$

for some $B&gt;0$. Now we only have to choose a $N_0\in\mathbb N$, $N_0\ge K_1$, such that

$$\tag{2}0\le t^{x}\sum_{k=1}^{K_1}b_kt^{k-1}=\dfrac{1}{n}\cdot t^xp\left(t,\dfrac{1}{n}\right)\le\dfrac{B}{n}\le\dfrac{\varepsilon}{2}$$

for $n\ge N_0$ and $t\in[0,t_0]$. The inequalities (1) and (2) prove now

$$\tag{3}{\left\vert e^{-t} t^{x-1}-\psi_n(t)\right\vert}&lt;\varepsilon$$

for $n\ge N_0$ and $t\in[0,t_0]$.

For $t&gt;0$ we get from

$${\lim_{n\to\infty}\left\vert e^{-t} t^{x-1}-\psi_n(t)\right\vert}\to 0$$

the well-known limit

$$\tag{4}\lim_{n\to\infty}\left(1-\dfrac{t}{n}\right)^n=e^{-t}$$

that holds for $t=0$ too. Next we prove for $t\in[0,n)$

$$\tag{5}0\le\left(1-\dfrac{t}{n}\right)^n\le\left(1-\dfrac{t}{n+1}\right)^{n+1}.$$

The first inequality is obvious. The derivitave of the function $h(z):=\left(1-\dfrac{t}{z}\right)^z$ is

$$h'(z)={\dfrac{d}{dz}e^{z \ln\left(1-\dfrac{t}{z}\right)}}=
\left\{\ln\left(1-\dfrac{t}{z}\right)+\dfrac{\dfrac{t}{z}}{1-\dfrac{t}{z}}\right\}\left(1-\dfrac{t}{z}\right)^z.$$

With

$${\ln\left(1-\dfrac{t}{z}\right)+\dfrac{\dfrac{t}{z}}{1-\dfrac{t}{z}}}=
-\sum_{k=1}^\infty\dfrac{1}{k}\left(\dfrac{t}{z}\right)^k+\sum_{k=1}^\infty\left(\dfrac{t}{z}\right)^k$$

we see that $h'(z)\ge0$ for $0\le t&lt;z$ and the inequality (5) is proven. Because the function $\psi_n(t)$ is not negative, grows with n and converges to $e^{-t}t^{x-1}$ we obtain

$$\tag{6}0\le{\int_0^\infty\left\{e^{-t}t^{x-1}-\psi_n(t)\right\}dt}\le
{\int_0^\infty e^{-t}t^{x-1}dt}=:\Gamma(x).$$

Because the integral on the right side converges, for any $\epsilon&gt;0$ there exists a $T_0&gt;0$ such that 

$$0\le{\int_{T_0}^\infty\left\{e^{-t}t^{x-1}-\psi_n(t)\right\}dt}\le
{\int_{T_0}^\infty e^{-t}t^{x-1}dt}&lt;\dfrac{\varepsilon}{2}.$$

We choose a $N_1\in\mathbb N$ such that 

$${\left\vert e^{-t} t^{x-1}-\psi_n(t)\right\vert}\le\dfrac{\varepsilon}{2T_0}$$

for $n\ge N_1$ in the inequality (3). The last $2$ inequalities render for $n\ge N_1$

$$0\le{\int_0^\infty\left\{e^{-t}t^{x-1}-\psi_n(t)\right\}dt}=
{\int_0^{T_0}\left\{e^{-t}t^{x-1}-\psi_n(t)\right\}dt+\int_{T_0}^\infty\left\{e^{-t}t^{x-1}-\psi_n(t)\right\}dt}\le
\int_0^{T_0}\dfrac{\varepsilon}{2T_0}dt+\dfrac{\varepsilon}{2}\le\varepsilon$$

or

$$\Gamma(x)={\lim_{n\to\infty}\int_0^\infty\psi_n(t)dt}=
{\lim_{n\to\infty}\int_0^n\left(1-\dfrac{t}{n}\right)^n t^{x-1}dt}.$$

You get your formula by integration by parts:

$${\int_\alpha^n\left(1-\dfrac{t}{n}\right)^n t^{x-1}dt}=
{\left[\left(1-\dfrac{t}{n}\right)^n \dfrac{t^x}{x}\right]_\alpha^n}+{\dfrac{1}{x}\int_\alpha^n\left(1-\dfrac{t}{n}\right)^{n-1} t^{x}dt}.$$

For $\alpha\to 0+$ we obtain

$${\int_0^n\left(1-\dfrac{t}{n}\right)^n t^{x-1}dt}=
{\dfrac{1}{x}\int_0^n\left(1-\dfrac{t}{n}\right)^{n-1} t^{x}dt}.$$

Consecutively integrating by parts renders

$${\int_0^n\left(1-\dfrac{t}{n}\right)^n t^{x-1}dt}=
{\dfrac{1}{x}\int_0^n\left(1-\dfrac{t}{n}\right)^{n-1} t^{x}dt}=
{\dfrac{n-1}{n}\cdot\dfrac{1}{x(x+1)}\int_0^n\left(1-\dfrac{t}{n}\right)^{n-2} t^{x+1}dt}=\cdots=
{\dfrac{(n-1)(n-2)\cdots 1}{n^{n-1}}\cdot\dfrac{1}{x(x+1)\cdots(x+n-1)}\int_0^nt^{x+n-1}dt}={\dfrac{n-1!}{n^{n-1}}\cdot\dfrac{1}{x(x+1)\cdots(x+n-1)}\cdot\dfrac{n^{x+n}}{x+n}}=
{\dfrac{n!n^x}{x(x+1)\cdots(x+n)}}.$$

In an alternative proof we could also have taken the well-known limit (4) in the case $x\ge1$. For any $\varepsilon&gt;$ and $t\in[0,t_0]$ we would have determined a $N_0(t_0)\in
\mathbb N$ such that 

$$\left\vert \left(1-\dfrac{t_0}{N_0(t_0)}\right)^{N_0(t_0)} t_0^{x-1}-e^{-t_0}t_0^{x-1} \right\vert\le\dfrac{\varepsilon}{2}.$$

Because the functions $\left(1-\dfrac{t}{N_0(t_0)}\right)^{N_0(t_0)} t^{x-1}$ and $e^{-t}t^{x-1}$ are continuous there exists an open interval $\mathbb M(t_0)$ around the point $t_0$ such that 

$$\tag{7}\left\vert\left(1-\dfrac{t}{N_0(t_0)}\right)^{N_0(t_0)} t^{x-1}-e^{-t}t^{x-1}\right\vert\le\varepsilon$$

for $t\in\mathbb M(t_0)$. Because the function $\left(1-\dfrac{t}{n}\right)^n t^{x-1}$ grows with n, the inequality (7) is true then for any $n\ge N_0(t_0)$. For an interval $[0,t_0]$ there exists a finite number of open intervals $M(t_k)$ such that $[0,T_0]\subset\bigcup_k\mathbb M(t_k)$ with the theorem of Heine-Borel. We determine for such a finite choice

$$N_0=\max_k N_0(t_0)$$

and obtain the inequality once again. However for $0&lt;x&lt;1$ we have

$$\lim_{t\to 0+}e^{-t}t^{x-1}\to +\infty$$

and the inequality is not immediate for $t$ about $0$. Hence in this case it is not enough that the function $\psi_n(t)$ is growing and converging to the function $e^{-t}t^{x-1}$. However we can find a $\delta&gt;0$ such that

$$0\le\int_0^\delta\left\{e^{-t}t^{x-1}-\psi_n(t)\right\}dt\le
\int_0^\delta e^{-t}t^{x-1}dt\le\varepsilon,$$

and then we restrict our attention to the interval $[\delta,T_0]$.

For calculating the gamma function you better go to Mathmatics:Algorithm to compute Gamma function. The product formula of the function $\Gamma(x)$ is mostly used in order to prove other formulas like

$$\dfrac{\pi}{\sin\pi x}=\Gamma(x)\Gamma(1-x)$$

or

$$\dfrac{\sqrt{\pi}}{2^{2x-1}}\Gamma(2x)=\Gamma(x)\Gamma(x+\dfrac{1}{2}).$$