n <- 1:1500
a <- n * (1 - 1/365)^n * 365
plot(a, type = "l")
polyroot(c(365, -1, -1)) # 365 - n - n^2 = 0
## [1] 18.61151+0i -19.61151+0i
n <- 19
sample(1:365, size = n, replace = TRUE)
## [1] 30 120 107 297 287 266 307 11 118 115 183 324 246 183 364 340 344 242 359