To calculate sums (or integrals, if the spins are continuous-valued) over such a large number of degrees of freedom, we resort to Monte Carlo techniques.
As with the Monte Carlo integration mentioned earlier, we could just
generate configurations at random, and approximate the real thermal
averages, such as 
, by Monte Carlo averages:
The problem here is that because of the rapidly varying exponential function in the Boltzmann distribution, most randomly chosen configurations will make a negligible contribution to the sum, since E will be relatively large.