How to generate geometric random numbers?

[Robert Kern]

Paul Rubin wrote:
> Paul Rubin <http://phr.cx@NOSPAM.invalid> writes:
>>    n = log(c, 1-p) - 1
> 
> I meant    n = log(c/p, 1-p) - 1
> sorry.

import random
from math import ceil, log

def geometric(p):
   """ Geometric distribution per Devroye, Luc. _Non-Uniform Random Variate
   Generation_, 1986, p 500. http://cg.scs.carleton.ca/~luc/rnbookindex.html
   """

   # p should be in (0.0, 1.0].
   if p <= 0.0 or p > 1.0:
     raise ValueError("p must be in the interval (0.0, 1.0]")
   elif p == 1.0:
     # If p is exactly 1.0, then the only possible generated value is 1.
     # Recognizing this case early means that we can avoid a log(0.0) later.
     # The exact floating point comparison should be fine. log(eps) works just
     # dandy.
     return 1

   # random() returns a number in [0, 1). The log() function does not
   # like 0.
   U = 1.0 - random.random()

   # Find the corresponding geometric variate by inverting the uniform variate.
   G = int(ceil(log(U) / log(1.0 - p)))
   return G

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
  that is made terrible by our own mad attempt to interpret it as though it had
  an underlying truth."
   -- Umberto Eco