You can use generating functions. Let me switch notations to $a = n_1$ and $b = n_2$, noticing that empty boxes make no difference. The generating function for the total number of coins is
$$
\left(\frac{1+x}{2}\right)^a \left(\frac{1+x^2}{2}\right)^b,
$$
that is, the coefficient of $x^m$ is the probability that the boxes contain a total of $m$ coins. The trick is now to substitute cubic roots of $3$ for $x$, using the following observation, where $\omega = e^{2\pi i/3}$:
$$
\frac{1^m + \omega^m + (\omega^2)^m}{3} = \begin{cases} 1 & \text{if } 3 \mid m, \\ 0 & \text{otherwise}. \end{cases}
$$
Using this, we find that the probability that you're after is
$$
\frac{1}{3} + \frac{1}{3} \left(\frac{1+\omega}{2}\right)^a \left(\frac{1+\omega^2}{2}\right)^b + \frac{1}{3} \left(\frac{1+\omega^2}{2}\right)^a \left(\frac{1+\omega}{2}\right)^b.
$$
This formula is not necessarily so helpful. Instead, we can use dynamic programming, which can be succinctly expressed using matrices. Suppose that we open the boxes one by one (each one with probability $1/2$), and keep track of the vector of probabilities of the various residues modulo $3$, a column vectors whose entries correspond to the residues $0,1,2$. An opened box with one coin affects the vector according to the matrix (multiplying the vector on the left)
$$
\begin{pmatrix}
0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0
\end{pmatrix}.
$$
If instead we open the box with probability $1/2$, we get the average of this matrix with the identity matrix:
$$
\begin{pmatrix}
\tfrac{1}{2} & 0 & \tfrac{1}{2} \\ \tfrac{1}{2} & \tfrac{1}{2} & 0 \\ 0 & \tfrac{1}{2} & \tfrac{1}{2}
\end{pmatrix}.
$$
We can analyze boxes with two coins similarly, and over all the vector of probabilities is
$$
\begin{pmatrix}
\tfrac{1}{2} & 0 & \tfrac{1}{2} \\ \tfrac{1}{2} & \tfrac{1}{2} & 0 \\ 0 & \tfrac{1}{2} & \tfrac{1}{2}
\end{pmatrix}^a
\begin{pmatrix}
\tfrac{1}{2} & \tfrac{1}{2} & 0 \\ 0 & \tfrac{1}{2} & \tfrac{1}{2} \\ \tfrac{1}{2} & 0 & \tfrac{1}{2}
\end{pmatrix}^b
\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}
$$
You can get your probability by taking the first entry. The matrix expression above can be evaluated in many ways – by multiplying the vector matrix-by-matrix, using repeated squaring to compute the matrix powers, or by diagonalizing the matrices to obtain an even more rapid algorithm (which is identical to the formula obtained above using generating functions).
Finally, let me spell out the dynamic programming algorithm, which just implements the matrix-by-matrix approach hinted to above:
- Initially, $(p_0,p_1,p_2) \gets (1,0,0)$.
- Repeat $n_1$ times: $(p_0,p_1,p_2) \gets (\frac{p_0+p_2}{2},\frac{p_0+p_1}{2},\frac{p_1+p_2}{2})$.
- Repeat $n_2$ times: $(p_0,p_1,p_2) \gets (\frac{p_0+p_1}{2},\frac{p_1+p_2}{2},\frac{p_0+p_2}{2})$.
