If the interaction is negligible, the partition function becomes a product of one-particle partition functions. 

Stam´s previous answer is only valid for classical systems, it is not true for quantum systems.

You didn't mention if your particles are bosons or fermions, the result will be different in each case. In the grand-canonical ensemble the partition function factorises as a product of "partition functions" associated to the occupation number of each state. Check out any stat mech book on the chapter on ideal quantum gases. For the problem you propose the grand canonical partition function is: 

 $(\sum_{n0} \exp( -\beta( E_0-\mu) n0) ) (\sum_{n1} \exp( -\beta( E_1-\mu) n1) )(\sum_{n2} \exp( -\beta(E_2-\mu) n2))$,

where mu is the chemical potential, E_0=0, E_1=E_2=W and the occupation numbers n0, n1, n2 run from 0 to infinity for bosons, or from 0 to 1 for fermions. In the canonical ensemble (take out the chemical potential mu from the equation above) the sums cannot be factorized because you have the restriction n1+n2+n0=N. For fermions, if N>=4, the partition function will be zero since the sum of the occupation numbers is at most 3 (you cannot fit 4 or more particles into only three states without violating Pauli principle). If N

That's why it's more useful, for technical convenience, in the general case, to work in the path integral formalism and introduce the appropriate chemical potentials, and boundary conditions, one doesn't need to worry explicitly about the combinatorics. For it's the path integral that computes the partition function. Incidentally, the entanglement that is the result of Bose-Einstein or Fermi-Dirac statistics is an interaction, albeit a long range one-and that's relevant only if the particles are indistinguishable.  So one does need to be careful and not confuse history with physics.

Hi Piyush,

In case you still need some help with the boson case, here are a few tips to get the calculations done. As I suggested earlier, compute first the grand canonical partition function, which can be factorized. You should find Qg(z)=1/[(1-z)(1-z exp(-beta W))^2], where z=exp(beta mu). To obtain from there the canonical partition function Qc(N), remember that Qg(z)=sum_N Qc(N) z^N : the canonical partition functions for each N are the coefficients of the Taylor expansion of the grand canonical partition function at z=0. So, Qc(N) can be computed as the residue of Qg(z)/z^{N+1} at z=0. This residue is not simple to compute because this pole is of order N+1, but since Qg(z) decays faster than 1/z^2 when z->infinity, you can prove that the sum of all residues of Qg(z)/z^{N+1} is zero. Then you only need to compute the other residues at z=1 (order 1) and z=exp(beta W)  (of order 2) to obtain your answer.

Hope this helps!  This is a nice exercise, thanks for bringing it up.  You have given me some material for a homework or exam for my students ;-)

By the way, the technique above can also be used for the fermion case, if one does not want to deal with the combinatorics explicitly (although they are trivial in this case). The grand canonical partition function for the fermion case is Qg(z)=(1+z)(1+z exp(-beta W))^2. Simply expand this order 3 polynomial and read off the canonical partition functions for N=0,1,2,3 from the coefficients of this polynomial.

That´s the beauty of the generating functions, isn't it?

$Q=\frac{1}{N!}  \left( 1+2\exp(-\beta W)\right)^N