See https://twitter.com/fermatslibrary/status/1388491536640487428?s=20

rNx.R

rNx <- function(k = 1E3, x = 1) {
  out <- numeric()
  for (i in 1:k) {
    nx <- 1
    u <- runif(1)
    while(u <= x) {
      u <- u + runif(1)
      nx <- nx + 1
    }
    out <- append(out, nx)
  }
  return(out)
}

N_1 <- rNx() # random sample of size 1000 from the distribution of N_1
mean(N_1) # estimate of 'e' by Monte Carlo
exp(1) # value of 'e'