This is related to the Coupon collector's problem as noted in the comments. Building off of [this][1] post, the probability of observing $k$ unique letters in $m$ random uniform samples from an alphabet of size $n$ is: $$\bigg\{\!{m\!\atop{k}}\bigg\}\binom{n}{k}\frac{k!}{n^m}=\bigg\{\!{m\!\atop{k}}\bigg\}\frac{n!}{n^m(n-k)!}$$ Where $\big\{\!{m\!\atop{k}}\big\}$ is the Stirling number of the second kind. For large $m$, $\ln\!\big(\big\{\!{m\!\atop{k}}\big\}\big)$ can be [approximated][2]. Here is an R function that returns the probability of every $k$: library(copula) coupons1 <- function(n, m) { l <- min(m, n) k <- 1:l if (m < 200) { logS <- log(Stirling2.all(m)[k]) } else { # estimate the log Stirling numbers v <- m/k G <- 1/v vexpv <- v/exp(v) for (i in 1:5) G <- G - (G - (vexpG <- vexpv*exp(G)))/(1 - vexpG) # Newton's method if (l == m) G[l] <- 1 logS <- (log(v - 1) - log(v*(1 - G)))/2 + (m - k)*(log(v - 1) - log(v - G)) + m*log(k) - k*log(m) + k*(1 - G) + lgamma(m + 1) - lgamma(k + 1) - lgamma(m - k + 1) } exp(logS + lgamma(n + 1) - lgamma(n - k + 1) - m*log(n)) } The probability of $k=1,2...43000$ for $n=m=43000$: system.time(k1 <- coupons1(43e3, 43e3)) #> user system elapsed #> 0.02 0.00 0.01 plot(26900:27500, k1[26900:27500], xlab = "k", ylab = "p(k)", col = "blue") [![enter image description here][3]][3] Comparing that result to a brute-force approach: Rcpp::cppFunction(" NumericVector coupons2(const int& n, const int& m) { int maxk; int n1 = n - 1; if (n > m) { maxk = m; } else { maxk = n; } NumericVector k (maxk); k(0) = 1; for (int i = 1; i < maxk; i++) { for (int j = i; j > 0; j--) { k(j) = (k(j - 1)*(n - j) + k(j)*(j + 1))/n; } k(0) = k(0)/n; } for (int i = maxk; i < m; i++) { for (int j = n1; j > 0; j--) { k(j) = (k(j - 1)*(n - j) + k(j)*(j + 1))/n; } k(0) = k(0)/n; } return k; } ") system.time(k2 <- coupons2(43e3, 43e3)) #> user system elapsed #> 12.29 0.00 12.31 The relative error from using the Stirling number approximation is small for a large $m$. max(abs(k1[26900:27500] - k2[26900:27500])/k2[26900:27500]) #> [1] 8.288009e-07 points(26900:27500, k2[26900:27500], col = "orange", pch = 20) legend("topright", legend = c("k1", "k2"), col = c("blue", "orange"), pch = c(1, 20)) [![enter image description here][4]][4] [1]: https://math.stackexchange.com/a/693254/697491 [2]: https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind#Asymptotic_approximation [3]: https://i.sstatic.net/gZYAZ.png [4]: https://i.sstatic.net/1wGdh.png