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$, $\big\{\!{m\!\atop{k}}\big\}$ can be [approximated][2] ([original paper][3]): $$\left\{{n \atop k}\right\} \sim \sqrt{\frac{v-1}{v(1-G)}} \left(\frac{v-1}{v-G}\right)^{n-k} \frac{k^n}{n^k} e^{k(1-G)} \left({n \atop k}\right)$$ where $v=n/k$, and $G\in(0,1)$ is the unique solution to $G=ve^{G-v}$. 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 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) if (l == m) logS[m] <- 0 } 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][4]][4] 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][5]][5] A note on the `occupancy` R package: ------------------------------------ The `occupancy` package features the distribution for this problem, with an option for exact or approximate calculations. `docc` returns the PMF of `k`. However, it is much slower than `coupons1` and `coupons2` above, with the computation time growing approximately geometrically for $m=n$. Additionally, the approximation provided is poor compared to the one provided by `coupons1`: n <- 43e3L system.time(k3 <- occupancy::docc(1:n, n, n, approx = TRUE)) #> user system elapsed #> 1.61 2.00 3.61 max(abs(k3[26900:27500] - k2[26900:27500])/k2[26900:27500]) #> [1] 0.01156602 Time `occupancy::docc` for $m=n=4300$ with `approx = FALSE`: n <- 4300L system.time(occupancy::docc(1:n, n, n)) #> user system elapsed #> 35.58 0.14 35.75 [1]: https://math.stackexchange.com/a/693254/697491 [2]: https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind#Asymptotic_approximation [3]: https://core.ac.uk/download/pdf/301651745.pdf [4]: https://i.sstatic.net/gZYAZ.png [5]: https://i.sstatic.net/1wGdh.png