Generating unique binomial random variates Let the set $\{X\}$ be a collection of $n$ iid random samples of $X_i \sim \operatorname{Binomial}(m,p)$. Define $\{Y\}$ to be the set of unique values of $\{X\}$. Is there an efficient method for generating a random sample of $\{Y\}$ directly, without having to construct $\{X\}$ and finding the uniques? 
Example: For $m=64$,  $p=1/2$, and $n=10$ we could have
$X=\{28, 34, 32, 35, 29, 32, 36, 34, 30, 30\}$
$Y=\{28, 34, 32, 35, 29, 36, 30\}$
For small $n$ finding the uniques of $\{X\}$ is trivial, but if $n=10^9$ (say) we have $|X|=10^9$ but $E[|Y|]=47.05$ as
$$E[|Y|] = \sum_{i=0}^m 1- {\left(1 - p^m_i \right)}^n$$
where $p^m_i$ is the Binomial probability mass function
 A: I think you have the right idea: take a single sample from the multinomial distribution defined by the vector of $m+1$ Binomial$(m,p)$ probabilities and the integer $n.$
In general, you can obtain a sample from a multinomial distribution of size $n$ with $k$ probabilities $p_1,p_2, \ldots, p_k$ recursively : draw a Binomial$(n,p_1)$ value $Y_1$ to determine the count in the first box and then repeat with the probability vector $(p_2,p_3,\ldots,p_k)$ (renormalized to sum to unity) and the new size $n = n - Y_1.$
Eliminating the tail recursion gives the following reasonably efficient algorithm to draw one such multinomial variate:
s <- p[1] + p[2] + ... + p[k]
For i from 1 to k:
    p[i] <- p[i] / s
    x[i] <- random binomial draw of size n and probability p[i]
    n <- n - x[i]
    s <- s * (1 - p[i])
Return x

Its inputs are (1) the value of $n$ and (2) a possibly unnormalized probability vector $p$ of length $k.$  The running variables are $n,$ the updated size, and $s,$ the value needed to normalize the remaining probabilities.  By updating $s$ we have to adjust each component of $p$ just once, making this a $O(k)$ algorithm.
Note that with $n \gt 52$ or so, some of the binomial probabilities $p_i$ will be so small that $1-p_i = 1$ in double precision floating point arithmetic. I don't think this causes much harm, though, because the relative error it creates is tiny.  But once $m$ exceeds $1024$ (approximately), double-precision arithmetic will treat the tiniest probabilities (which get as low as $2^{-m}$) as true zeros.  Thus, this algorithm requires more precision when $n 2^{-m}$ is a large enough chance to affect the number of unique values.
Here is the result of $5000$ iterations of this algorithm:

The mean is $47.053 \pm 0.012.$

This is the R code I used to generate the figure.
rmult <- function(n, prob) {
  k <- length(prob)
  x <- rep(as.integer(0), k)
  s <- sum(prob)
  for (i in 1:k) {
    prob[i] <- prob[i] / s
    #
    # `rbinom` fails when n * prob[i] is too large.  Replace it with a Normal
    # approximation in such cases.
    #
    if (n * prob[i] > 1e8) {
      x[i] <- max(0, min(n, round(rnorm(1, n * prob[i], sqrt(n * prob[i] * (1-prob[i]))))))
    } else {
      x[i] <- rbinom(1, n, prob[i])
    }
    if (!is.na(x[i])) n <- n - x[i]
    s <- s * (1 - prob[i])
  }
  x
}

n <- 64
p <- dbinom(0:n, n, 1/2)

set.seed(17)
x <- replicate(5e3, sum(rmult(1e9, prob=p) > 0))
plot(table(x)/length(x), xlab="# Unique", ylab="Frequency",
     main=paste0(length(x), " Simulated Experiments"))

(round(c(Mean=mean(x), SE=sd(x) / sqrt(length(x))), 3))

