The "appropriate exponential distribution" for the number of bins $K$ is not exponential by any means. It actually is the discrete distribution of parameter $n$ given by
$$\Pr(K=k)\ \propto\ \frac{k^n}{k!}\tag{1.4}$$
for $k=1, 2, \ldots.$
In the original paper, A. J. Stam argued as follows. When a partition $\pi$ consists of $|\pi|\ge 1$ nonempty subsets of $\{1,2,\ldots,n\},$ the conditional probability of realizing $\pi$ under this sample scheme after randomly choosing $K=k$ bins is
$$\Pr(\pi\mid K=k) = k^{(|\pi|)}k^{-n} = \frac{k(k-1)\cdots(k-|\pi|+1)}{k^n}.\tag{1.5}$$
(The equation numbering is Stam's, but I have slightly changed the notation of the variables.)
Let's justify this formula. The denominator counts the number of sequences of bin choices: each random toss can land in one of the $k$ bins and all $n$ tosses are independent. Thus each sequence of tosses has chance $k^{-n}.$ However, two sequences of tosses determine the same partition when the bins for one can be re-ordered to correspond to the bins of the other. If we order the bins so that the first partition occupies bins $1,2,\ldots, |\pi|,$ then the possible equivalent re-orderings are determined by putting bin $1$ at one of the $k$ positions, bin $2$ at one of the $k-1$ remaining positions, and so on, until we have repositioned bin $|\pi|.$ Thus, any partition $\pi$ of $|\pi|$ pieces created by tossing $n$ numbers into $k$ bins shows up multiple times as counted by the numerator of $(1.5).$
In conjunction with $(1.4),$ the unconditional probability may be found by summing over all possible values of $K.$ Evidently values $K=0, 1, \ldots, |\pi|-1$ cannot give $\pi,$ justifying the change in the summation index in the penultimate step:
$$\begin{aligned}
\Pr(\pi) &= \sum_{k=0}^\infty\Pr(\pi\mid K=k)\Pr(K=k) \\
&\propto\sum_{k=0}^\infty\frac{k(k-1)\cdots(k-|\pi|+1)}{k^n}\,\frac{k^n}{k!}\\
&= \sum_{k=|\pi|}^\infty\frac{1}{(k-|\pi|)!} = \sum_{i=0}^\infty \frac{1}{i!} = e.
\end{aligned}$$
Because this probability does not depend on $\pi,$ the distribution is uniform, QED.
As a consequence of this argument we recover Dobinski's formula for the number of partitions $T_n$ of $\{1,2,\ldots, n\},$ because it must be the normalizing constant of this distribution, whence (employing $(1.4)$)
$$T_n = \Pr(\pi)^{-1} = \left(\frac{e}{\sum_{k=1}^\infty \Pr(K=k)}\right)^{-1} = \frac{1}{e}\sum_{k=1} \frac{k^n}{k!}.\tag{1.1}$$
For ways to sample from discrete distributions like $(1.4),$ see our thread on this topic. The alias method is efficient.
