I want to know about the generating function of the Riemann zeta function which is related with the Laurent expansion at $z=0$.
$f(z) := \dfrac{d}{dz} \log(\sin\pi z)$
$f(z) = \dfrac{1}{z} -2\sum \zeta(2m+2)z^{2m+1}$
I want to know how to prove the above.
Assuming Euler's product for $\sin$:
$$\sin \pi z= \pi z \prod_{n=1}^\infty \left(1-\left(\frac zn\right)^2\right)$$
Take the log and you get:
$$\log\sin\pi z = \log \pi +\log z + \sum_{n} \log\left(1-\left(\frac zn\right)^2\right)$$
Taking the derivative termwise, we get:
$$\begin{align} f(z)&=\frac{1}{z} + \sum_{n=1}^\infty \left(\frac{-2z}{n^2}\right)\frac{1}{1-\left(\frac zn\right)^2}\\ &=\frac{1}{z}-2\sum_{n=1}^\infty\sum_{m=0}^\infty \frac{z^{2m+1}}{n^{2m+2}}\\ &=\frac{1}{z} -2\sum_{m=0}^\infty z^{2m+1}\sum_{n=1}^\infty \frac{1}{n^{2m+2}}\\ &=\frac{1}{z}-2\sum_{m=0}^\infty \zeta(2m+2)z^{2m+1} \end{align}$$
It all really comes from $$\frac{1}{1-w^2}=\sum_{m=0}^\infty w^{2m},$$ applied in each case to $w=\frac{z}{n}$, then switching the orders of the sums. We can switch the orders of the sums by showing that the sum is absolutely convergent (when $|z|<1$, at least.) We also needed absolute convergence to do the derivatives term-wise, and, of course $|w|<1$ for the power series substitutions.
Aside: Your $f(z)=\pi \cot \pi z$, and the formula only works for $|z|<1$ since $\zeta(2m+2)>1$ so the right hand side does not converge when $|z|\geq 1$.
There are some elementary proofs of Euler's product formula in this question. They are not pretty. Easier to prove with more advanced techniques like complex analysis or Fourier series.