The state $\rho$ for a quantum mechanical Hamiltonian $\hat{H}$ at temperature $T = \frac{1}{k_B\beta}$ is:
So let's consider the cost of a matrix exponential (what you call "exact diagonalization"). This is a very active field of research even today. The algorithm MATLAB R2016b uses comes from a 2010 thesis and 2016 has seen at least 3 new algorithms (Ruiz et al. 2016, Grebrimedhin et al. 2016, Guttel et al., 2016). The choice of algorithm, and the cost of that algorithm depends on the properties of $H$ and on how accurately you are aiming to get the answer, but $15n^3$ floating point operations (FLOPs) has been given as a crude estimate here, where $n$ is the dimension of the matrix, which for a Hamiltonian is $M^2$ where $M$ is the number of levels included in your quantum system.
- $M=2$ for a two-level-system such as the spin of an electron,
- $M=39$ for the vibrational levels of the ground electronic state of $^{6,6}$Li$_2$,
- $M=2^Q$ for $Q$ qubits or $Q$ spin-1/2 particles,
- $M$ is countably infinitely large for a quantum harmonic oscillator,
- $M$ is uncountably infinitely large for a continuous-variable system such as position $\hat{x}$.
So a rough estimate of the deterministic algorithm cost is $15M^6$ FLOPs for an $M$-level system, which you are right, is exponential with respect to the number of spin-1/2 particles.
As for Monte Carlo methods, you will have to be more specific about what you are after. The number of macro-iterations required to obtain your answer with a precision of $\pm\epsilon$ is $\mathcal{O}\left(\sqrt{1/\epsilon}\right)$ and each of these macro-iterations will depend on $M$, but to give you the FLOP count you will have to be more specific in your question.