Another common approach is to derive it from the limit definition of the gamma function. (See below.)
The multiplication formula can be written in the form
$$n^{nz-1/2} \prod_{k=0}^{n-1} \Gamma \left(z+\frac{k}{n} \right) = (\sqrt{2 \pi})^{n-1} \Gamma(nz)$$
Using the limit definition of the gamma function, we have
$$ \Gamma \left(z +\frac{k}{n} \right) = \lim_{m \to \infty} \frac{m! \ m^{z+\frac{k}{n}-1}}{(z+\frac{k}{n})(z+\frac{k}{n}+1) \cdots (z+\frac{k}{n} + m -1)}$$
Then using Stirling's formula, we get
$$ \begin{align} \Gamma \left(z+\frac{k}{n} \right) &= \lim_{m \to \infty} \frac{\sqrt{2 \pi m} (\frac{m}{e})^m m^{z+\frac{k}{n}-1}}{(z+\frac{k}{n})(z+\frac{k}{n}+1) \cdots (z+\frac{k}{n} + m -1)} \\ &=\lim_{m \to \infty} \frac{\sqrt{2 \pi} (\frac{mn}{e})^m m^{z+\frac{k}{n}-1/2}}{(nz+k)(nz+k+n) \cdots (nz+k + mn -n)} \end{align}$$
So
$$ n^{nz-1/2} \prod_{k=0}^{n-1} \Gamma \left(z+\frac{k}{n} \right) $$
$$ = n^{nz-1/2}\lim_{m \to \infty}\frac{(\sqrt{2 \pi})^{n} (\frac{mn}{e})^{mn} m^{nz-n/2} m^{\frac{1}{n} \sum_{k=1}^{n-1} k}}{(nz)(nz+1)\cdots (nz+n-1)(nz+n) \cdots (nz+mn-n)\cdots(nz+mn-1)} $$
$$ = \lim_{m \to \infty} \frac{(\sqrt{2 \pi})^{n} (\frac{mn}{e})^{mn} (mn)^{nz-1/2}}{(nz)(nz+1)\cdots (nz+n-1)(nz+n) \cdots (nz+mn-n) \cdots(nz+mn-1)} $$
Replacing $mn$ with $m$ shouldn't change the value of the limit (I think).
Therefore,
$$ \begin{align} n^{nz-1/2} \prod_{k=0}^{n-1} \Gamma \left(z+\frac{k}{n} \right) &=\lim_{m \to \infty} \frac{(\sqrt{2 \pi})^{n} (\frac{m}{e})^{m} m^{nz-1/2}}{(nz)(nz+1)\cdots(nz+m-1)} \\ &=\lim_{m \to \infty} \frac{(\sqrt{2 \pi})^{n-1}m! \ m^{nz-1}}{(nz)(nz+1)\cdots(nz+m-1)} \\ &= (\sqrt{2\pi})^{n-1}\Gamma(nz) \end{align}$$
EDIT:
Wikipedia states that the limit definition is
$$ \Gamma(t) = \lim_{n \to \infty} \frac{n! \ n^{t}}{t(t+1) \cdots (t+n)}$$
But notice that
$$ \Gamma(t-1) = \frac{\Gamma(t)}{t-1} = \lim_{n \to \infty} \frac{n! \ n^{t-1}}{(t-1)t \ldots (t+n-1)}$$
$$\implies \Gamma(t) = \lim_{n \to \infty} \frac{n! \ n^{t-1}}{t(t+1) \ldots (t+n-1)}$$
Best Answer
We have
\begin{align*} |S|^2 &= \sum_{k=0}^{m-1}\sum_{l=0}^{m-1} \exp\left(\frac{2\pi i}{m}(l^2-k^2)\right) \\ &= \sum_{k=0}^{m-1}\sum_{d=0}^{m-1} \exp\left(\frac{2\pi i}{m}(2kd+d^2)\right) \tag{$l\equiv k+d \pmod{m}$} \\ &= \sum_{d=0}^{m-1} \Biggl[ \sum_{k=0}^{m-1} \exp\left(\frac{4\pi i d}{m} k\right) \Biggr] \exp\left(\frac{2\pi i d^2}{m}\right). \end{align*}
As pointed out by @Cade Reinberger, the geometric sum formula shows that the inner sum only survives with the value $m$ when $d=0$, and hence the desired claim follows.