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