The following exercise is a recap on probability and maths for statistical mechanics:
Maximize $$f(n_1,n_2,\, ...,\, n_M)=\dfrac{N!}{\prod_{j=1}^M n_j!}$$ subject to $$\sum_{j=1}^M n_j=N\ \text{ and }\ \sum_{j=1}^M e_jn_j=E,$$ with $e_j$ and $E$ constants.
Physical background: This problem can be understood as a problem of $M$ states and $N$ independent particles. The aim is to determine what is the distribution of particles among the different states that maximizes the number of microstates compatible with fixed values of $N$ and $E$.
Attempt:
I should impose $$\frac{\partial}{\partial n_i}(f(\mathbf n)-\lambda A(\mathbf n)-\beta B(\mathbf n))=0,\text{ for } i=\{1,\,...,M\},$$ where $A(\mathbf n)\equiv\sum_{j=1}^M n_j-N$ and $B(\mathbf n)\equiv \sum_{j=1}^M e_jn_j-E$, and so, $$\frac{\partial}{\partial n_i}\left(\dfrac{N!}{\prod_{j=1}^M n_j!}-\sum_{j=1}^M(\lambda+\beta e_j)n_j -\lambda N-\beta E\right)=0$$ $$\dfrac{N!n_i!}{\prod_{j=1}^M n_j!}\dfrac{\partial}{\partial n_i}\left(\frac{1}{n_i!}\right)-(\lambda+\beta e_i)=0,$$ but now I'm stuck here as I don't really know how to calculate the derivative of $1/n_i!$ or even how to proceed once I do it because the derivative should include digamma or gamma functions...