Residue of the Fermi Distribution Function

The easiest way to see this is to assume that the poles $z_n$ are order one, at which point you can find the residue via $$ R = lim_{z rightarrow z_n} (z - z_n) f(z) = lim_{z rightarrow z_n} frac{z - z_n}{e^{beta z} + 1} = lim_{z rightarrow z_n} frac{1}{beta e^{beta z}} = - frac{1}{beta} $$ where I've used L'Hopital's rule along with $e^{beta z_n} = -1$ for fermionic Matsubara frequencies. More generally in complex analysis, if you had an $n$th order pole located at $z = c$ instead of a 1st order pole, then the above limit will diverge. You could then calculate the residue using $$ R = frac{1}{(n-1)!} lim_{z rightarrow c} frac{d^{n-1}}{dz^{n-1}} [(z - c)^n f(z)] $$ In practice, you could either find an explicit expression for $f(z)$ as a Laurent series and just read off the value of the residue, or you can just keep trying the above formula for $n = 1, 2, ldots$ until one of them doesn't diverge. Often you can stare at a function and get a reasonable feeling for whether the pole will be first order, second order, or something crazier.