Efficient random generation for generalized birthday paradox problem?

The probability of getting $k$ unique values from $[0, n)$ when choosing $m$ times is given by:

$$P(V = k) = \binom{n}{k}\displaystyle\sum_{i=0}^k (-1)^i \binom{k}{i} \left(\frac{k-i}{n}\right)^m$$

where $V$ is a random variable giving the number of unique outcomes and $\binom{\cdot}{\cdot}$ is the binomial coefficient.

Sampling from such distribution can be easily achieved by sampling from discrete uniform distribution and then counting the number of unique values. Is there any better way from random generation from this distribution? I call the direct simulation inefficient because it needs generating each time $m$ samples from the discrete uniform distribution.

For a given value of $(n,m)$ [to be read as (days,draws)], one can compute exactly $$\mathbb{P}(V = k) = \binom{n}{k}\displaystyle\sum_{i=0}^k (-1)^i \binom{k}{i} \left(\frac{k-i}{n}\right)^m$$for all values of $k$ at a cost of $\text{O}(n²)$ and then simulate exactly from this distribution by generating $U(0,1)$ variates and comparing them with the cumulated sums of the $\mathbb{P}(V = k)$'s.
As noted by W. Huber, this computation may and should run into instability issues. Using the recurrence formula provided by Random $(1\le j\le m\wedge n)$ \begin{align} \mathbb{P}_{n,m=1}(V=1)=p_{n,1}(1)&=1\\ \mathbb{P}_{n,m+1}(V=j)=p_{n,m+1}(j)&=\frac{j}{n}p_{n,m}(j)+\frac{n-j+1}{n}p_{n,m}(j-1)\\ \end{align} leads to a manageable (if presumably approximate) probability distribution.
It would be nice to exploit the summation representation through a latent variable connected with a Binomial $\mathfrak{B}(k,\rho_k)$ distribution, but the negative terms lead to a possibly negative expression I do not know how to handle.
• Beware: this computation is unstable due to the alternating series. For realistic values of $n,m,k$ (such as those involved in the original Birthday Problem), all accuracy is lost if you use only double-precision floating point arithmetic. For instance, with $n=365, m=50, k=47$ the size of many of the terms exceeds $10^{19}$, but the result must lie between $0$ and $1$. Therefore (in this particular case) at least $20$ decimal places must be retained throughout the calculation in order to get anywhere near an approximate answer. – whuber Oct 31 '16 at 19:15