For the partition sum, you have so sum $e^{-E}$ ($T=1$) over all possible eigenstates of the system where $E$ is the energy of the corresponding state.
Two bosons can be in the 10 states
$|kl\rangle$, with $1\leq k \leq l \leq 4$ where we accounted for the degeneracy by introducing an additional state with $E_4 =2E$. The corresponding partition sum reads (we assume the particles to be noninteracting)
$$ Z_B = \sum_{k\leq l} e^{-E_k- E_l} = 1+ e^{-E} + 3 e^{-2E} +2 e^{-3 E} +3 e^{-4E}.$$
Similarly, for fermions we have 6 states
$|kl\rangle$, with $1\leq k < l \leq 4$ with the partition sum
$$ Z_F = e^{-E} + 2 e^{-2E} +2 e^{-3 E} + e^{-4E}.$$
So the difference of the partition functions of a pair of two bosons and that of a pair of two fermions is ;-)
$$ Z_B - Z_F = 1 + e^{-2E} +2 e^{-4E}.$$
You ask
how can you formally go from the summation over $s$ to the double sum over $s_1$ and $s_2$?
As we'll see in a moment, passing to the double sum relies on the mathematical fact (about tensor products of Hilbert spaces) that if $s_1$ labels a basis of states for system $1$, and if $s_2$ labels a basis of states for system $2$, then the set of all pairs $(s_1, s_2)$ labels a basis of states for the composite system.
Let's assume that the system at hand is a quantum system described by a state space (Hilbert space) $\mathcal H$. Let the (state) labels $s$ index an orthonormal basis of $\mathcal H$ consisting of eigenvectors $|s\rangle$ of the total system's Hamiltonian $H$;
\begin{align}
H|s\rangle = E(s)|s\rangle,
\end{align}
then the partition function is obtained by summing over these states;
\begin{align}
Z = \sum_s e^{-\beta E(s)}.
\end{align}
Notice that the state labels $s$ don't have to be numbers. They could, for example, be pairs of numbers or triples of numbers, whatever is most convenient to label the possible states for the system at hand.
Now, suppose that the system consists of a pair of subsystems $1$ and $2$, then the Hilbert space of the combined system can be written as a tensor product of the Hilbert spaces for the individual subsystems;
\begin{align}
\mathcal H = \mathcal H_1\otimes \mathcal H_2.
\end{align}
Let $H_1$ denote the hamiltonian for subsystem $1$, and let $H_2$ denote the hamiltonian for subsystem $2$. Let $\{|1, s_1\rangle\}$ be an orthonormal basis for $\mathcal H_1$ consisting of eigenvectors of $H_1$, and let $\{|2, s_2\rangle\}$ be an orthonormal basis for $\mathcal H_2$ consisting of eigenvectors of $H_2$, then we have the following basic fact;
The tensor product states
\begin{align}
|s_1, s_2\rangle := |1,s_1\rangle\otimes|2,s_2\rangle
\end{align}
yield an orthonormal basis for the full state space $\mathcal H = \mathcal H_1\otimes \mathcal H_2$.
In particular, the set of states one needs to sum over in the partition function can be enumerated by pairs $s=(s_1, s_2)$. Moreover, if the systems are non-interacting, then the Hamiltonian of the full system is essentially the sum of the Hamiltonians of the individual subsystems (with appropriate identity operator factors);
\begin{align}
H = H_1\otimes I_2 + I_1\otimes H_2,
\end{align}
so that if $E_1(s_1)$ and $E_2(s_2)$ denote the energy eigenvalues of $H_1$ and $H_2$ corresponding to the states $|1,s_1\rangle$ and $|2, s_2\rangle$, then the energy $E(s_1, s_2)$ of a tensor product basis state is their sum;
\begin{align}
H|s_1,s_2\rangle = E(s_1, s_2)|s_1, s_2\rangle = \big(E_1(s_1) + E_2(s_2)\big)|s_1,s_2\rangle,
\end{align}
and the partition function can be written as a sum over the tensor product basis states;
\begin{align}
Z = \sum_{(s_1,s_2)} e^{-\beta E(s_1, s_2)} = \sum_{s_1}\sum_{s_2} e^{-\beta\big(E_1(s_2) + E_2(s_2)\big)} = \sum_{s_1} e^{-\beta E_1(s_1)} \sum_{s_2} e^{-\beta E(s_2)} = Z_1Z_2
\end{align}
where in the second equality we have used the fact that a single sum over all possible pairs $(s_1, s_2)$ is equivalent to iterated sums over the possible values of $s_1$ and $s_2$.
Best Answer
I'll show you how the book derives that solution.
First, we start with the definition for energy fluctuation: $$ \langle(\Delta E)^2\rangle =\frac{\partial^2}{\partial \beta^2}\ln(Z) $$ where $Z$ is the Canonical Partition Function. The heat capacity, $C$, can be defined as: $$ C=\frac{1}{k_BT^2}\langle (\Delta E)^2\rangle $$
Next, let's calculate $Z$ for your question: \begin{align} Z&=\sum_i g_i e^{-\beta \epsilon_i}\\ &=1\cdot e^{-\beta \epsilon}+2\cdot e^{-2\beta \epsilon}+1\cdot e^{-3\beta \epsilon} \end{align} Where the degeneracies are the factors in front of each term. The little "trick" now is to factorise this, because later we plan on taking the logarithm and so we want to write this expression as a product of factors. \begin{align} Z&=e^{-\beta \epsilon}(1+2\cdot e^{-\beta \epsilon}+e^{-2\beta \epsilon})\\ &=e^{-\beta \epsilon}(e^{-\beta \epsilon}+1)(e^{-\beta \epsilon}+1)\\ &=e^{-\beta \epsilon}(e^{-\beta \epsilon}+1)^2 \end{align}
Okay, now that we have the canonical partition function written as a product (rather than sum) we can take the natural logarithm. \begin{align} \ln(Z)&=\ln\left(e^{-\beta \epsilon}(e^{-\beta \epsilon}+1)^2\right)\\ &=\ln\left(e^{-\beta \epsilon}\right)+\ln\left((e^{-\beta \epsilon}+1)^2\right)\\ &=-\beta\epsilon+2\ln\left(e^{-\beta \epsilon}+1\right)\\ \end{align}
Next step is to calculate the energy fluctuation by taking the second derivative of this function: \begin{align} \langle(\Delta E)^2\rangle&=\frac{\partial^2}{\partial \beta^2}\ln(Z)\\ &=\frac{\partial^2}{\partial \beta^2}\left(-\beta\epsilon+2\ln\left(e^{-\beta \epsilon}+1\right)\right)\\ &=\frac{\partial}{\partial \beta}\left(-\epsilon +\frac{2(-\epsilon e^{-\beta\epsilon})}{e^{-\beta\epsilon}+1}\right)\\ &=\frac{\partial}{\partial \beta}\left(-\epsilon +\frac{-2\epsilon}{1+e^{\beta\epsilon}}\right)\\ &=-(-1)2\epsilon(\epsilon e^{\beta\epsilon})\left(1+e^{\beta\epsilon}\right)^{-2}\\ &=\frac{2\epsilon^2e^{\beta\epsilon}}{\left(1+e^{\beta\epsilon}\right)^2} \end{align} The final step is to calculate the heat capacity using the definition above: \begin{align} C&=\frac{1}{k_BT^2} \langle(\Delta E)^2\rangle\\ &=\frac{1}{k_BT^2}\frac{2\epsilon^2e^{\beta\epsilon}}{\left(1+e^{\beta\epsilon}\right)^2}\\ &=\beta^2k_B\frac{2\epsilon^2e^{\beta\epsilon}}{\left(1+e^{\beta\epsilon}\right)^2}\\ &=k_B\frac{2(\beta\epsilon)^2e^{\beta\epsilon}}{\left(1+e^{\beta\epsilon}\right)^2} \end{align} Hence, with $x=\beta\epsilon$ this can be written as: \begin{align} C=k_B\frac{2x^2 e^{x}}{\left(1+e^{x}\right)^2} \end{align} If any steps don't make sense let me know and I'll explain them!