Lib/random.py

   1 """Random variable generators.
   2
   3     integers
   4     --------
   5            uniform within range
   6
   7     sequences
   8     ---------
   9            pick random element
  10            pick random sample
  11            generate random permutation
  12
  13     distributions on the real line:
  14     ------------------------------
  15            uniform
  16            triangular
  17            normal (Gaussian)
  18            lognormal
  19            negative exponential
  20            gamma
  21            beta
  22            pareto
  23            Weibull
  24
  25     distributions on the circle (angles 0 to 2pi)
  26     ---------------------------------------------
  27            circular uniform
  28            von Mises
  29
  30 General notes on the underlying Mersenne Twister core generator:
  31
  32 * The period is 2**19937-1.
  33 * It is one of the most extensively tested generators in existence.
  34 * Without a direct way to compute N steps forward, the semantics of
  35   jumpahead(n) are weakened to simply jump to another distant state and rely
  36   on the large period to avoid overlapping sequences.
  37 * The random() method is implemented in C, executes in a single Python step,
  38   and is, therefore, threadsafe.
  39
  40 """
  41
  42 from __future__ import division
  43 from warnings import warn as _warn
  44 from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
  45 from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
  46 from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
  47 from os import urandom as _urandom
  48 from binascii import hexlify as _hexlify
  49
  50 __all__ = ["Random","seed","random","uniform","randint","choice","sample",
  51            "randrange","shuffle","normalvariate","lognormvariate",
  52            "expovariate","vonmisesvariate","gammavariate","triangular",
  53            "gauss","betavariate","paretovariate","weibullvariate",
  54            "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
  55            "SystemRandom"]
  56
  57 NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
  58 TWOPI = 2.0*_pi
  59 LOG4 = _log(4.0)
  60 SG_MAGICCONST = 1.0 + _log(4.5)
  61 BPF = 53        # Number of bits in a float
  62 RECIP_BPF = 2**-BPF
  63
  64
  65 # Translated by Guido van Rossum from C source provided by
  66 # Adrian Baddeley.  Adapted by Raymond Hettinger for use with
  67 # the Mersenne Twister  and os.urandom() core generators.
  68
  69 import _random
  70
  71 class Random(_random.Random):
  72     """Random number generator base class used by bound module functions.
  73
  74     Used to instantiate instances of Random to get generators that don't
  75     share state.  Especially useful for multi-threaded programs, creating
  76     a different instance of Random for each thread, and using the jumpahead()
  77     method to ensure that the generated sequences seen by each thread don't
  78     overlap.
  79
  80     Class Random can also be subclassed if you want to use a different basic
  81     generator of your own devising: in that case, override the following
  82     methods: random(), seed(), getstate(), setstate() and jumpahead().
  83     Optionally, implement a getrandbits() method so that randrange() can cover
  84     arbitrarily large ranges.
  85
  86     """
  87
  88     VERSION = 3     # used by getstate/setstate
  89
  90     def __init__(self, x=None):
  91         """Initialize an instance.
  92
  93         Optional argument x controls seeding, as for Random.seed().
  94         """
  95
  96         self.seed(x)
  97         self.gauss_next = None
  98
  99     def seed(self, a=None):
 100         """Initialize internal state from hashable object.
 101
 102         None or no argument seeds from current time or from an operating
 103         system specific randomness source if available.
 104
 105         If a is not None or an int or long, hash(a) is used instead.
 106         """
 107
 108         if a is None:
 109             try:
 110                 a = long(_hexlify(_urandom(16)), 16)
 111             except NotImplementedError:
 112                 import time
 113                 a = long(time.time() * 256) # use fractional seconds
 114
 115         super(Random, self).seed(a)
 116         self.gauss_next = None
 117
 118     def getstate(self):
 119         """Return internal state; can be passed to setstate() later."""
 120         return self.VERSION, super(Random, self).getstate(), self.gauss_next
 121
 122     def setstate(self, state):
 123         """Restore internal state from object returned by getstate()."""
 124         version = state[0]
 125         if version == 3:
 126             version, internalstate, self.gauss_next = state
 127             super(Random, self).setstate(internalstate)
 128         elif version == 2:
 129             version, internalstate, self.gauss_next = state
 130             # In version 2, the state was saved as signed ints, which causes
 131             #   inconsistencies between 32/64-bit systems. The state is
 132             #   really unsigned 32-bit ints, so we convert negative ints from
 133             #   version 2 to positive longs for version 3.
 134             try:
 135                 internalstate = tuple( long(x) % (2**32) for x in internalstate )
 136             except ValueError, e:
 137                 raise TypeError, e
 138             super(Random, self).setstate(internalstate)
 139         else:
 140             raise ValueError("state with version %s passed to "
 141                              "Random.setstate() of version %s" %
 142                              (version, self.VERSION))
 143
 144 ## ---- Methods below this point do not need to be overridden when
 145 ## ---- subclassing for the purpose of using a different core generator.
 146
 147 ## -------------------- pickle support  -------------------
 148
 149     def __getstate__(self): # for pickle
 150         return self.getstate()
 151
 152     def __setstate__(self, state):  # for pickle
 153         self.setstate(state)
 154
 155     def __reduce__(self):
 156         return self.__class__, (), self.getstate()
 157
 158 ## -------------------- integer methods  -------------------
 159
 160     def randrange(self, start, stop=None, step=1, int=int, default=None,
 161                   maxwidth=1L<<BPF):
 162         """Choose a random item from range(start, stop[, step]).
 163
 164         This fixes the problem with randint() which includes the
 165         endpoint; in Python this is usually not what you want.
 166         Do not supply the 'int', 'default', and 'maxwidth' arguments.
 167         """
 168
 169         # This code is a bit messy to make it fast for the
 170         # common case while still doing adequate error checking.
 171         istart = int(start)
 172         if istart != start:
 173             raise ValueError, "non-integer arg 1 for randrange()"
 174         if stop is default:
 175             if istart > 0:
 176                 if istart >= maxwidth:
 177                     return self._randbelow(istart)
 178                 return int(self.random() * istart)
 179             raise ValueError, "empty range for randrange()"
 180
 181         # stop argument supplied.
 182         istop = int(stop)
 183         if istop != stop:
 184             raise ValueError, "non-integer stop for randrange()"
 185         width = istop - istart
 186         if step == 1 and width > 0:
 187             # Note that
 188             #     int(istart + self.random()*width)
 189             # instead would be incorrect.  For example, consider istart
 190             # = -2 and istop = 0.  Then the guts would be in
 191             # -2.0 to 0.0 exclusive on both ends (ignoring that random()
 192             # might return 0.0), and because int() truncates toward 0, the
 193             # final result would be -1 or 0 (instead of -2 or -1).
 194             #     istart + int(self.random()*width)
 195             # would also be incorrect, for a subtler reason:  the RHS
 196             # can return a long, and then randrange() would also return
 197             # a long, but we're supposed to return an int (for backward
 198             # compatibility).
 199
 200             if width >= maxwidth:
 201                 return int(istart + self._randbelow(width))
 202             return int(istart + int(self.random()*width))
 203         if step == 1:
 204             raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)
 205
 206         # Non-unit step argument supplied.
 207         istep = int(step)
 208         if istep != step:
 209             raise ValueError, "non-integer step for randrange()"
 210         if istep > 0:
 211             n = (width + istep - 1) // istep
 212         elif istep < 0:
 213             n = (width + istep + 1) // istep
 214         else:
 215             raise ValueError, "zero step for randrange()"
 216
 217         if n <= 0:
 218             raise ValueError, "empty range for randrange()"
 219
 220         if n >= maxwidth:
 221             return istart + istep*self._randbelow(n)
 222         return istart + istep*int(self.random() * n)
 223
 224     def randint(self, a, b):
 225         """Return random integer in range [a, b], including both end points.
 226         """
 227
 228         return self.randrange(a, b+1)
 229
 230     def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF,
 231                    _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
 232         """Return a random int in the range [0,n)
 233
 234         Handles the case where n has more bits than returned
 235         by a single call to the underlying generator.
 236         """
 237
 238         try:
 239             getrandbits = self.getrandbits
 240         except AttributeError:
 241             pass
 242         else:
 243             # Only call self.getrandbits if the original random() builtin method
 244             # has not been overridden or if a new getrandbits() was supplied.
 245             # This assures that the two methods correspond.
 246             if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
 247                 k = int(1.00001 + _log(n-1, 2.0))   # 2**k > n-1 > 2**(k-2)
 248                 r = getrandbits(k)
 249                 while r >= n:
 250                     r = getrandbits(k)
 251                 return r
 252         if n >= _maxwidth:
 253             _warn("Underlying random() generator does not supply \n"
 254                 "enough bits to choose from a population range this large")
 255         return int(self.random() * n)
 256
 257 ## -------------------- sequence methods  -------------------
 258
 259     def choice(self, seq):
 260         """Choose a random element from a non-empty sequence."""
 261         return seq[int(self.random() * len(seq))]  # raises IndexError if seq is empty
 262
 263     def shuffle(self, x, random=None, int=int):
 264         """x, random=random.random -> shuffle list x in place; return None.
 265
 266         Optional arg random is a 0-argument function returning a random
 267         float in [0.0, 1.0); by default, the standard random.random.
 268         """
 269
 270         if random is None:
 271             random = self.random
 272         for i in reversed(xrange(1, len(x))):
 273             # pick an element in x[:i+1] with which to exchange x[i]
 274             j = int(random() * (i+1))
 275             x[i], x[j] = x[j], x[i]
 276
 277     def sample(self, population, k):
 278         """Chooses k unique random elements from a population sequence.
 279
 280         Returns a new list containing elements from the population while
 281         leaving the original population unchanged.  The resulting list is
 282         in selection order so that all sub-slices will also be valid random
 283         samples.  This allows raffle winners (the sample) to be partitioned
 284         into grand prize and second place winners (the subslices).
 285
 286         Members of the population need not be hashable or unique.  If the
 287         population contains repeats, then each occurrence is a possible
 288         selection in the sample.
 289
 290         To choose a sample in a range of integers, use xrange as an argument.
 291         This is especially fast and space efficient for sampling from a
 292         large population:   sample(xrange(10000000), 60)
 293         """
 294
 295         # XXX Although the documentation says `population` is "a sequence",
 296         # XXX attempts are made to cater to any iterable with a __len__
 297         # XXX method.  This has had mixed success.  Examples from both
 298         # XXX sides:  sets work fine, and should become officially supported;
 299         # XXX dicts are much harder, and have failed in various subtle
 300         # XXX ways across attempts.  Support for mapping types should probably
 301         # XXX be dropped (and users should pass mapping.keys() or .values()
 302         # XXX explicitly).
 303
 304         # Sampling without replacement entails tracking either potential
 305         # selections (the pool) in a list or previous selections in a set.
 306
 307         # When the number of selections is small compared to the
 308         # population, then tracking selections is efficient, requiring
 309         # only a small set and an occasional reselection.  For
 310         # a larger number of selections, the pool tracking method is
 311         # preferred since the list takes less space than the
 312         # set and it doesn't suffer from frequent reselections.
 313
 314         n = len(population)
 315         if not 0 <= k <= n:
 316             raise ValueError, "sample larger than population"
 317         random = self.random
 318         _int = int
 319         result = [None] * k
 320         setsize = 21        # size of a small set minus size of an empty list
 321         if k > 5:
 322             setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
 323         if n <= setsize or hasattr(population, "keys"):
 324             # An n-length list is smaller than a k-length set, or this is a
 325             # mapping type so the other algorithm wouldn't work.
 326             pool = list(population)
 327             for i in xrange(k):         # invariant:  non-selected at [0,n-i)
 328                 j = _int(random() * (n-i))
 329                 result[i] = pool[j]
 330                 pool[j] = pool[n-i-1]   # move non-selected item into vacancy
 331         else:
 332             try:
 333                 selected = set()
 334                 selected_add = selected.add
 335                 for i in xrange(k):
 336                     j = _int(random() * n)
 337                     while j in selected:
 338                         j = _int(random() * n)
 339                     selected_add(j)
 340                     result[i] = population[j]
 341             except (TypeError, KeyError):   # handle (at least) sets
 342                 if isinstance(population, list):
 343                     raise
 344                 return self.sample(tuple(population), k)
 345         return result
 346
 347 ## -------------------- real-valued distributions  -------------------
 348
 349 ## -------------------- uniform distribution -------------------
 350
 351     def uniform(self, a, b):
 352         "Get a random number in the range [a, b) or [a, b] depending on rounding."
 353         return a + (b-a) * self.random()
 354
 355 ## -------------------- triangular --------------------
 356
 357     def triangular(self, low=0.0, high=1.0, mode=None):
 358         """Triangular distribution.
 359
 360         Continuous distribution bounded by given lower and upper limits,
 361         and having a given mode value in-between.
 362
 363         http://en.wikipedia.org/wiki/Triangular_distribution
 364
 365         """
 366         u = self.random()
 367         c = 0.5 if mode is None else (mode - low) / (high - low)
 368         if u > c:
 369             u = 1.0 - u
 370             c = 1.0 - c
 371             low, high = high, low
 372         return low + (high - low) * (u * c) ** 0.5
 373
 374 ## -------------------- normal distribution --------------------
 375
 376     def normalvariate(self, mu, sigma):
 377         """Normal distribution.
 378
 379         mu is the mean, and sigma is the standard deviation.
 380
 381         """
 382         # mu = mean, sigma = standard deviation
 383
 384         # Uses Kinderman and Monahan method. Reference: Kinderman,
 385         # A.J. and Monahan, J.F., "Computer generation of random
 386         # variables using the ratio of uniform deviates", ACM Trans
 387         # Math Software, 3, (1977), pp257-260.
 388
 389         random = self.random
 390         while 1:
 391             u1 = random()
 392             u2 = 1.0 - random()
 393             z = NV_MAGICCONST*(u1-0.5)/u2
 394             zz = z*z/4.0
 395             if zz <= -_log(u2):
 396                 break
 397         return mu + z*sigma
 398
 399 ## -------------------- lognormal distribution --------------------
 400
 401     def lognormvariate(self, mu, sigma):
 402         """Log normal distribution.
 403
 404         If you take the natural logarithm of this distribution, you'll get a
 405         normal distribution with mean mu and standard deviation sigma.
 406         mu can have any value, and sigma must be greater than zero.
 407
 408         """
 409         return _exp(self.normalvariate(mu, sigma))
 410
 411 ## -------------------- exponential distribution --------------------
 412
 413     def expovariate(self, lambd):
 414         """Exponential distribution.
 415
 416         lambd is 1.0 divided by the desired mean.  It should be
 417         nonzero.  (The parameter would be called "lambda", but that is
 418         a reserved word in Python.)  Returned values range from 0 to
 419         positive infinity if lambd is positive, and from negative
 420         infinity to 0 if lambd is negative.
 421
 422         """
 423         # lambd: rate lambd = 1/mean
 424         # ('lambda' is a Python reserved word)
 425
 426         random = self.random
 427         u = random()
 428         while u <= 1e-7:
 429             u = random()
 430         return -_log(u)/lambd
 431
 432 ## -------------------- von Mises distribution --------------------
 433
 434     def vonmisesvariate(self, mu, kappa):
 435         """Circular data distribution.
 436
 437         mu is the mean angle, expressed in radians between 0 and 2*pi, and
 438         kappa is the concentration parameter, which must be greater than or
 439         equal to zero.  If kappa is equal to zero, this distribution reduces
 440         to a uniform random angle over the range 0 to 2*pi.
 441
 442         """
 443         # mu:    mean angle (in radians between 0 and 2*pi)
 444         # kappa: concentration parameter kappa (>= 0)
 445         # if kappa = 0 generate uniform random angle
 446
 447         # Based upon an algorithm published in: Fisher, N.I.,
 448         # "Statistical Analysis of Circular Data", Cambridge
 449         # University Press, 1993.
 450
 451         # Thanks to Magnus Kessler for a correction to the
 452         # implementation of step 4.
 453
 454         random = self.random
 455         if kappa <= 1e-6:
 456             return TWOPI * random()
 457
 458         a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
 459         b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
 460         r = (1.0 + b * b)/(2.0 * b)
 461
 462         while 1:
 463             u1 = random()
 464
 465             z = _cos(_pi * u1)
 466             f = (1.0 + r * z)/(r + z)
 467             c = kappa * (r - f)
 468
 469             u2 = random()
 470
 471             if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
 472                 break
 473
 474         u3 = random()
 475         if u3 > 0.5:
 476             theta = (mu % TWOPI) + _acos(f)
 477         else:
 478             theta = (mu % TWOPI) - _acos(f)
 479
 480         return theta
 481
 482 ## -------------------- gamma distribution --------------------
 483
 484     def gammavariate(self, alpha, beta):
 485         """Gamma distribution.  Not the gamma function!
 486
 487         Conditions on the parameters are alpha > 0 and beta > 0.
 488
 489         """
 490
 491         # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
 492
 493         # Warning: a few older sources define the gamma distribution in terms
 494         # of alpha > -1.0
 495         if alpha <= 0.0 or beta <= 0.0:
 496             raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
 497
 498         random = self.random
 499         if alpha > 1.0:
 500
 501             # Uses R.C.H. Cheng, "The generation of Gamma
 502             # variables with non-integral shape parameters",
 503             # Applied Statistics, (1977), 26, No. 1, p71-74
 504
 505             ainv = _sqrt(2.0 * alpha - 1.0)
 506             bbb = alpha - LOG4
 507             ccc = alpha + ainv
 508
 509             while 1:
 510                 u1 = random()
 511                 if not 1e-7 < u1 < .9999999:
 512                     continue
 513                 u2 = 1.0 - random()
 514                 v = _log(u1/(1.0-u1))/ainv
 515                 x = alpha*_exp(v)
 516                 z = u1*u1*u2
 517                 r = bbb+ccc*v-x
 518                 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
 519                     return x * beta
 520
 521         elif alpha == 1.0:
 522             # expovariate(1)
 523             u = random()
 524             while u <= 1e-7:
 525                 u = random()
 526             return -_log(u) * beta
 527
 528         else:   # alpha is between 0 and 1 (exclusive)
 529
 530             # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
 531
 532             while 1:
 533                 u = random()
 534                 b = (_e + alpha)/_e
 535                 p = b*u
 536                 if p <= 1.0:
 537                     x = p ** (1.0/alpha)
 538                 else:
 539                     x = -_log((b-p)/alpha)
 540                 u1 = random()
 541                 if p > 1.0:
 542                     if u1 <= x ** (alpha - 1.0):
 543                         break
 544                 elif u1 <= _exp(-x):
 545                     break
 546             return x * beta
 547
 548 ## -------------------- Gauss (faster alternative) --------------------
 549
 550     def gauss(self, mu, sigma):
 551         """Gaussian distribution.
 552
 553         mu is the mean, and sigma is the standard deviation.  This is
 554         slightly faster than the normalvariate() function.
 555
 556         Not thread-safe without a lock around calls.
 557
 558         """
 559
 560         # When x and y are two variables from [0, 1), uniformly
 561         # distributed, then
 562         #
 563         #    cos(2*pi*x)*sqrt(-2*log(1-y))
 564         #    sin(2*pi*x)*sqrt(-2*log(1-y))
 565         #
 566         # are two *independent* variables with normal distribution
 567         # (mu = 0, sigma = 1).
 568         # (Lambert Meertens)
 569         # (corrected version; bug discovered by Mike Miller, fixed by LM)
 570
 571         # Multithreading note: When two threads call this function
 572         # simultaneously, it is possible that they will receive the
 573         # same return value.  The window is very small though.  To
 574         # avoid this, you have to use a lock around all calls.  (I
 575         # didn't want to slow this down in the serial case by using a
 576         # lock here.)
 577
 578         random = self.random
 579         z = self.gauss_next
 580         self.gauss_next = None
 581         if z is None:
 582             x2pi = random() * TWOPI
 583             g2rad = _sqrt(-2.0 * _log(1.0 - random()))
 584             z = _cos(x2pi) * g2rad
 585             self.gauss_next = _sin(x2pi) * g2rad
 586
 587         return mu + z*sigma
 588
 589 ## -------------------- beta --------------------
 590 ## See
 591 ## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
 592 ## for Ivan Frohne's insightful analysis of why the original implementation:
 593 ##
 594 ##    def betavariate(self, alpha, beta):
 595 ##        # Discrete Event Simulation in C, pp 87-88.
 596 ##
 597 ##        y = self.expovariate(alpha)
 598 ##        z = self.expovariate(1.0/beta)
 599 ##        return z/(y+z)
 600 ##
 601 ## was dead wrong, and how it probably got that way.
 602
 603     def betavariate(self, alpha, beta):
 604         """Beta distribution.
 605
 606         Conditions on the parameters are alpha > 0 and beta > 0.
 607         Returned values range between 0 and 1.
 608
 609         """
 610
 611         # This version due to Janne Sinkkonen, and matches all the std
 612         # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
 613         y = self.gammavariate(alpha, 1.)
 614         if y == 0:
 615             return 0.0
 616         else:
 617             return y / (y + self.gammavariate(beta, 1.))
 618
 619 ## -------------------- Pareto --------------------
 620
 621     def paretovariate(self, alpha):
 622         """Pareto distribution.  alpha is the shape parameter."""
 623         # Jain, pg. 495
 624
 625         u = 1.0 - self.random()
 626         return 1.0 / pow(u, 1.0/alpha)
 627
 628 ## -------------------- Weibull --------------------
 629
 630     def weibullvariate(self, alpha, beta):
 631         """Weibull distribution.
 632
 633         alpha is the scale parameter and beta is the shape parameter.
 634
 635         """
 636         # Jain, pg. 499; bug fix courtesy Bill Arms
 637
 638         u = 1.0 - self.random()
 639         return alpha * pow(-_log(u), 1.0/beta)
 640
 641 ## -------------------- Wichmann-Hill -------------------
 642
 643 class WichmannHill(Random):
 644
 645     VERSION = 1     # used by getstate/setstate
 646
 647     def seed(self, a=None):
 648         """Initialize internal state from hashable object.
 649
 650         None or no argument seeds from current time or from an operating
 651         system specific randomness source if available.
 652
 653         If a is not None or an int or long, hash(a) is used instead.
 654
 655         If a is an int or long, a is used directly.  Distinct values between
 656         0 and 27814431486575L inclusive are guaranteed to yield distinct
 657         internal states (this guarantee is specific to the default
 658         Wichmann-Hill generator).
 659         """
 660
 661         if a is None:
 662             try:
 663                 a = long(_hexlify(_urandom(16)), 16)
 664             except NotImplementedError:
 665                 import time
 666                 a = long(time.time() * 256) # use fractional seconds
 667
 668         if not isinstance(a, (int, long)):
 669             a = hash(a)
 670
 671         a, x = divmod(a, 30268)
 672         a, y = divmod(a, 30306)
 673         a, z = divmod(a, 30322)
 674         self._seed = int(x)+1, int(y)+1, int(z)+1
 675
 676         self.gauss_next = None
 677
 678     def random(self):
 679         """Get the next random number in the range [0.0, 1.0)."""
 680
 681         # Wichman-Hill random number generator.
 682         #
 683         # Wichmann, B. A. & Hill, I. D. (1982)
 684         # Algorithm AS 183:
 685         # An efficient and portable pseudo-random number generator
 686         # Applied Statistics 31 (1982) 188-190
 687         #
 688         # see also:
 689         #        Correction to Algorithm AS 183
 690         #        Applied Statistics 33 (1984) 123
 691         #
 692         #        McLeod, A. I. (1985)
 693         #        A remark on Algorithm AS 183
 694         #        Applied Statistics 34 (1985),198-200
 695
 696         # This part is thread-unsafe:
 697         # BEGIN CRITICAL SECTION
 698         x, y, z = self._seed
 699         x = (171 * x) % 30269
 700         y = (172 * y) % 30307
 701         z = (170 * z) % 30323
 702         self._seed = x, y, z
 703         # END CRITICAL SECTION
 704
 705         # Note:  on a platform using IEEE-754 double arithmetic, this can
 706         # never return 0.0 (asserted by Tim; proof too long for a comment).
 707         return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
 708
 709     def getstate(self):
 710         """Return internal state; can be passed to setstate() later."""
 711         return self.VERSION, self._seed, self.gauss_next
 712
 713     def setstate(self, state):
 714         """Restore internal state from object returned by getstate()."""
 715         version = state[0]
 716         if version == 1:
 717             version, self._seed, self.gauss_next = state
 718         else:
 719             raise ValueError("state with version %s passed to "
 720                              "Random.setstate() of version %s" %
 721                              (version, self.VERSION))
 722
 723     def jumpahead(self, n):
 724         """Act as if n calls to random() were made, but quickly.
 725
 726         n is an int, greater than or equal to 0.
 727
 728         Example use:  If you have 2 threads and know that each will
 729         consume no more than a million random numbers, create two Random
 730         objects r1 and r2, then do
 731             r2.setstate(r1.getstate())
 732             r2.jumpahead(1000000)
 733         Then r1 and r2 will use guaranteed-disjoint segments of the full
 734         period.
 735         """
 736
 737         if not n >= 0:
 738             raise ValueError("n must be >= 0")
 739         x, y, z = self._seed
 740         x = int(x * pow(171, n, 30269)) % 30269
 741         y = int(y * pow(172, n, 30307)) % 30307
 742         z = int(z * pow(170, n, 30323)) % 30323
 743         self._seed = x, y, z
 744
 745     def __whseed(self, x=0, y=0, z=0):
 746         """Set the Wichmann-Hill seed from (x, y, z).
 747
 748         These must be integers in the range [0, 256).
 749         """
 750
 751         if not type(x) == type(y) == type(z) == int:
 752             raise TypeError('seeds must be integers')
 753         if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
 754             raise ValueError('seeds must be in range(0, 256)')
 755         if 0 == x == y == z:
 756             # Initialize from current time
 757             import time
 758             t = long(time.time() * 256)
 759             t = int((t&0xffffff) ^ (t>>24))
 760             t, x = divmod(t, 256)
 761             t, y = divmod(t, 256)
 762             t, z = divmod(t, 256)
 763         # Zero is a poor seed, so substitute 1
 764         self._seed = (x or 1, y or 1, z or 1)
 765
 766         self.gauss_next = None
 767
 768     def whseed(self, a=None):
 769         """Seed from hashable object's hash code.
 770
 771         None or no argument seeds from current time.  It is not guaranteed
 772         that objects with distinct hash codes lead to distinct internal
 773         states.
 774
 775         This is obsolete, provided for compatibility with the seed routine
 776         used prior to Python 2.1.  Use the .seed() method instead.
 777         """
 778
 779         if a is None:
 780             self.__whseed()
 781             return
 782         a = hash(a)
 783         a, x = divmod(a, 256)
 784         a, y = divmod(a, 256)
 785         a, z = divmod(a, 256)
 786         x = (x + a) % 256 or 1
 787         y = (y + a) % 256 or 1
 788         z = (z + a) % 256 or 1
 789         self.__whseed(x, y, z)
 790
 791 ## --------------- Operating System Random Source  ------------------
 792
 793 class SystemRandom(Random):
 794     """Alternate random number generator using sources provided
 795     by the operating system (such as /dev/urandom on Unix or
 796     CryptGenRandom on Windows).
 797
 798      Not available on all systems (see os.urandom() for details).
 799     """
 800
 801     def random(self):
 802         """Get the next random number in the range [0.0, 1.0)."""
 803         return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
 804
 805     def getrandbits(self, k):
 806         """getrandbits(k) -> x.  Generates a long int with k random bits."""
 807         if k <= 0:
 808             raise ValueError('number of bits must be greater than zero')
 809         if k != int(k):
 810             raise TypeError('number of bits should be an integer')
 811         bytes = (k + 7) // 8                    # bits / 8 and rounded up
 812         x = long(_hexlify(_urandom(bytes)), 16)
 813         return x >> (bytes * 8 - k)             # trim excess bits
 814
 815     def _stub(self, *args, **kwds):
 816         "Stub method.  Not used for a system random number generator."
 817         return None
 818     seed = jumpahead = _stub
 819
 820     def _notimplemented(self, *args, **kwds):
 821         "Method should not be called for a system random number generator."
 822         raise NotImplementedError('System entropy source does not have state.')
 823     getstate = setstate = _notimplemented
 824
 825 ## -------------------- test program --------------------
 826
 827 def _test_generator(n, func, args):
 828     import time
 829     print n, 'times', func.__name__
 830     total = 0.0
 831     sqsum = 0.0
 832     smallest = 1e10
 833     largest = -1e10
 834     t0 = time.time()
 835     for i in range(n):
 836         x = func(*args)
 837         total += x
 838         sqsum = sqsum + x*x
 839         smallest = min(x, smallest)
 840         largest = max(x, largest)
 841     t1 = time.time()
 842     print round(t1-t0, 3), 'sec,',
 843     avg = total/n
 844     stddev = _sqrt(sqsum/n - avg*avg)
 845     print 'avg %g, stddev %g, min %g, max %g' % \
 846               (avg, stddev, smallest, largest)
 847
 848
 849 def _test(N=2000):
 850     _test_generator(N, random, ())
 851     _test_generator(N, normalvariate, (0.0, 1.0))
 852     _test_generator(N, lognormvariate, (0.0, 1.0))
 853     _test_generator(N, vonmisesvariate, (0.0, 1.0))
 854     _test_generator(N, gammavariate, (0.01, 1.0))
 855     _test_generator(N, gammavariate, (0.1, 1.0))
 856     _test_generator(N, gammavariate, (0.1, 2.0))
 857     _test_generator(N, gammavariate, (0.5, 1.0))
 858     _test_generator(N, gammavariate, (0.9, 1.0))
 859     _test_generator(N, gammavariate, (1.0, 1.0))
 860     _test_generator(N, gammavariate, (2.0, 1.0))
 861     _test_generator(N, gammavariate, (20.0, 1.0))
 862     _test_generator(N, gammavariate, (200.0, 1.0))
 863     _test_generator(N, gauss, (0.0, 1.0))
 864     _test_generator(N, betavariate, (3.0, 3.0))
 865     _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0))
 866
 867 # Create one instance, seeded from current time, and export its methods
 868 # as module-level functions.  The functions share state across all uses
 869 #(both in the user's code and in the Python libraries), but that's fine
 870 # for most programs and is easier for the casual user than making them
 871 # instantiate their own Random() instance.
 872
 873 _inst = Random()
 874 seed = _inst.seed
 875 random = _inst.random
 876 uniform = _inst.uniform
 877 triangular = _inst.triangular
 878 randint = _inst.randint
 879 choice = _inst.choice
 880 randrange = _inst.randrange
 881 sample = _inst.sample
 882 shuffle = _inst.shuffle
 883 normalvariate = _inst.normalvariate
 884 lognormvariate = _inst.lognormvariate
 885 expovariate = _inst.expovariate
 886 vonmisesvariate = _inst.vonmisesvariate
 887 gammavariate = _inst.gammavariate
 888 gauss = _inst.gauss
 889 betavariate = _inst.betavariate
 890 paretovariate = _inst.paretovariate
 891 weibullvariate = _inst.weibullvariate
 892 getstate = _inst.getstate
 893 setstate = _inst.setstate
 894 jumpahead = _inst.jumpahead
 895 getrandbits = _inst.getrandbits
 896
 897 if __name__ == '__main__':
 898     _test()