birthday <- function(n, b) exp(lfactorial(365) + (n-b)*log(365-b) - lfactorial(365-b) - n*log(365))
plot(1:80, mapply(birthday, 80, 1:80))