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