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