I'll start with a few definitions: $$\beta \equiv \frac{1}{k_bT}$$
Where T is the temperature of a system.
And the partition function:
$$Z \equiv \sum_{j}e^{-\beta \epsilon_j}=\int D(\epsilon)d\epsilon e^{-\beta\epsilon}$$
Where s is an index over all distinguishable microstates of the system (in the case of discrete states) and $\epsilon$ is the energy of each of the microstates. In the continuous case, $D(\epsilon)$ is the density of states.
So, if we consider the partition function of a photon gas of temperature $T$ where the energy of a photon w/ frequency $\omega$ is $\epsilon = \hbar\omega$, the total energy of an assemblage of photons is given by:
$$E = \sum_i s_i \hbar\omega_i$$
Where i is the index over all possible frequencies $\omega_i$ and $s_i$ is the number of photons with frequency $\omega_i$.
We can then see that the partition function for a photon gas can be expressed as:
$$Z = \sum_{k}e^{-\beta E_k} = \sum_{s_1,s_2,...=0}^{\infty} e^{-\beta\hbar(s_1\omega_1 + s_2\omega_2 +...)} = \sum_{s_1,s_2,...=0}^{\infty} e^{-\beta\hbar s_1\omega_1}e^{ -\beta \hbar s_2\omega_2}...$$
The second equivalence is established by noting that each state $\Psi_k$ is uniquely determined by the number of photons in each mode $M_i$ (each of which has a corresponding frequency $\omega_i$).
$$Z = \sum_{s_1=0}^{\infty}\sum_{s_2=0}^{\infty}... (e^{-\beta\hbar s_1\omega_1}e^{ -\beta \hbar s_2\omega_2}...) = \sum_{s_1=0}^{\infty}e^{-\beta\hbar s_1 \omega_1}\sum_{s_2=0}^{\infty}\sum_{s_3=0}^{\infty}... (e^{-\beta\hbar s_2\omega_2}e^{ -\beta \hbar s_3\omega_3}...)$$
By pulling out the $\omega_1$ term. We can continue this line to arrive at:
$$Z = (\sum_{s_1=0}^{\infty}e^{-\beta\hbar s_1\omega_1})*(\sum_{s_2=0}^{\infty}e^{-\beta\hbar s_2\omega_2})*... = \prod_{i=1}Z_{\omega_i}$$
Where $Z_{\omega}$ is the partition function $\sum_{s=0}^{\infty} e^{-\beta\hbar s \omega}$ for a system consisting of photons only of frequency $\omega$.
However, the partition function of thermal blackbody radiation ends up being very different. For some reason, it is given by:
$$Z = \sum_i Z_{\omega_i} = \int D(\omega) d\omega Z_{\omega}$$
Where $Z_{\omega}$ is the same as defined above. This seems to be a contradiction. So I suppose my question is: Which of these solutions is correct?