
my.rbinom <- function(n, p) {
  # MOLTO INEFFICIENTE!!
  fxm <- max(dbinom(0:n, n, p))
  max.prove <- 100
  prova = 1

  while(prova < max.prove) {
    x <- sample(0:n)[1]
    fx <- dbinom(x, n, p)
    rf <- runif(1)*fxm
    if (rf < fx) break()    
  }
  return(x)
}

# test
test = TRUE
if(test) {
  N = 10000
  x <- rep(0, N)
  for (i in 1:N)  x[i] <- my.rbinom(10, 0.5)
  plot(table(x))
}
