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