Efficient random generation for generalized birthday paradox problem?

Question

In 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.

Answer 1

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.