Solution (source code)

= Solution

Stirling and replacing the magnetization sum by an <integral> gives the stated form with $P(m)=-\beta(Bm+3Jm^2+2Km^3)$, up to subexponential constants.