Попытка понять метод обратного преобразования для генерации пуассоновской случайной величины

Question

Алгоритм Википедии дает для генерации распределенных по Пуассону случайных величин с использованием метода обратного преобразования:

init:
     Let x ← 0, p ← e^−λ, s ← p.
     Generate uniform random number u in [0,1].
while u > s do:
     x ← x + 1.
     p ← p * λ / x.
     s ← s + p.
return x.

Я реализовал это в R:

f<-function(lambda)
{
x<-0
p<-exp(-lambda)
s<-p    
u<-runif(1)    
while (u>s )
{
   x<-x+1    
   p<-p*lambda/x    
   s<-s+p
}    
return(x)
}

но я не понимаю, как это генерирует значения этой случайной величины, и я также думаю, что код неполон, потому что вывод всегда 1.

Может кто-нибудь объяснить, пожалуйста?