p = b*u                if p <= 1.0:                    x = p ** (1.0/alpha)                else:                    x = -_log((b-p)/alpha)                u1 = random()                if p > 1.0:                    if u1 <= x ** (alpha - 1.0):                        break                elif u1 <= _exp(-x):                    break            return x * beta ## -------------------- Gauss (faster alternative) --------------------     def gauss(self, mu, sigma):        """Gaussian distribution.         mu is the mean, and sigma is the standard deviation.  This is        slightly faster than the normalvariate() function.         Not thread-safe without a lock around calls.         """         # When x and y are two variables from [0, 1), uniformly        # distributed, then        #        #    cos(2*pi*x)*sqrt(-2*log(1-y))        #    sin(2*pi*x)*sqrt(-2*log(1-y))        #        # are two *independent* variables with normal distribution        # (mu = 0, sigma = 1).        # (Lambert Meertens)        # (corrected version; bug discovered by Mike Miller, fixed by LM)         # Multithreading note: When two threads call this function        # simultaneously, it is possible that they will receive the        # same return value.  The window is very small though.  To        # avoid this, you have to use a lock around all calls.  (I        # didn't want to slow this down in the serial case by using a        # lock here.)         random = self.random        z = self.gauss_next        self.gauss_next = None        if z is None:            x2pi = random() * TWOPI            g2rad = _sqrt(-2.0 * _log(1.0 - random()))            z = _cos(x2pi) * g2rad            self.gauss_next = _sin(x2pi) * g2rad         return mu + z*sigma ## -------------------- beta --------------------## See## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470## for Ivan Frohne's insightful analysis of why the original implementation:####    def betavariate(self, alpha, beta):##        # Discrete Event Simulation in C, pp 87-88.####        y = self.expovariate(alpha)##        z = self.expovariate(1.0/beta)##        return z/(y+z)#### was dead wrong, and how it probably got that way.     def betavariate(self, alpha, beta):        """Beta distribution.         Conditions on the parameters are alpha > 0 and beta > 0.        Returned values range between 0 and 1.         """         # This version due to Janne Sinkkonen, and matches all the std        # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").        y = self.gammavariate(alpha, 1.)        if y == 0:            return 0.0        else:            return y / (y + self.gammavariate(beta, 1.)) ## -------------------- Pareto --------------------     def paretovariate(self, alpha):        """Pareto distribution.  alpha is the shape parameter."""        # Jain, pg. 495         u = 1.0 - self.random()        return 1.0 / pow(u, 1.0/alpha) ## -------------------- Weibull --------------------     def weibullvariate(self, alpha, beta):        """Weibull distribution.         alpha is the scale parameter and beta is the shape parameter.         """        # Jain, pg. 499; bug fix courtesy Bill Arms         u = 1.0 - self.random()        return alpha * pow(-_log(u), 1.0/beta) ## -------------------- Wichmann-Hill ------------------- class WichmannHill(Random):     VERSION = 1     # used by getstate/setstate     def seed(self, a=None):        """Initialize internal state from hashable object.         None or no argument seeds from current time or from an operating        system specific randomness source if available.         If a is not None or an int or long, hash(a) is used instead.         If a is an int or long, a is used directly.  Distinct values between        0 and 27814431486575L inclusive are guaranteed to yield distinct        internal states (this guarantee is specific to the default        Wichmann-Hill generator).        """         if a is None:            try:                a = long(_hexlify(_urandom(16)), 16)            except NotImplementedError:                import time                a = long(time.time() * 256) # use fractional seconds         if not isinstance(a, (int, long)):            a = hash(a)         a, x = divmod(a, 30268)        a, y = divmod(a, 30306)        a, z = divmod(a, 30322)        self._seed = int(x)+1, int(y)+1, int(z)+1         self.gauss_next = None     def random(self):        """Get the next random number in the range [0.0, 1.0)."""         # Wichman-Hill random number generator.        #        # Wichmann, B. A. & Hill, I. D. (1982)        # Algorithm AS 183:        # An efficient and portable pseudo-random number generator        # Applied Statistics 31 (1982) 188-190        #        # see also:        #        Correction to Algorithm AS 183        #        Applied Statistics 33 (1984) 123        #        #        McLeod, A. I. (1985)        #        A remark on Algorithm AS 183        #        Applied Statistics 34 (1985),198-200         # This part is thread-unsafe:        # BEGIN CRITICAL SECTION        x, y, z = self._seed        x = (171 * x) % 30269        y = (172 * y) % 30307        z = (170 * z) % 30323        self._seed = x, y, z        # END CRITICAL SECTION         # Note:  on a platform using IEEE-754 double arithmetic, this can        # never return 0.0 (asserted by Tim; proof too long for a comment).        return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0     def getstate(self):        """Return internal state; can be passed to setstate() later."""        return self.VERSION, self._seed, self.gauss_next     def setstate(self, state):        """Restore internal state from object returned by getstate()."""        version = state[0]        if version == 1:            version, self._seed, self.gauss_next = state        else:            raise ValueError("state with version %s passed to "                             "Random.setstate() of version %s" %                             (version, self.VERSION))     def jumpahead(self, n):        """Act as if n calls to random() were made, but quickly.         n is an int, greater than or equal to 0.         Example use:  If you have 2 threads and know that each will        consume no more than a million random numbers, create two Random        objects r1 and r2, then do            r2.setstate(r1.getstate())            r2.jumpahead(1000000)        Then r1 and r2 will use guaranteed-disjoint segments of the full        period.        """         if not n >= 0:            raise ValueError("n must be >= 0")        x, y, z = self._seed        x = int(x * pow(171, n, 30269)) % 30269        y = int(y * pow(172, n, 30307)) % 30307        z = int(z * pow(170, n, 30323)) % 30323        self._seed = x, y, z     def __whseed(self, x=0, y=0, z=0):        """Set the Wichmann-Hill seed from (x, y, z).         These must be integers in the range [0, 256).        """         if not type(x) == type(y) == type(z) == int:            raise TypeError('seeds must be integers')        if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):            raise ValueError('seeds must be in range(0, 256)')        if 0 == x == y == z:            # Initialize from current time            import time            t = long(time.time() * 256)            t = int((t&0xffffff) ^ (t>>24))            t, x = divmod(t, 256)            t, y = divmod(t, 256)            t, z = divmod(t, 256)        # Zero is a poor seed, so substitute 1        self._seed = (x or 1, y or 1, z or 1)         self.gauss_next = None     def whseed(self, a=None):        """Seed from hashable object's hash code.         None or no argument seeds from current time.  It is not guaranteed        that objects with distinct hash codes lead to distinct internal        states.         This is obsolete, provided for compatibility with the seed routine        used prior to Python 2.1.  Use the .seed() method instead.        """         if a is None:            self.__whseed()            return        a = hash(a)        a, x = divmod(a, 256)        a, y = divmod(a, 256)        a, z = divmod(a, 256)        x = (x + a) % 256 or 1        y = (y + a) % 256 or 1        z = (z + a) % 256 or 1        self.__whseed(x, y, z) ## --------------- Operating System Random Source  ------------------ class SystemRandom(Random):    """Alternate random number generator using sources provided    by the operating system (such as /dev/urandom on Unix or    CryptGenRandom on Windows).      Not available on all systems (see os.urandom() for details).    """     def random(self):        """Get the next random number in the range [0.0, 1.0)."""        return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF     def getrandbits(self, k):        """getrandbits(k) -> x.  Generates a long int with k random bits."""        if k <= 0:            raise ValueError('number of bits must be greater than zero')        if k != int(k):            raise TypeError('number of bits should be an integer')        bytes = (k + 7) // 8                    # bits / 8 and rounded up        x = long(_hexlify(_urandom(bytes)), 16)        return x >> (bytes * 8 - k)             # trim excess bits     def _stub(self, *args, **kwds):        "Stub method.  Not used for a system random number generator."        return None    seed = jumpahead = _stub     def _notimplemented(self, *args, **kwds):        "Method should not be called for a system random number generator."        raise NotImplementedError('System entropy source does not have state.')    getstate = setstate = _notimplemented ## -------------------- test program -------------------- def _test_generator(n, func, args):    import time    print n, 'times', func.__name__    total = 0.0    sqsum = 0.0    smallest = 1e10    largest = -1e10    t0 = time.time()    for i in range(n):        x = func(*args)        total += x        sqsum = sqsum + x*x        smallest = min(x, smallest)        largest = max(x, largest)    t1 = time.time()    print round(t1-t0, 3), 'sec,',    avg = total/n    stddev = _sqrt(sqsum/n - avg*avg)    print 'avg %g, stddev %g, min %g, max %g' % \              (avg, stddev, smallest, largest)  def _test(N=2000):    _test_generator(N, random, ())    _test_generator(N, normalvariate, (0.0, 1.0))    _test_generator(N, lognormvariate, (0.0, 1.0))    _test_generator(N, vonmisesvariate, (0.0, 1.0))    _test_generator(N, gammavariate, (0.01, 1.0))    _test_generator(N, gammavariate, (0.1, 1.0))    _test_generator(N, gammavariate, (0.1, 2.0))    _test_generator(N, gammavariate, (0.5, 1.0))    _test_generator(N, gammavariate, (0.9, 1.0))    _test_generator(N, gammavariate, (1.0, 1.0))    _test_generator(N, gammavariate, (2.0, 1.0))    _test_generator(N, gammavariate, (20.0, 1.0))    _test_generator(N, gammavariate, (200.0, 1.0))    _test_generator(N, gauss, (0.0, 1.0))    _test_generator(N, betavariate, (3.0, 3.0)) # Create one instance, seeded from current time, and export its methods# as module-level functions.  The functions share state across all uses#(both in the user's code and in the Python libraries), but that's fine# for most programs and is easier for the casual user than making them# instantiate their own Random() instance. _inst = Random()seed = _inst.seedrandom = _inst.randomuniform = _inst.uniformrandint = _inst.randintchoice = _inst.choicerandrange = _inst.randrangesample = _inst.sampleshuffle = _inst.shufflenormalvariate = _inst.normalvariatelognormvariate = _inst.lognormvariateexpovariate = _inst.expovariatevonmisesvariate = _inst.vonmisesvariategammavariate = _inst.gammavariategauss = _inst.gaussbetavariate = _inst.betavariateparetovariate = _inst.paretovariateweibullvariate = _inst.weibullvariategetstate = _inst.getstatesetstate = _inst.setstatejumpahead = _inst.jumpaheadgetrandbits = _inst.getrandbits if __name__ == '__main__':    _test()

Lines 501–862 of 862. Previous lines