| """Random variable generators. | 
 |  | 
 |     integers | 
 |     -------- | 
 |            uniform within range | 
 |  | 
 |     sequences | 
 |     --------- | 
 |            pick random element | 
 |            pick random sample | 
 |            generate random permutation | 
 |  | 
 |     distributions on the real line: | 
 |     ------------------------------ | 
 |            uniform | 
 |            triangular | 
 |            normal (Gaussian) | 
 |            lognormal | 
 |            negative exponential | 
 |            gamma | 
 |            beta | 
 |            pareto | 
 |            Weibull | 
 |  | 
 |     distributions on the circle (angles 0 to 2pi) | 
 |     --------------------------------------------- | 
 |            circular uniform | 
 |            von Mises | 
 |  | 
 | General notes on the underlying Mersenne Twister core generator: | 
 |  | 
 | * The period is 2**19937-1. | 
 | * It is one of the most extensively tested generators in existence. | 
 | * The random() method is implemented in C, executes in a single Python step, | 
 |   and is, therefore, threadsafe. | 
 |  | 
 | """ | 
 |  | 
 | from __future__ import division | 
 | from warnings import warn as _warn | 
 | from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType | 
 | from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil | 
 | from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin | 
 | from os import urandom as _urandom | 
 | from collections import Set as _Set, Sequence as _Sequence | 
 | from hashlib import sha512 as _sha512 | 
 |  | 
 | __all__ = ["Random","seed","random","uniform","randint","choice","sample", | 
 |            "randrange","shuffle","normalvariate","lognormvariate", | 
 |            "expovariate","vonmisesvariate","gammavariate","triangular", | 
 |            "gauss","betavariate","paretovariate","weibullvariate", | 
 |            "getstate","setstate", "getrandbits", | 
 |            "SystemRandom"] | 
 |  | 
 | NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) | 
 | TWOPI = 2.0*_pi | 
 | LOG4 = _log(4.0) | 
 | SG_MAGICCONST = 1.0 + _log(4.5) | 
 | BPF = 53        # Number of bits in a float | 
 | RECIP_BPF = 2**-BPF | 
 |  | 
 |  | 
 | # Translated by Guido van Rossum from C source provided by | 
 | # Adrian Baddeley.  Adapted by Raymond Hettinger for use with | 
 | # the Mersenne Twister  and os.urandom() core generators. | 
 |  | 
 | import _random | 
 |  | 
 | class Random(_random.Random): | 
 |     """Random number generator base class used by bound module functions. | 
 |  | 
 |     Used to instantiate instances of Random to get generators that don't | 
 |     share state. | 
 |  | 
 |     Class Random can also be subclassed if you want to use a different basic | 
 |     generator of your own devising: in that case, override the following | 
 |     methods:  random(), seed(), getstate(), and setstate(). | 
 |     Optionally, implement a getrandbits() method so that randrange() | 
 |     can cover arbitrarily large ranges. | 
 |  | 
 |     """ | 
 |  | 
 |     VERSION = 3     # used by getstate/setstate | 
 |  | 
 |     def __init__(self, x=None): | 
 |         """Initialize an instance. | 
 |  | 
 |         Optional argument x controls seeding, as for Random.seed(). | 
 |         """ | 
 |  | 
 |         self.seed(x) | 
 |         self.gauss_next = None | 
 |  | 
 |     def seed(self, a=None, version=2): | 
 |         """Initialize internal state from hashable object. | 
 |  | 
 |         None or no argument seeds from current time or from an operating | 
 |         system specific randomness source if available. | 
 |  | 
 |         For version 2 (the default), all of the bits are used if *a *is a str, | 
 |         bytes, or bytearray.  For version 1, the hash() of *a* is used instead. | 
 |  | 
 |         If *a* is an int, all bits are used. | 
 |  | 
 |         """ | 
 |  | 
 |         if a is None: | 
 |             try: | 
 |                 a = int.from_bytes(_urandom(32), 'big') | 
 |             except NotImplementedError: | 
 |                 import time | 
 |                 a = int(time.time() * 256) # use fractional seconds | 
 |  | 
 |         if version == 2: | 
 |             if isinstance(a, (str, bytes, bytearray)): | 
 |                 if isinstance(a, str): | 
 |                     a = a.encode() | 
 |                 a += _sha512(a).digest() | 
 |                 a = int.from_bytes(a, 'big') | 
 |  | 
 |         super().seed(a) | 
 |         self.gauss_next = None | 
 |  | 
 |     def getstate(self): | 
 |         """Return internal state; can be passed to setstate() later.""" | 
 |         return self.VERSION, super().getstate(), self.gauss_next | 
 |  | 
 |     def setstate(self, state): | 
 |         """Restore internal state from object returned by getstate().""" | 
 |         version = state[0] | 
 |         if version == 3: | 
 |             version, internalstate, self.gauss_next = state | 
 |             super().setstate(internalstate) | 
 |         elif version == 2: | 
 |             version, internalstate, self.gauss_next = state | 
 |             # In version 2, the state was saved as signed ints, which causes | 
 |             #   inconsistencies between 32/64-bit systems. The state is | 
 |             #   really unsigned 32-bit ints, so we convert negative ints from | 
 |             #   version 2 to positive longs for version 3. | 
 |             try: | 
 |                 internalstate = tuple(x % (2**32) for x in internalstate) | 
 |             except ValueError as e: | 
 |                 raise TypeError from e | 
 |             super().setstate(internalstate) | 
 |         else: | 
 |             raise ValueError("state with version %s passed to " | 
 |                              "Random.setstate() of version %s" % | 
 |                              (version, self.VERSION)) | 
 |  | 
 | ## ---- Methods below this point do not need to be overridden when | 
 | ## ---- subclassing for the purpose of using a different core generator. | 
 |  | 
 | ## -------------------- pickle support  ------------------- | 
 |  | 
 |     def __getstate__(self): # for pickle | 
 |         return self.getstate() | 
 |  | 
 |     def __setstate__(self, state):  # for pickle | 
 |         self.setstate(state) | 
 |  | 
 |     def __reduce__(self): | 
 |         return self.__class__, (), self.getstate() | 
 |  | 
 | ## -------------------- integer methods  ------------------- | 
 |  | 
 |     def randrange(self, start, stop=None, step=1, int=int): | 
 |         """Choose a random item from range(start, stop[, step]). | 
 |  | 
 |         This fixes the problem with randint() which includes the | 
 |         endpoint; in Python this is usually not what you want. | 
 |  | 
 |         Do not supply the 'int' argument. | 
 |         """ | 
 |  | 
 |         # This code is a bit messy to make it fast for the | 
 |         # common case while still doing adequate error checking. | 
 |         istart = int(start) | 
 |         if istart != start: | 
 |             raise ValueError("non-integer arg 1 for randrange()") | 
 |         if stop is None: | 
 |             if istart > 0: | 
 |                 return self._randbelow(istart) | 
 |             raise ValueError("empty range for randrange()") | 
 |  | 
 |         # stop argument supplied. | 
 |         istop = int(stop) | 
 |         if istop != stop: | 
 |             raise ValueError("non-integer stop for randrange()") | 
 |         width = istop - istart | 
 |         if step == 1 and width > 0: | 
 |             return istart + self._randbelow(width) | 
 |         if step == 1: | 
 |             raise ValueError("empty range for randrange() (%d,%d, %d)" % (istart, istop, width)) | 
 |  | 
 |         # Non-unit step argument supplied. | 
 |         istep = int(step) | 
 |         if istep != step: | 
 |             raise ValueError("non-integer step for randrange()") | 
 |         if istep > 0: | 
 |             n = (width + istep - 1) // istep | 
 |         elif istep < 0: | 
 |             n = (width + istep + 1) // istep | 
 |         else: | 
 |             raise ValueError("zero step for randrange()") | 
 |  | 
 |         if n <= 0: | 
 |             raise ValueError("empty range for randrange()") | 
 |  | 
 |         return istart + istep*self._randbelow(n) | 
 |  | 
 |     def randint(self, a, b): | 
 |         """Return random integer in range [a, b], including both end points. | 
 |         """ | 
 |  | 
 |         return self.randrange(a, b+1) | 
 |  | 
 |     def _randbelow(self, n, int=int, maxsize=1<<BPF, type=type, | 
 |                    Method=_MethodType, BuiltinMethod=_BuiltinMethodType): | 
 |         "Return a random int in the range [0,n).  Raises ValueError if n==0." | 
 |  | 
 |         getrandbits = self.getrandbits | 
 |         # Only call self.getrandbits if the original random() builtin method | 
 |         # has not been overridden or if a new getrandbits() was supplied. | 
 |         if type(self.random) is BuiltinMethod or type(getrandbits) is Method: | 
 |             k = n.bit_length()  # don't use (n-1) here because n can be 1 | 
 |             r = getrandbits(k)          # 0 <= r < 2**k | 
 |             while r >= n: | 
 |                 r = getrandbits(k) | 
 |             return r | 
 |         # There's an overriden random() method but no new getrandbits() method, | 
 |         # so we can only use random() from here. | 
 |         random = self.random | 
 |         if n >= maxsize: | 
 |             _warn("Underlying random() generator does not supply \n" | 
 |                 "enough bits to choose from a population range this large.\n" | 
 |                 "To remove the range limitation, add a getrandbits() method.") | 
 |             return int(random() * n) | 
 |         rem = maxsize % n | 
 |         limit = (maxsize - rem) / maxsize   # int(limit * maxsize) % n == 0 | 
 |         r = random() | 
 |         while r >= limit: | 
 |             r = random() | 
 |         return int(r*maxsize) % n | 
 |  | 
 | ## -------------------- sequence methods  ------------------- | 
 |  | 
 |     def choice(self, seq): | 
 |         """Choose a random element from a non-empty sequence.""" | 
 |         try: | 
 |             i = self._randbelow(len(seq)) | 
 |         except ValueError: | 
 |             raise IndexError('Cannot choose from an empty sequence') | 
 |         return seq[i] | 
 |  | 
 |     def shuffle(self, x, random=None, int=int): | 
 |         """x, random=random.random -> shuffle list x in place; return None. | 
 |  | 
 |         Optional arg random is a 0-argument function returning a random | 
 |         float in [0.0, 1.0); by default, the standard random.random. | 
 |         """ | 
 |  | 
 |         randbelow = self._randbelow | 
 |         for i in reversed(range(1, len(x))): | 
 |             # pick an element in x[:i+1] with which to exchange x[i] | 
 |             j = randbelow(i+1) if random is None else int(random() * (i+1)) | 
 |             x[i], x[j] = x[j], x[i] | 
 |  | 
 |     def sample(self, population, k): | 
 |         """Chooses k unique random elements from a population sequence or set. | 
 |  | 
 |         Returns a new list containing elements from the population while | 
 |         leaving the original population unchanged.  The resulting list is | 
 |         in selection order so that all sub-slices will also be valid random | 
 |         samples.  This allows raffle winners (the sample) to be partitioned | 
 |         into grand prize and second place winners (the subslices). | 
 |  | 
 |         Members of the population need not be hashable or unique.  If the | 
 |         population contains repeats, then each occurrence is a possible | 
 |         selection in the sample. | 
 |  | 
 |         To choose a sample in a range of integers, use range as an argument. | 
 |         This is especially fast and space efficient for sampling from a | 
 |         large population:   sample(range(10000000), 60) | 
 |         """ | 
 |  | 
 |         # Sampling without replacement entails tracking either potential | 
 |         # selections (the pool) in a list or previous selections in a set. | 
 |  | 
 |         # When the number of selections is small compared to the | 
 |         # population, then tracking selections is efficient, requiring | 
 |         # only a small set and an occasional reselection.  For | 
 |         # a larger number of selections, the pool tracking method is | 
 |         # preferred since the list takes less space than the | 
 |         # set and it doesn't suffer from frequent reselections. | 
 |  | 
 |         if isinstance(population, _Set): | 
 |             population = tuple(population) | 
 |         if not isinstance(population, _Sequence): | 
 |             raise TypeError("Population must be a sequence or set.  For dicts, use list(d).") | 
 |         randbelow = self._randbelow | 
 |         n = len(population) | 
 |         if not 0 <= k <= n: | 
 |             raise ValueError("Sample larger than population") | 
 |         result = [None] * k | 
 |         setsize = 21        # size of a small set minus size of an empty list | 
 |         if k > 5: | 
 |             setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets | 
 |         if n <= setsize: | 
 |             # An n-length list is smaller than a k-length set | 
 |             pool = list(population) | 
 |             for i in range(k):         # invariant:  non-selected at [0,n-i) | 
 |                 j = randbelow(n-i) | 
 |                 result[i] = pool[j] | 
 |                 pool[j] = pool[n-i-1]   # move non-selected item into vacancy | 
 |         else: | 
 |             selected = set() | 
 |             selected_add = selected.add | 
 |             for i in range(k): | 
 |                 j = randbelow(n) | 
 |                 while j in selected: | 
 |                     j = randbelow(n) | 
 |                 selected_add(j) | 
 |                 result[i] = population[j] | 
 |         return result | 
 |  | 
 | ## -------------------- real-valued distributions  ------------------- | 
 |  | 
 | ## -------------------- uniform distribution ------------------- | 
 |  | 
 |     def uniform(self, a, b): | 
 |         "Get a random number in the range [a, b) or [a, b] depending on rounding." | 
 |         return a + (b-a) * self.random() | 
 |  | 
 | ## -------------------- triangular -------------------- | 
 |  | 
 |     def triangular(self, low=0.0, high=1.0, mode=None): | 
 |         """Triangular distribution. | 
 |  | 
 |         Continuous distribution bounded by given lower and upper limits, | 
 |         and having a given mode value in-between. | 
 |  | 
 |         http://en.wikipedia.org/wiki/Triangular_distribution | 
 |  | 
 |         """ | 
 |         u = self.random() | 
 |         c = 0.5 if mode is None else (mode - low) / (high - low) | 
 |         if u > c: | 
 |             u = 1.0 - u | 
 |             c = 1.0 - c | 
 |             low, high = high, low | 
 |         return low + (high - low) * (u * c) ** 0.5 | 
 |  | 
 | ## -------------------- normal distribution -------------------- | 
 |  | 
 |     def normalvariate(self, mu, sigma): | 
 |         """Normal distribution. | 
 |  | 
 |         mu is the mean, and sigma is the standard deviation. | 
 |  | 
 |         """ | 
 |         # mu = mean, sigma = standard deviation | 
 |  | 
 |         # Uses Kinderman and Monahan method. Reference: Kinderman, | 
 |         # A.J. and Monahan, J.F., "Computer generation of random | 
 |         # variables using the ratio of uniform deviates", ACM Trans | 
 |         # Math Software, 3, (1977), pp257-260. | 
 |  | 
 |         random = self.random | 
 |         while 1: | 
 |             u1 = random() | 
 |             u2 = 1.0 - random() | 
 |             z = NV_MAGICCONST*(u1-0.5)/u2 | 
 |             zz = z*z/4.0 | 
 |             if zz <= -_log(u2): | 
 |                 break | 
 |         return mu + z*sigma | 
 |  | 
 | ## -------------------- lognormal distribution -------------------- | 
 |  | 
 |     def lognormvariate(self, mu, sigma): | 
 |         """Log normal distribution. | 
 |  | 
 |         If you take the natural logarithm of this distribution, you'll get a | 
 |         normal distribution with mean mu and standard deviation sigma. | 
 |         mu can have any value, and sigma must be greater than zero. | 
 |  | 
 |         """ | 
 |         return _exp(self.normalvariate(mu, sigma)) | 
 |  | 
 | ## -------------------- exponential distribution -------------------- | 
 |  | 
 |     def expovariate(self, lambd): | 
 |         """Exponential distribution. | 
 |  | 
 |         lambd is 1.0 divided by the desired mean.  It should be | 
 |         nonzero.  (The parameter would be called "lambda", but that is | 
 |         a reserved word in Python.)  Returned values range from 0 to | 
 |         positive infinity if lambd is positive, and from negative | 
 |         infinity to 0 if lambd is negative. | 
 |  | 
 |         """ | 
 |         # lambd: rate lambd = 1/mean | 
 |         # ('lambda' is a Python reserved word) | 
 |  | 
 |         # we use 1-random() instead of random() to preclude the | 
 |         # possibility of taking the log of zero. | 
 |         return -_log(1.0 - self.random())/lambd | 
 |  | 
 | ## -------------------- von Mises distribution -------------------- | 
 |  | 
 |     def vonmisesvariate(self, mu, kappa): | 
 |         """Circular data distribution. | 
 |  | 
 |         mu is the mean angle, expressed in radians between 0 and 2*pi, and | 
 |         kappa is the concentration parameter, which must be greater than or | 
 |         equal to zero.  If kappa is equal to zero, this distribution reduces | 
 |         to a uniform random angle over the range 0 to 2*pi. | 
 |  | 
 |         """ | 
 |         # mu:    mean angle (in radians between 0 and 2*pi) | 
 |         # kappa: concentration parameter kappa (>= 0) | 
 |         # if kappa = 0 generate uniform random angle | 
 |  | 
 |         # Based upon an algorithm published in: Fisher, N.I., | 
 |         # "Statistical Analysis of Circular Data", Cambridge | 
 |         # University Press, 1993. | 
 |  | 
 |         # Thanks to Magnus Kessler for a correction to the | 
 |         # implementation of step 4. | 
 |  | 
 |         random = self.random | 
 |         if kappa <= 1e-6: | 
 |             return TWOPI * random() | 
 |  | 
 |         a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) | 
 |         b = (a - _sqrt(2.0 * a))/(2.0 * kappa) | 
 |         r = (1.0 + b * b)/(2.0 * b) | 
 |  | 
 |         while 1: | 
 |             u1 = random() | 
 |  | 
 |             z = _cos(_pi * u1) | 
 |             f = (1.0 + r * z)/(r + z) | 
 |             c = kappa * (r - f) | 
 |  | 
 |             u2 = random() | 
 |  | 
 |             if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c): | 
 |                 break | 
 |  | 
 |         u3 = random() | 
 |         if u3 > 0.5: | 
 |             theta = (mu % TWOPI) + _acos(f) | 
 |         else: | 
 |             theta = (mu % TWOPI) - _acos(f) | 
 |  | 
 |         return theta | 
 |  | 
 | ## -------------------- gamma distribution -------------------- | 
 |  | 
 |     def gammavariate(self, alpha, beta): | 
 |         """Gamma distribution.  Not the gamma function! | 
 |  | 
 |         Conditions on the parameters are alpha > 0 and beta > 0. | 
 |  | 
 |         The probability distribution function is: | 
 |  | 
 |                     x ** (alpha - 1) * math.exp(-x / beta) | 
 |           pdf(x) =  -------------------------------------- | 
 |                       math.gamma(alpha) * beta ** alpha | 
 |  | 
 |         """ | 
 |  | 
 |         # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 | 
 |  | 
 |         # Warning: a few older sources define the gamma distribution in terms | 
 |         # of alpha > -1.0 | 
 |         if alpha <= 0.0 or beta <= 0.0: | 
 |             raise ValueError('gammavariate: alpha and beta must be > 0.0') | 
 |  | 
 |         random = self.random | 
 |         if alpha > 1.0: | 
 |  | 
 |             # Uses R.C.H. Cheng, "The generation of Gamma | 
 |             # variables with non-integral shape parameters", | 
 |             # Applied Statistics, (1977), 26, No. 1, p71-74 | 
 |  | 
 |             ainv = _sqrt(2.0 * alpha - 1.0) | 
 |             bbb = alpha - LOG4 | 
 |             ccc = alpha + ainv | 
 |  | 
 |             while 1: | 
 |                 u1 = random() | 
 |                 if not 1e-7 < u1 < .9999999: | 
 |                     continue | 
 |                 u2 = 1.0 - random() | 
 |                 v = _log(u1/(1.0-u1))/ainv | 
 |                 x = alpha*_exp(v) | 
 |                 z = u1*u1*u2 | 
 |                 r = bbb+ccc*v-x | 
 |                 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): | 
 |                     return x * beta | 
 |  | 
 |         elif alpha == 1.0: | 
 |             # expovariate(1) | 
 |             u = random() | 
 |             while u <= 1e-7: | 
 |                 u = random() | 
 |             return -_log(u) * beta | 
 |  | 
 |         else:   # alpha is between 0 and 1 (exclusive) | 
 |  | 
 |             # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle | 
 |  | 
 |             while 1: | 
 |                 u = random() | 
 |                 b = (_e + alpha)/_e | 
 |                 p = b*u | 
 |                 if p <= 1.0: | 
 |                     x = p ** (1.0/alpha) | 
 |                 else: | 
 |                     x = -_log((b-p)/alpha) | 
 |                 u1 = random() | 
 |                 if p > 1.0: | 
 |                     if u1 <= x ** (alpha - 1.0): | 
 |                         break | 
 |                 elif u1 <= _exp(-x): | 
 |                     break | 
 |             return x * beta | 
 |  | 
 | ## -------------------- Gauss (faster alternative) -------------------- | 
 |  | 
 |     def gauss(self, mu, sigma): | 
 |         """Gaussian distribution. | 
 |  | 
 |         mu is the mean, and sigma is the standard deviation.  This is | 
 |         slightly faster than the normalvariate() function. | 
 |  | 
 |         Not thread-safe without a lock around calls. | 
 |  | 
 |         """ | 
 |  | 
 |         # When x and y are two variables from [0, 1), uniformly | 
 |         # distributed, then | 
 |         # | 
 |         #    cos(2*pi*x)*sqrt(-2*log(1-y)) | 
 |         #    sin(2*pi*x)*sqrt(-2*log(1-y)) | 
 |         # | 
 |         # are two *independent* variables with normal distribution | 
 |         # (mu = 0, sigma = 1). | 
 |         # (Lambert Meertens) | 
 |         # (corrected version; bug discovered by Mike Miller, fixed by LM) | 
 |  | 
 |         # Multithreading note: When two threads call this function | 
 |         # simultaneously, it is possible that they will receive the | 
 |         # same return value.  The window is very small though.  To | 
 |         # avoid this, you have to use a lock around all calls.  (I | 
 |         # didn't want to slow this down in the serial case by using a | 
 |         # lock here.) | 
 |  | 
 |         random = self.random | 
 |         z = self.gauss_next | 
 |         self.gauss_next = None | 
 |         if z is None: | 
 |             x2pi = random() * TWOPI | 
 |             g2rad = _sqrt(-2.0 * _log(1.0 - random())) | 
 |             z = _cos(x2pi) * g2rad | 
 |             self.gauss_next = _sin(x2pi) * g2rad | 
 |  | 
 |         return mu + z*sigma | 
 |  | 
 | ## -------------------- beta -------------------- | 
 | ## See | 
 | ## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html | 
 | ## for Ivan Frohne's insightful analysis of why the original implementation: | 
 | ## | 
 | ##    def betavariate(self, alpha, beta): | 
 | ##        # Discrete Event Simulation in C, pp 87-88. | 
 | ## | 
 | ##        y = self.expovariate(alpha) | 
 | ##        z = self.expovariate(1.0/beta) | 
 | ##        return z/(y+z) | 
 | ## | 
 | ## was dead wrong, and how it probably got that way. | 
 |  | 
 |     def betavariate(self, alpha, beta): | 
 |         """Beta distribution. | 
 |  | 
 |         Conditions on the parameters are alpha > 0 and beta > 0. | 
 |         Returned values range between 0 and 1. | 
 |  | 
 |         """ | 
 |  | 
 |         # This version due to Janne Sinkkonen, and matches all the std | 
 |         # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). | 
 |         y = self.gammavariate(alpha, 1.) | 
 |         if y == 0: | 
 |             return 0.0 | 
 |         else: | 
 |             return y / (y + self.gammavariate(beta, 1.)) | 
 |  | 
 | ## -------------------- Pareto -------------------- | 
 |  | 
 |     def paretovariate(self, alpha): | 
 |         """Pareto distribution.  alpha is the shape parameter.""" | 
 |         # Jain, pg. 495 | 
 |  | 
 |         u = 1.0 - self.random() | 
 |         return 1.0 / u ** (1.0/alpha) | 
 |  | 
 | ## -------------------- Weibull -------------------- | 
 |  | 
 |     def weibullvariate(self, alpha, beta): | 
 |         """Weibull distribution. | 
 |  | 
 |         alpha is the scale parameter and beta is the shape parameter. | 
 |  | 
 |         """ | 
 |         # Jain, pg. 499; bug fix courtesy Bill Arms | 
 |  | 
 |         u = 1.0 - self.random() | 
 |         return alpha * (-_log(u)) ** (1.0/beta) | 
 |  | 
 | ## --------------- Operating System Random Source  ------------------ | 
 |  | 
 | class SystemRandom(Random): | 
 |     """Alternate random number generator using sources provided | 
 |     by the operating system (such as /dev/urandom on Unix or | 
 |     CryptGenRandom on Windows). | 
 |  | 
 |      Not available on all systems (see os.urandom() for details). | 
 |     """ | 
 |  | 
 |     def random(self): | 
 |         """Get the next random number in the range [0.0, 1.0).""" | 
 |         return (int.from_bytes(_urandom(7), 'big') >> 3) * RECIP_BPF | 
 |  | 
 |     def getrandbits(self, k): | 
 |         """getrandbits(k) -> x.  Generates a long int with k random bits.""" | 
 |         if k <= 0: | 
 |             raise ValueError('number of bits must be greater than zero') | 
 |         if k != int(k): | 
 |             raise TypeError('number of bits should be an integer') | 
 |         numbytes = (k + 7) // 8                       # bits / 8 and rounded up | 
 |         x = int.from_bytes(_urandom(numbytes), 'big') | 
 |         return x >> (numbytes * 8 - k)                # trim excess bits | 
 |  | 
 |     def seed(self, *args, **kwds): | 
 |         "Stub method.  Not used for a system random number generator." | 
 |         return None | 
 |  | 
 |     def _notimplemented(self, *args, **kwds): | 
 |         "Method should not be called for a system random number generator." | 
 |         raise NotImplementedError('System entropy source does not have state.') | 
 |     getstate = setstate = _notimplemented | 
 |  | 
 | ## -------------------- test program -------------------- | 
 |  | 
 | def _test_generator(n, func, args): | 
 |     import time | 
 |     print(n, 'times', func.__name__) | 
 |     total = 0.0 | 
 |     sqsum = 0.0 | 
 |     smallest = 1e10 | 
 |     largest = -1e10 | 
 |     t0 = time.time() | 
 |     for i in range(n): | 
 |         x = func(*args) | 
 |         total += x | 
 |         sqsum = sqsum + x*x | 
 |         smallest = min(x, smallest) | 
 |         largest = max(x, largest) | 
 |     t1 = time.time() | 
 |     print(round(t1-t0, 3), 'sec,', end=' ') | 
 |     avg = total/n | 
 |     stddev = _sqrt(sqsum/n - avg*avg) | 
 |     print('avg %g, stddev %g, min %g, max %g' % \ | 
 |               (avg, stddev, smallest, largest)) | 
 |  | 
 |  | 
 | def _test(N=2000): | 
 |     _test_generator(N, random, ()) | 
 |     _test_generator(N, normalvariate, (0.0, 1.0)) | 
 |     _test_generator(N, lognormvariate, (0.0, 1.0)) | 
 |     _test_generator(N, vonmisesvariate, (0.0, 1.0)) | 
 |     _test_generator(N, gammavariate, (0.01, 1.0)) | 
 |     _test_generator(N, gammavariate, (0.1, 1.0)) | 
 |     _test_generator(N, gammavariate, (0.1, 2.0)) | 
 |     _test_generator(N, gammavariate, (0.5, 1.0)) | 
 |     _test_generator(N, gammavariate, (0.9, 1.0)) | 
 |     _test_generator(N, gammavariate, (1.0, 1.0)) | 
 |     _test_generator(N, gammavariate, (2.0, 1.0)) | 
 |     _test_generator(N, gammavariate, (20.0, 1.0)) | 
 |     _test_generator(N, gammavariate, (200.0, 1.0)) | 
 |     _test_generator(N, gauss, (0.0, 1.0)) | 
 |     _test_generator(N, betavariate, (3.0, 3.0)) | 
 |     _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0)) | 
 |  | 
 | # Create one instance, seeded from current time, and export its methods | 
 | # as module-level functions.  The functions share state across all uses | 
 | #(both in the user's code and in the Python libraries), but that's fine | 
 | # for most programs and is easier for the casual user than making them | 
 | # instantiate their own Random() instance. | 
 |  | 
 | _inst = Random() | 
 | seed = _inst.seed | 
 | random = _inst.random | 
 | uniform = _inst.uniform | 
 | triangular = _inst.triangular | 
 | randint = _inst.randint | 
 | choice = _inst.choice | 
 | randrange = _inst.randrange | 
 | sample = _inst.sample | 
 | shuffle = _inst.shuffle | 
 | normalvariate = _inst.normalvariate | 
 | lognormvariate = _inst.lognormvariate | 
 | expovariate = _inst.expovariate | 
 | vonmisesvariate = _inst.vonmisesvariate | 
 | gammavariate = _inst.gammavariate | 
 | gauss = _inst.gauss | 
 | betavariate = _inst.betavariate | 
 | paretovariate = _inst.paretovariate | 
 | weibullvariate = _inst.weibullvariate | 
 | getstate = _inst.getstate | 
 | setstate = _inst.setstate | 
 | getrandbits = _inst.getrandbits | 
 |  | 
 | if __name__ == '__main__': | 
 |     _test() |