pb <- function(n) 1-exp(lfactorial(365)-lfactorial(365-n)-n*log(365))
plot(1:100, pb(1:100),
  xlab="Number of students in the class", ylab="Probability of a common birthday")