Lib/random.py

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