pb <- function(n) 1-exp(lfactorial(365)-lfactorial(365-n)-n*log(365))
pb(c(20,22,25,30))