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