Solution (source code)

= Solution

If exactly $k$ spins equal $+1$, then $M=(2k-n)/n$. Therefore its possible values are $E_n$, with $k=n(1+m)/2$ for $m\in E_n$. Choosing the $k$ positive spins gives the <binomial coefficient>
$$
\binom nk=\frac{n!}{[n(1+m)/2]!\,[n(1-m)/2]!}.
$$
Every such configuration has the same <Curie–Weiss model> weight because $\sum_{i,j}s_is_j=(\sum_is_i)^2=n^2m^2$. The <probability mass function> of this <spin magnetization> is consequently
$$
\boxed{\mathbb P(M=m)=\frac1{Z_{n,\beta}}
\binom{n}{n(1+m)/2}e^{\beta nm^2/2},\qquad m\in E_n,}
$$
and zero elsewhere. Grouping the original <partition function> by the same count also gives
$$
Z_{n,\beta}=\sum_{k=0}^n\binom nk
\exp\!\left[\frac{\beta(2k-n)^2}{2n}\right],
$$
which verifies normalization without counting any configuration twice.