random.py 30.3 KB
Newer Older
1 2
"""Random variable generators.

3 4 5 6 7 8 9
    integers
    --------
           uniform within range

    sequences
    ---------
           pick random element
10
           pick random sample
11 12
           generate random permutation

13 14
    distributions on the real line:
    ------------------------------
15
           uniform
16 17 18 19 20
           normal (Gaussian)
           lognormal
           negative exponential
           gamma
           beta
21 22
           pareto
           Weibull
23 24 25 26 27 28

    distributions on the circle (angles 0 to 2pi)
    ---------------------------------------------
           circular uniform
           von Mises

29 30 31
General notes on the underlying Mersenne Twister core generator:

* The period is 2**19937-1.
32 33 34 35 36 37
* It is one of the most extensively tested generators in existence.
* Without a direct way to compute N steps forward, the semantics of
  jumpahead(n) are weakened to simply jump to another distant state and rely
  on the large period to avoid overlapping sequences.
* The random() method is implemented in C, executes in a single Python step,
  and is, therefore, threadsafe.
38

39
"""
40

41 42
from warnings import warn as _warn
from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
43
from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
44
from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
45 46
from os import urandom as _urandom
from binascii import hexlify as _hexlify
47

48
__all__ = ["Random","seed","random","uniform","randint","choice","sample",
49
           "randrange","shuffle","normalvariate","lognormvariate",
50 51
           "expovariate","vonmisesvariate","gammavariate",
           "gauss","betavariate","paretovariate","weibullvariate",
52
           "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
53
           "SystemRandom"]
54

55 56 57 58
NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
TWOPI = 2.0*_pi
LOG4 = _log(4.0)
SG_MAGICCONST = 1.0 + _log(4.5)
59
BPF = 53        # Number of bits in a float
60
RECIP_BPF = 2**-BPF
61

62

63
# Translated by Guido van Rossum from C source provided by
64
# Adrian Baddeley.  Adapted by Raymond Hettinger for use with
65
# the Mersenne Twister  and os.urandom() core generators.
66

67
import _random
68

69
class Random(_random.Random):
70 71 72 73 74 75 76 77 78 79 80
    """Random number generator base class used by bound module functions.

    Used to instantiate instances of Random to get generators that don't
    share state.  Especially useful for multi-threaded programs, creating
    a different instance of Random for each thread, and using the jumpahead()
    method to ensure that the generated sequences seen by each thread don't
    overlap.

    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(), setstate() and jumpahead().
81 82
    Optionally, implement a getrandombits() method so that randrange()
    can cover arbitrarily large ranges.
83

84
    """
85

86
    VERSION = 3     # used by getstate/setstate
87

88 89
    def __init__(self, x=None):
        """Initialize an instance.
90

91 92
        Optional argument x controls seeding, as for Random.seed().
        """
93

94
        self.seed(x)
95
        self.gauss_next = None
96

97 98
    def seed(self, a=None):
        """Initialize internal state from hashable object.
99

100 101
        None or no argument seeds from current time or from an operating
        system specific randomness source if available.
102

Tim Peters's avatar
Tim Peters committed
103
        If a is not None or an int or long, hash(a) is used instead.
104 105
        """

106
        if a is None:
107
            try:
108
                a = int(_hexlify(_urandom(16)), 16)
109
            except NotImplementedError:
110
                import time
111
                a = int(time.time() * 256) # use fractional seconds
112

113
        super().seed(a)
114 115
        self.gauss_next = None

116 117
    def getstate(self):
        """Return internal state; can be passed to setstate() later."""
118
        return self.VERSION, super().getstate(), self.gauss_next
119 120 121 122

    def setstate(self, state):
        """Restore internal state from object returned by getstate()."""
        version = state[0]
123
        if version == 3:
124
            version, internalstate, self.gauss_next = state
125
            super().setstate(internalstate)
126 127 128 129 130 131 132 133 134 135 136
        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(Random, self).setstate(internalstate)
137 138 139 140 141
        else:
            raise ValueError("state with version %s passed to "
                             "Random.setstate() of version %s" %
                             (version, self.VERSION))

142 143
## ---- Methods below this point do not need to be overridden when
## ---- subclassing for the purpose of using a different core generator.
144

145
## -------------------- pickle support  -------------------
146

147 148
    def __getstate__(self): # for pickle
        return self.getstate()
149

150 151 152
    def __setstate__(self, state):  # for pickle
        self.setstate(state)

153 154 155
    def __reduce__(self):
        return self.__class__, (), self.getstate()

156
## -------------------- integer methods  -------------------
157

158
    def randrange(self, start, stop=None, step=1, int=int, default=None,
159
                  maxwidth=1<<BPF):
160 161 162 163
        """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.
164
        Do not supply the 'int', 'default', and 'maxwidth' arguments.
165 166 167
        """

        # This code is a bit messy to make it fast for the
168
        # common case while still doing adequate error checking.
169 170
        istart = int(start)
        if istart != start:
171
            raise ValueError("non-integer arg 1 for randrange()")
172 173
        if stop is default:
            if istart > 0:
174 175
                if istart >= maxwidth:
                    return self._randbelow(istart)
176
                return int(self.random() * istart)
177
            raise ValueError("empty range for randrange()")
178 179

        # stop argument supplied.
180 181
        istop = int(stop)
        if istop != stop:
182
            raise ValueError("non-integer stop for randrange()")
183 184
        width = istop - istart
        if step == 1 and width > 0:
185
            # Note that
186
            #     int(istart + self.random()*width)
187 188 189 190 191
            # instead would be incorrect.  For example, consider istart
            # = -2 and istop = 0.  Then the guts would be in
            # -2.0 to 0.0 exclusive on both ends (ignoring that random()
            # might return 0.0), and because int() truncates toward 0, the
            # final result would be -1 or 0 (instead of -2 or -1).
192
            #     istart + int(self.random()*width)
193 194 195 196
            # would also be incorrect, for a subtler reason:  the RHS
            # can return a long, and then randrange() would also return
            # a long, but we're supposed to return an int (for backward
            # compatibility).
197 198

            if width >= maxwidth:
Tim Peters's avatar
Tim Peters committed
199
                return int(istart + self._randbelow(width))
200
            return int(istart + int(self.random()*width))
201
        if step == 1:
202
            raise ValueError("empty range for randrange() (%d,%d, %d)" % (istart, istop, width))
203 204

        # Non-unit step argument supplied.
205 206
        istep = int(step)
        if istep != step:
207
            raise ValueError("non-integer step for randrange()")
208
        if istep > 0:
209
            n = (width + istep - 1) // istep
210
        elif istep < 0:
211
            n = (width + istep + 1) // istep
212
        else:
213
            raise ValueError("zero step for randrange()")
214 215

        if n <= 0:
216
            raise ValueError("empty range for randrange()")
217 218

        if n >= maxwidth:
219
            return istart + istep*self._randbelow(n)
220 221 222
        return istart + istep*int(self.random() * n)

    def randint(self, a, b):
223
        """Return random integer in range [a, b], including both end points.
224 225 226 227
        """

        return self.randrange(a, b+1)

228
    def _randbelow(self, n, _log=_log, int=int, _maxwidth=1<<BPF,
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
                   _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
        """Return a random int in the range [0,n)

        Handles the case where n has more bits than returned
        by a single call to the underlying generator.
        """

        try:
            getrandbits = self.getrandbits
        except AttributeError:
            pass
        else:
            # Only call self.getrandbits if the original random() builtin method
            # has not been overridden or if a new getrandbits() was supplied.
            # This assures that the two methods correspond.
            if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
                k = int(1.00001 + _log(n-1, 2.0))   # 2**k > n-1 > 2**(k-2)
                r = getrandbits(k)
                while r >= n:
                    r = getrandbits(k)
                return r
        if n >= _maxwidth:
            _warn("Underlying random() generator does not supply \n"
                "enough bits to choose from a population range this large")
        return int(self.random() * n)

255 256
## -------------------- sequence methods  -------------------

257 258
    def choice(self, seq):
        """Choose a random element from a non-empty sequence."""
Raymond Hettinger's avatar
Raymond Hettinger committed
259
        return seq[int(self.random() * len(seq))]  # raises IndexError if seq is empty
260 261 262 263 264 265 266 267 268 269

    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.
        """

        if random is None:
            random = self.random
270
        for i in reversed(range(1, len(x))):
271
            # pick an element in x[:i+1] with which to exchange x[i]
272 273 274
            j = int(random() * (i+1))
            x[i], x[j] = x[j], x[i]

275
    def sample(self, population, k):
276 277
        """Chooses k unique random elements from a population sequence.

278 279 280 281 282
        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).
283

284 285 286
        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.
287

288
        To choose a sample in a range of integers, use range as an argument.
289
        This is especially fast and space efficient for sampling from a
290
        large population:   sample(range(10000000), 60)
291 292
        """

293 294 295 296 297 298 299 300 301
        # XXX Although the documentation says `population` is "a sequence",
        # XXX attempts are made to cater to any iterable with a __len__
        # XXX method.  This has had mixed success.  Examples from both
        # XXX sides:  sets work fine, and should become officially supported;
        # XXX dicts are much harder, and have failed in various subtle
        # XXX ways across attempts.  Support for mapping types should probably
        # XXX be dropped (and users should pass mapping.keys() or .values()
        # XXX explicitly).

302
        # Sampling without replacement entails tracking either potential
303
        # selections (the pool) in a list or previous selections in a set.
304

Jeremy Hylton's avatar
Jeremy Hylton committed
305 306
        # When the number of selections is small compared to the
        # population, then tracking selections is efficient, requiring
307
        # only a small set and an occasional reselection.  For
Jeremy Hylton's avatar
Jeremy Hylton committed
308 309
        # a larger number of selections, the pool tracking method is
        # preferred since the list takes less space than the
310
        # set and it doesn't suffer from frequent reselections.
311

312 313
        n = len(population)
        if not 0 <= k <= n:
314
            raise ValueError("sample larger than population")
315
        random = self.random
316
        _int = int
317
        result = [None] * k
318 319
        setsize = 21        # size of a small set minus size of an empty list
        if k > 5:
320
            setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
321 322 323
        if n <= setsize or hasattr(population, "keys"):
            # An n-length list is smaller than a k-length set, or this is a
            # mapping type so the other algorithm wouldn't work.
324
            pool = list(population)
325
            for i in range(k):         # invariant:  non-selected at [0,n-i)
326
                j = _int(random() * (n-i))
327
                result[i] = pool[j]
328
                pool[j] = pool[n-i-1]   # move non-selected item into vacancy
329
        else:
330
            try:
331 332
                selected = set()
                selected_add = selected.add
333
                for i in range(k):
334
                    j = _int(random() * n)
335 336 337 338 339 340 341 342
                    while j in selected:
                        j = _int(random() * n)
                    selected_add(j)
                    result[i] = population[j]
            except (TypeError, KeyError):   # handle (at least) sets
                if isinstance(population, list):
                    raise
                return self.sample(tuple(population), k)
343
        return result
344

345 346 347
## -------------------- real-valued distributions  -------------------

## -------------------- uniform distribution -------------------
348 349 350 351

    def uniform(self, a, b):
        """Get a random number in the range [a, b)."""
        return a + (b-a) * self.random()
352

353
## -------------------- normal distribution --------------------
354

355
    def normalvariate(self, mu, sigma):
356 357 358
        """Normal distribution.

        mu is the mean, and sigma is the standard deviation.
359

360
        """
361 362 363 364 365 366 367 368
        # 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
369
        while 1:
370
            u1 = random()
371
            u2 = 1.0 - random()
372 373 374 375 376
            z = NV_MAGICCONST*(u1-0.5)/u2
            zz = z*z/4.0
            if zz <= -_log(u2):
                break
        return mu + z*sigma
377

378
## -------------------- lognormal distribution --------------------
379

380
    def lognormvariate(self, mu, sigma):
381 382 383 384 385
        """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.
386

387
        """
388
        return _exp(self.normalvariate(mu, sigma))
389

390
## -------------------- exponential distribution --------------------
391

392
    def expovariate(self, lambd):
393 394 395 396 397
        """Exponential distribution.

        lambd is 1.0 divided by the desired mean.  (The parameter would be
        called "lambda", but that is a reserved word in Python.)  Returned
        values range from 0 to positive infinity.
398

399
        """
400 401
        # lambd: rate lambd = 1/mean
        # ('lambda' is a Python reserved word)
402

403
        random = self.random
Tim Peters's avatar
Tim Peters committed
404
        u = random()
405 406 407
        while u <= 1e-7:
            u = random()
        return -_log(u)/lambd
408

409
## -------------------- von Mises distribution --------------------
410

411
    def vonmisesvariate(self, mu, kappa):
412
        """Circular data distribution.
413

414 415 416 417
        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.
418

419
        """
420 421 422
        # mu:    mean angle (in radians between 0 and 2*pi)
        # kappa: concentration parameter kappa (>= 0)
        # if kappa = 0 generate uniform random angle
423

424 425 426
        # Based upon an algorithm published in: Fisher, N.I.,
        # "Statistical Analysis of Circular Data", Cambridge
        # University Press, 1993.
427

428 429
        # Thanks to Magnus Kessler for a correction to the
        # implementation of step 4.
430

431 432 433
        random = self.random
        if kappa <= 1e-6:
            return TWOPI * random()
434

435 436 437
        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)
438

439
        while 1:
440
            u1 = random()
441

442 443 444
            z = _cos(_pi * u1)
            f = (1.0 + r * z)/(r + z)
            c = kappa * (r - f)
445

446
            u2 = random()
447

448
            if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
449
                break
450

451 452 453 454 455
        u3 = random()
        if u3 > 0.5:
            theta = (mu % TWOPI) + _acos(f)
        else:
            theta = (mu % TWOPI) - _acos(f)
456

457
        return theta
458

459
## -------------------- gamma distribution --------------------
460

461
    def gammavariate(self, alpha, beta):
462 463 464 465 466
        """Gamma distribution.  Not the gamma function!

        Conditions on the parameters are alpha > 0 and beta > 0.

        """
Tim Peters's avatar
Tim Peters committed
467

468
        # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters's avatar
Tim Peters committed
469

470 471 472
        # Warning: a few older sources define the gamma distribution in terms
        # of alpha > -1.0
        if alpha <= 0.0 or beta <= 0.0:
473
            raise ValueError('gammavariate: alpha and beta must be > 0.0')
Tim Peters's avatar
Tim Peters committed
474

475 476 477 478 479 480 481
        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

482 483 484
            ainv = _sqrt(2.0 * alpha - 1.0)
            bbb = alpha - LOG4
            ccc = alpha + ainv
Tim Peters's avatar
Tim Peters committed
485

486
            while 1:
487
                u1 = random()
488 489 490
                if not 1e-7 < u1 < .9999999:
                    continue
                u2 = 1.0 - random()
491 492 493 494 495
                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):
496
                    return x * beta
497 498 499

        elif alpha == 1.0:
            # expovariate(1)
Tim Peters's avatar
Tim Peters committed
500
            u = random()
501 502
            while u <= 1e-7:
                u = random()
503
            return -_log(u) * beta
504 505 506 507 508

        else:   # alpha is between 0 and 1 (exclusive)

            # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

509
            while 1:
510 511 512 513
                u = random()
                b = (_e + alpha)/_e
                p = b*u
                if p <= 1.0:
514
                    x = p ** (1.0/alpha)
515 516 517
                else:
                    x = -_log((b-p)/alpha)
                u1 = random()
518 519 520 521
                if p > 1.0:
                    if u1 <= x ** (alpha - 1.0):
                        break
                elif u1 <= _exp(-x):
522
                    break
523 524
            return x * beta

525
## -------------------- Gauss (faster alternative) --------------------
526

527
    def gauss(self, mu, sigma):
528 529 530 531 532 533
        """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.
534

535
        """
536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564

        # 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
565

566
## -------------------- beta --------------------
567 568 569 570 571 572 573 574 575 576 577 578
## See
## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
## 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.
579

580
    def betavariate(self, alpha, beta):
581 582
        """Beta distribution.

583
        Conditions on the parameters are alpha > 0 and beta > 0.
584
        Returned values range between 0 and 1.
585

586
        """
587

588 589 590 591 592 593 594
        # 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.))
595

596
## -------------------- Pareto --------------------
597

598
    def paretovariate(self, alpha):
599
        """Pareto distribution.  alpha is the shape parameter."""
600
        # Jain, pg. 495
601

602
        u = 1.0 - self.random()
603
        return 1.0 / pow(u, 1.0/alpha)
604

605
## -------------------- Weibull --------------------
606

607
    def weibullvariate(self, alpha, beta):
608 609 610
        """Weibull distribution.

        alpha is the scale parameter and beta is the shape parameter.
611

612
        """
613
        # Jain, pg. 499; bug fix courtesy Bill Arms
614

615
        u = 1.0 - self.random()
616
        return alpha * pow(-_log(u), 1.0/beta)
617

618 619 620 621 622 623 624 625 626
## -------------------- Wichmann-Hill -------------------

class WichmannHill(Random):

    VERSION = 1     # used by getstate/setstate

    def seed(self, a=None):
        """Initialize internal state from hashable object.

627 628
        None or no argument seeds from current time or from an operating
        system specific randomness source if available.
629 630 631 632 633 634 635 636 637 638

        If a is not None or an int or long, hash(a) is used instead.

        If a is an int or long, a is used directly.  Distinct values between
        0 and 27814431486575L inclusive are guaranteed to yield distinct
        internal states (this guarantee is specific to the default
        Wichmann-Hill generator).
        """

        if a is None:
639
            try:
640
                a = int(_hexlify(_urandom(16)), 16)
641
            except NotImplementedError:
642
                import time
643
                a = int(time.time() * 256) # use fractional seconds
644

645
        if not isinstance(a, int):
646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734
            a = hash(a)

        a, x = divmod(a, 30268)
        a, y = divmod(a, 30306)
        a, z = divmod(a, 30322)
        self._seed = int(x)+1, int(y)+1, int(z)+1

        self.gauss_next = None

    def random(self):
        """Get the next random number in the range [0.0, 1.0)."""

        # Wichman-Hill random number generator.
        #
        # Wichmann, B. A. & Hill, I. D. (1982)
        # Algorithm AS 183:
        # An efficient and portable pseudo-random number generator
        # Applied Statistics 31 (1982) 188-190
        #
        # see also:
        #        Correction to Algorithm AS 183
        #        Applied Statistics 33 (1984) 123
        #
        #        McLeod, A. I. (1985)
        #        A remark on Algorithm AS 183
        #        Applied Statistics 34 (1985),198-200

        # This part is thread-unsafe:
        # BEGIN CRITICAL SECTION
        x, y, z = self._seed
        x = (171 * x) % 30269
        y = (172 * y) % 30307
        z = (170 * z) % 30323
        self._seed = x, y, z
        # END CRITICAL SECTION

        # Note:  on a platform using IEEE-754 double arithmetic, this can
        # never return 0.0 (asserted by Tim; proof too long for a comment).
        return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0

    def getstate(self):
        """Return internal state; can be passed to setstate() later."""
        return self.VERSION, self._seed, self.gauss_next

    def setstate(self, state):
        """Restore internal state from object returned by getstate()."""
        version = state[0]
        if version == 1:
            version, self._seed, self.gauss_next = state
        else:
            raise ValueError("state with version %s passed to "
                             "Random.setstate() of version %s" %
                             (version, self.VERSION))

    def jumpahead(self, n):
        """Act as if n calls to random() were made, but quickly.

        n is an int, greater than or equal to 0.

        Example use:  If you have 2 threads and know that each will
        consume no more than a million random numbers, create two Random
        objects r1 and r2, then do
            r2.setstate(r1.getstate())
            r2.jumpahead(1000000)
        Then r1 and r2 will use guaranteed-disjoint segments of the full
        period.
        """

        if not n >= 0:
            raise ValueError("n must be >= 0")
        x, y, z = self._seed
        x = int(x * pow(171, n, 30269)) % 30269
        y = int(y * pow(172, n, 30307)) % 30307
        z = int(z * pow(170, n, 30323)) % 30323
        self._seed = x, y, z

    def __whseed(self, x=0, y=0, z=0):
        """Set the Wichmann-Hill seed from (x, y, z).

        These must be integers in the range [0, 256).
        """

        if not type(x) == type(y) == type(z) == int:
            raise TypeError('seeds must be integers')
        if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
            raise ValueError('seeds must be in range(0, 256)')
        if 0 == x == y == z:
            # Initialize from current time
            import time
735
            t = int(time.time() * 256)
736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767
            t = int((t&0xffffff) ^ (t>>24))
            t, x = divmod(t, 256)
            t, y = divmod(t, 256)
            t, z = divmod(t, 256)
        # Zero is a poor seed, so substitute 1
        self._seed = (x or 1, y or 1, z or 1)

        self.gauss_next = None

    def whseed(self, a=None):
        """Seed from hashable object's hash code.

        None or no argument seeds from current time.  It is not guaranteed
        that objects with distinct hash codes lead to distinct internal
        states.

        This is obsolete, provided for compatibility with the seed routine
        used prior to Python 2.1.  Use the .seed() method instead.
        """

        if a is None:
            self.__whseed()
            return
        a = hash(a)
        a, x = divmod(a, 256)
        a, y = divmod(a, 256)
        a, z = divmod(a, 256)
        x = (x + a) % 256 or 1
        y = (y + a) % 256 or 1
        z = (z + a) % 256 or 1
        self.__whseed(x, y, z)

768
## --------------- Operating System Random Source  ------------------
769

770 771 772 773
class SystemRandom(Random):
    """Alternate random number generator using sources provided
    by the operating system (such as /dev/urandom on Unix or
    CryptGenRandom on Windows).
774 775 776 777 778 779

     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)."""
780
        return (int(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
781 782 783 784 785 786 787 788

    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')
        bytes = (k + 7) // 8                    # bits / 8 and rounded up
789
        x = int(_hexlify(_urandom(bytes)), 16)
790 791 792
        return x >> (bytes * 8 - k)             # trim excess bits

    def _stub(self, *args, **kwds):
793
        "Stub method.  Not used for a system random number generator."
794 795 796 797
        return None
    seed = jumpahead = _stub

    def _notimplemented(self, *args, **kwds):
798 799
        "Method should not be called for a system random number generator."
        raise NotImplementedError('System entropy source does not have state.')
800 801
    getstate = setstate = _notimplemented

802
## -------------------- test program --------------------
803

804
def _test_generator(n, func, args):
Tim Peters's avatar
Tim Peters committed
805
    import time
806
    print(n, 'times', func.__name__)
807
    total = 0.0
Tim Peters's avatar
Tim Peters committed
808 809 810 811 812
    sqsum = 0.0
    smallest = 1e10
    largest = -1e10
    t0 = time.time()
    for i in range(n):
813
        x = func(*args)
814
        total += x
Tim Peters's avatar
Tim Peters committed
815 816 817 818
        sqsum = sqsum + x*x
        smallest = min(x, smallest)
        largest = max(x, largest)
    t1 = time.time()
819
    print(round(t1-t0, 3), 'sec,', end=' ')
820
    avg = total/n
821
    stddev = _sqrt(sqsum/n - avg*avg)
822 823
    print('avg %g, stddev %g, min %g, max %g' % \
              (avg, stddev, smallest, largest))
824

825 826

def _test(N=2000):
827 828 829 830 831 832 833 834 835 836 837 838 839 840 841
    _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))
842

843
# Create one instance, seeded from current time, and export its methods
844 845 846 847 848
# 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.

849 850 851 852 853 854 855
_inst = Random()
seed = _inst.seed
random = _inst.random
uniform = _inst.uniform
randint = _inst.randint
choice = _inst.choice
randrange = _inst.randrange
856
sample = _inst.sample
857 858 859 860 861 862 863 864 865 866 867 868
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
869
jumpahead = _inst.jumpahead
870
getrandbits = _inst.getrandbits
871

872
if __name__ == '__main__':
873
    _test()