Here is another way to prove it. Surprisingly, for $n$ integral and $k$ real, the integral in question can be written down as an indefinite integral. This gives a direct proof for non-integer $k$, though obviously less clear than the contour method. (In fact, it is convenient to avoid integer $k$ in this method, and extend to integer $k$ by continuity.)

Writing $y=x/n$ and $l=n-k$, we have for example for $n=2$:
$$\int\frac{\sin^2(2y)}{\sin^k(ky)\sin^l(ly)}dy=
\frac{\frac{2}{l-k}\sin((l-k)y)+\sin(2y)}
{kl\sin^{k-1}(ky)\sin^{l-1}(ly)}$$

In general for $n$ even ($n$ odd is similar with cosines):
$$I_{n,k}(y)=\int\frac{\sin^n(ny)}{\sin^k(ky)\sin^l(ly)}dy=
\frac{\sum_{r=0}^{n-1}\sum_{s=0}^{n-1}\lambda_{r,s}\sin(((n-1-2r)k+(n-1-2s)l)y)}
{kl\sin^{k-1}(ky)\sin^{l-1}(ly)}$$

where
$$\lambda_{r,s}=\begin{cases}
\displaystyle
\frac{(-1)^r(n-1)^{\underline{r}}}{(r-l)^{\underline{\smash{r-s}}}s!}\lambda_{0,0},\;\;r\ge s\\
\displaystyle
\frac{(-1)^s(n-1)^{\underline{s}}}{(s-k)^{\underline{\smash{s-r}}}r!}\lambda_{0,0},\;\;s\ge r
\end{cases}
$$
$$\lambda_{0,0}=(-1)^{n/2+1}2^{1-n},$$
and $x^{\underline{r}}$ denotes the falling power $x(x-1)\ldots(x-r+1)$.

It is easy to check the derivative, $I_{n,k}'(y)$ is correct by considering the coefficient
of $\cos((ak+bl)y)/(\sin^k(ky)\sin^l(ly))$ for each $a$, $b$. If $a\neq b$ then you get zero,
otherwise for $a=b=n-2r$, $(ak+bl)y=(n-2r)ny$ and you get $\frac{1}{2}(-1)^r\binom{n}{r}\lambda_{0,0}\cos((n-2r)ny)$. Then $\sin^n(ny)$ arises from the binomial expansion:
$$(-1)^{n/2}2^{-n}\sum_{r=0}^n (-1)^r\binom{n}{r}\cos((n-2r)ny)=\sin^n(ny).$$

Note that $I_{n,k}(0)=0$ because, being the integral of something well-behaved at $0$,
$I_{n,k}(y)$ must be continuous at $0$, so its numerator must vanish to order $n-2$
like its denominator. Using L'H\^{o}pital, taking $n-2$ derivatives of the numerator gives only sines, which themselves vanish at 0.
To evaluate $I_{n,k}(\pi/n)$, note that $\sin(k\pi/n)=\sin(l\pi/n)$ and
$\sin(((n-1-2r)k+(n-1-2s)l)y)=\sin(2(r-s)k\pi/n)$.
Conditioning on $r-s=d>0$, you get (e.g., by considering partial fractions in $k$)
$$\sum_{s=0}^{n-1-d}\lambda_{s+d,s}=
\frac{(-1)^d(n-1)!\binom{n-2}{d-1}}{(k-1)^{\underline{\smash{n-1}}}}\lambda_{0,0},$$
and similarly for $d<0$ with the opposite sign, and using $-d$ in place of $d$.
Using the binomial expansion of $(1-e^{2\pi ik/n})^{n-2}$, you get
$$\sum_{d=1}^{n-1}(-1)^d\binom{n-2}{d-1}\sin\left(\frac{2dk\pi}{n}\right)=
(-1)^{n/2}2^{n-2}\sin^{n-2}\left(\frac{k\pi}{n}\right)\sin(k\pi)$$
So finally, putting the pieces together,
$$\frac{n}{\pi}I_{n,k}\left(\frac{\pi}{n}\right)=\frac{n\sum_{r=0}^{n-1}\sum_{s=0}^{n-1}\lambda_{r,s}\sin\left(\frac{2(r-s)k\pi}{n}\right)}
{\pi kl\sin^{k-1}\left(\frac{k\pi}{n}\right)\sin^{l-1}\left(\frac{l\pi}{n}\right)}=\frac{n!\sin(k\pi)}{\pi k^{\underline{\smash{n+1}}}}$$
which (for even $n$) we recognise as $\binom{n}{k}$ by the reflection formula for factorials.

21more comments