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