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