Indeed, the problem was already solved. It's a generalization of Birthday paradox. The formula for probability mass function is
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.
as quoted by George V. Williams on math.stackexchange.com who refers to http://www.randomservices.org/random/urn/Birthday.html#General
set.seed(123)
m <- 100
n <- 40
f <- function(x, m, n) {
vapply(x, function(k) {
choose(n, k) * sum( ((-1)^(0:k) * choose(k, (0:k)) * ((k - (0:k))/n)^m) )
}, numeric(1))
}
out <- replicate(1e6, {
length(unique(sample.int(n, m, replace = TRUE)))
})
plot(prop.table(table(out)), xlab = "k", ylab = "P(V=k)")
lines(1:m, f(1:m, m, n), col = "red")