Here is a solution to the problem. Suppose that there are $N>0$ boxes to start with. Look at the boxes just before the last round, and color one black and the other white. If we move back in time, we can partition the boxes into those giving rise to the black box, which we also color black, and those giving rise to the white box, which we color white. For $n\ge 0$ and $k\in\{0,\dots,n\}$, let $p_{n,k}$ be the probability that, $n$ rounds before the last round, we have $k+1$ white boxes and $n+1-k$ black boxes. We prove by induction on $n$ that $p_{n,k}=1/(n+1)$, regardless of $k$. If $n=0$, this is obvious. Otherwise, observe that if we have $k+1$ white boxes $n$ rounds before the last, we must have had either $k$ white boxes and $n+1-k$ black boxes $n-1$ rounds before the last (case A), or $k+1$ white boxes and $n-k$ black boxes (case B). In case A, we must have chosen to combine two boxes into one of the $k$ white boxes in the preceding round, which will happen with probability $k/(n+1)$. In case B, we must have chosen to combine two boxes into one of the $n-k$ black boxes, which will happen with probability $(n-k)/(n+1)$. Therefore, $ p_{n,k}=\frac{k}{n+1} p_{n-1,k-1} + \frac{n-k}{n+1} p_{n-1,k}, $ and, by the induction hypothesis, this equals $1/(n+1)$, as desired.
Now, if we condition on a given partition of the boxes into white and black, the behavior of the white boxes is the same as if we had started the problem with only the white boxes, and similarly for the black boxes. Also, in the final round:
- If we combine two boxes with 1 coin, we get a box with 2 coins.
- If we combine a box with 1 coin with one with 2 or more coins, we get a box with 1 coin.
- If we combine two boxes with 2 coins, we get a box with 3 or more coins.
- If we combine a box with 2 coins with a box with 3 or more coins, we get a box with 2 coins.
- If we combine two boxes with 3 or more coins, we get a box with 3 or more coins.
Therefore, if we let $q_N$, $r_N$, and $s_N$ be the probabilities that, starting with $N$ coins, we end up with a box with 1 coin, 2 coins, or 3 or more coins, we have the recurrence relations
$q_N=\delta_{N1}+\frac{1}{N-1} \sum_{i+j=N, i\ge 1, j\ge 1} q_i (r_j+s_j)+(r_i+s_i)q_j, \qquad (1)$ $r_N=\frac{1}{N-1} \sum_{i+j=N, i\ge 1, j\ge 1} q_i q_j + r_i s_j+s_i r_j,\qquad (2)$ $s_N=\frac{1}{N-1} \sum_{i+j=N, i\ge 1, j\ge 1} r_i r_j + s_i s_j,\qquad (3)$
where the sums are over the numbers $i$ and $j$ of white and black boxes we start with at the beginning.
It immediately follows by induction that, if $N$ is odd, $q_N=1$ and $r_N=s_N=0$, and, if $N$ is even, $q_N=0$ and $r_N+s_N=1$. This is because, as already mentioned, the parity of the number of boxes which contain 1 coin is preserved by each round. Therefore, from (2), for all $n>0$, $r_{2n}=\frac{n}{2n-1}+\frac{1}{2n-1} \sum_{i+j=n, i\ge 1, j\ge 1} r_{2i} (1-r_{2j})+(1-r_{2i}) r_{2j}.\qquad (4)$ It follows from this that, if $\lim_{n\to\infty} r_{2n}$ exists and equals $L$, say, then $L=\frac{1}{2}+L(1-L)$, so $L^2=\frac12$ and, since $L\ge 0$, $L=1/{\sqrt{2}}$. It only remains to prove that $\lim_{n\to\infty} r_{2n}$ exists.
If we multiply each side of (4) by $(2n-1)x^{2n}$ and sum over $n$, we get $ x \frac{\partial \phi}{\partial x}-\phi=\frac{x^2}{(1-x^2)^2}+ 2\phi(\frac{x^2}{1-x^2}-\phi), $ where $\phi=\sum_{n\ge 1} r_{2n} x^{2n}.$ This differential equation can be simplified by setting $ x=\frac{1-t}{1+t},\qquad \phi = \theta (t-t^{-1}), $ and then becomes, after dividing through by $(t-t^{-1})^2$, $ \frac{t}{2} \frac{\partial \theta}{\partial t}=\frac{1}{16}-2\theta^2. $ This is separable and can be integrated to $ \theta=-\frac{1}{4\sqrt{2}} \frac{1-Ct^{\sqrt{2}}}{1+Ct^{\sqrt{2}}} $ for a constant $C$, so $ \phi=\frac{x}{\sqrt{2}(1-x^2)} \frac{ (1+x)^{\sqrt{2}}-C(1-x)^{\sqrt{2}} }{ (1+x)^{\sqrt{2}}+C(1-x)^{\sqrt{2}}.} $ Since the coefficients of 1 and $x$ in $\phi$ must vanish, we must have $C=1$. Let $z_0:=(1-\exp \pi i/\sqrt{2})/(1+\exp \pi i/\sqrt{2})\approx -2.018i $. $\phi$ is then analytic in the complex plane with $(-\infty,-1]\cup[1,\infty)\cup\{\pm z_0\}$ removed, and we have $ \phi=\frac{1}{2\sqrt{2}(1-x)}(1+O(x-1)) \qquad \hbox{near $x=1$,} $ $ \phi=\frac{1}{2\sqrt{2}(1+x)}(1+O(x+1)) \qquad \hbox{near $x=-1$.} $ Then, by singularity analysis, $ r_{2n}=[x^{2n}]\phi=\frac{1}{\sqrt{2}}+O(\frac{1}{n}). $ This completes the solution.