random.py 27 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 32 33 34 35 36 37 38 39
General notes on the underlying Mersenne Twister core generator:

* The period is 2**19937-1.
* It is one of the most extensively tested generators in existence
* 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.

40
"""
41

42 43
from math import log as _log, exp as _exp, pi as _pi, e as _e
from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
44
from math import floor as _floor
45

46
__all__ = ["Random","seed","random","uniform","randint","choice","sample",
47 48 49
           "randrange","shuffle","normalvariate","lognormvariate",
           "cunifvariate","expovariate","vonmisesvariate","gammavariate",
           "stdgamma","gauss","betavariate","paretovariate","weibullvariate",
50
           "getstate","setstate","jumpahead"]
51

52 53 54 55
NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
TWOPI = 2.0*_pi
LOG4 = _log(4.0)
SG_MAGICCONST = 1.0 + _log(4.5)
56

57
# Translated by Guido van Rossum from C source provided by
58 59
# Adrian Baddeley.  Adapted by Raymond Hettinger for use with
# the Mersenne Twister core generator.
60

61
import _random
62

63
class Random(_random.Random):
64 65 66 67 68 69 70 71 72 73 74
    """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().
75

76
    """
77

78
    VERSION = 2     # used by getstate/setstate
79

80 81
    def __init__(self, x=None):
        """Initialize an instance.
82

83 84
        Optional argument x controls seeding, as for Random.seed().
        """
85

86
        self.seed(x)
87
        self.gauss_next = None
88

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

92 93
        None or no argument seeds from current time.

Tim Peters's avatar
Tim Peters committed
94
        If a is not None or an int or long, hash(a) is used instead.
95 96
        """

97
        super(Random, self).seed(a)
98 99
        self.gauss_next = None

100 101
    def getstate(self):
        """Return internal state; can be passed to setstate() later."""
102
        return self.VERSION, super(Random, self).getstate(), self.gauss_next
103 104 105 106

    def setstate(self, state):
        """Restore internal state from object returned by getstate()."""
        version = state[0]
107 108
        if version == 2:
            version, internalstate, self.gauss_next = state
109
            super(Random, self).setstate(internalstate)
110 111 112 113 114
        else:
            raise ValueError("state with version %s passed to "
                             "Random.setstate() of version %s" %
                             (version, self.VERSION))

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

118
## -------------------- pickle support  -------------------
119

120 121
    def __getstate__(self): # for pickle
        return self.getstate()
122

123 124 125
    def __setstate__(self, state):  # for pickle
        self.setstate(state)

126 127 128
    def __reduce__(self):
        return self.__class__, (), self.getstate()

129
## -------------------- integer methods  -------------------
130 131 132 133 134 135 136 137 138 139

    def randrange(self, start, stop=None, step=1, int=int, default=None):
        """Choose a random item from range(start, stop[, step]).

        This fixes the problem with randint() which includes the
        endpoint; in Python this is usually not what you want.
        Do not supply the 'int' and 'default' arguments.
        """

        # This code is a bit messy to make it fast for the
140
        # common case while still doing adequate error checking.
141 142 143 144 145 146 147
        istart = int(start)
        if istart != start:
            raise ValueError, "non-integer arg 1 for randrange()"
        if stop is default:
            if istart > 0:
                return int(self.random() * istart)
            raise ValueError, "empty range for randrange()"
148 149

        # stop argument supplied.
150 151 152
        istop = int(stop)
        if istop != stop:
            raise ValueError, "non-integer stop for randrange()"
153
        if step == 1 and istart < istop:
154 155 156 157 158 159 160 161 162 163 164 165 166
            # Note that
            #     int(istart + self.random()*(istop - istart))
            # 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).
            #     istart + int(self.random()*(istop - istart))
            # 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).
            return int(istart + int(self.random()*(istop - istart)))
167 168
        if step == 1:
            raise ValueError, "empty range for randrange()"
169 170

        # Non-unit step argument supplied.
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
        istep = int(step)
        if istep != step:
            raise ValueError, "non-integer step for randrange()"
        if istep > 0:
            n = (istop - istart + istep - 1) / istep
        elif istep < 0:
            n = (istop - istart + istep + 1) / istep
        else:
            raise ValueError, "zero step for randrange()"

        if n <= 0:
            raise ValueError, "empty range for randrange()"
        return istart + istep*int(self.random() * n)

    def randint(self, a, b):
186
        """Return random integer in range [a, b], including both end points.
187 188 189 190
        """

        return self.randrange(a, b+1)

191 192
## -------------------- sequence methods  -------------------

193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
    def choice(self, seq):
        """Choose a random element from a non-empty sequence."""
        return seq[int(self.random() * len(seq))]

    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.

        Note that for even rather small len(x), the total number of
        permutations of x is larger than the period of most random number
        generators; this implies that "most" permutations of a long
        sequence can never be generated.
        """

        if random is None:
            random = self.random
        for i in xrange(len(x)-1, 0, -1):
212
            # pick an element in x[:i+1] with which to exchange x[i]
213 214 215
            j = int(random() * (i+1))
            x[i], x[j] = x[j], x[i]

216
    def sample(self, population, k):
217 218
        """Chooses k unique random elements from a population sequence.

219 220 221 222 223
        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).
224

225 226 227
        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.
228

229 230 231
        To choose a sample in a range of integers, use xrange as an argument.
        This is especially fast and space efficient for sampling from a
        large population:   sample(xrange(10000000), 60)
232 233
        """

234
        # Sampling without replacement entails tracking either potential
235 236
        # selections (the pool) in a list or previous selections in a
        # dictionary.
237

238 239 240 241 242 243
        # When the number of selections is small compared to the population,
        # then tracking selections is efficient, requiring only a small
        # dictionary and an occasional reselection.  For a larger number of
        # selections, the pool tracking method is preferred since the list takes
        # less space than the dictionary and it doesn't suffer from frequent
        # reselections.
244

245 246 247
        n = len(population)
        if not 0 <= k <= n:
            raise ValueError, "sample larger than population"
248
        random = self.random
249
        _int = int
250
        result = [None] * k
251
        if n < 6 * k:     # if n len list takes less space than a k len dict
252 253
            pool = list(population)
            for i in xrange(k):         # invariant:  non-selected at [0,n-i)
254
                j = _int(random() * (n-i))
255
                result[i] = pool[j]
256
                pool[j] = pool[n-i-1]   # move non-selected item into vacancy
257
        else:
258
            selected = {}
259
            for i in xrange(k):
260
                j = _int(random() * n)
261
                while j in selected:
262
                    j = _int(random() * n)
263
                result[i] = selected[j] = population[j]
264
        return result
265

266 267 268
## -------------------- real-valued distributions  -------------------

## -------------------- uniform distribution -------------------
269 270 271 272

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

274
## -------------------- normal distribution --------------------
275

276
    def normalvariate(self, mu, sigma):
277 278 279
        """Normal distribution.

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

281
        """
282 283 284 285 286 287 288 289
        # 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
290
        while True:
291
            u1 = random()
292
            u2 = 1.0 - random()
293 294 295 296 297
            z = NV_MAGICCONST*(u1-0.5)/u2
            zz = z*z/4.0
            if zz <= -_log(u2):
                break
        return mu + z*sigma
298

299
## -------------------- lognormal distribution --------------------
300

301
    def lognormvariate(self, mu, sigma):
302 303 304 305 306
        """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.
307

308
        """
309
        return _exp(self.normalvariate(mu, sigma))
310

311
## -------------------- circular uniform --------------------
312

313
    def cunifvariate(self, mean, arc):
314 315 316 317 318 319 320 321 322
        """Circular uniform distribution.

        mean is the mean angle, and arc is the range of the distribution,
        centered around the mean angle.  Both values must be expressed in
        radians.  Returned values range between mean - arc/2 and
        mean + arc/2 and are normalized to between 0 and pi.

        Deprecated in version 2.3.  Use:
            (mean + arc * (Random.random() - 0.5)) % Math.pi
323

324
        """
325 326
        # mean: mean angle (in radians between 0 and pi)
        # arc:  range of distribution (in radians between 0 and pi)
327 328 329 330
        import warnings
        warnings.warn("The cunifvariate function is deprecated; Use (mean "
                      "+ arc * (Random.random() - 0.5)) % Math.pi instead",
                      DeprecationWarning)
331

332
        return (mean + arc * (self.random() - 0.5)) % _pi
333

334
## -------------------- exponential distribution --------------------
335

336
    def expovariate(self, lambd):
337 338 339 340 341
        """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.
342

343
        """
344 345
        # lambd: rate lambd = 1/mean
        # ('lambda' is a Python reserved word)
346

347
        random = self.random
Tim Peters's avatar
Tim Peters committed
348
        u = random()
349 350 351
        while u <= 1e-7:
            u = random()
        return -_log(u)/lambd
352

353
## -------------------- von Mises distribution --------------------
354

355
    def vonmisesvariate(self, mu, kappa):
356
        """Circular data distribution.
357

358 359 360 361
        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.
362

363
        """
364 365 366
        # mu:    mean angle (in radians between 0 and 2*pi)
        # kappa: concentration parameter kappa (>= 0)
        # if kappa = 0 generate uniform random angle
367

368 369 370
        # Based upon an algorithm published in: Fisher, N.I.,
        # "Statistical Analysis of Circular Data", Cambridge
        # University Press, 1993.
371

372 373
        # Thanks to Magnus Kessler for a correction to the
        # implementation of step 4.
374

375 376 377
        random = self.random
        if kappa <= 1e-6:
            return TWOPI * random()
378

379 380 381
        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)
382

383
        while True:
384
            u1 = random()
385

386 387 388
            z = _cos(_pi * u1)
            f = (1.0 + r * z)/(r + z)
            c = kappa * (r - f)
389

390
            u2 = random()
391

392 393
            if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
                break
394

395 396 397 398 399
        u3 = random()
        if u3 > 0.5:
            theta = (mu % TWOPI) + _acos(f)
        else:
            theta = (mu % TWOPI) - _acos(f)
400

401
        return theta
402

403
## -------------------- gamma distribution --------------------
404

405
    def gammavariate(self, alpha, beta):
406 407 408 409 410
        """Gamma distribution.  Not the gamma function!

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

        """
Tim Peters's avatar
Tim Peters committed
411

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

414 415 416 417
        # Warning: a few older sources define the gamma distribution in terms
        # of alpha > -1.0
        if alpha <= 0.0 or beta <= 0.0:
            raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
Tim Peters's avatar
Tim Peters committed
418

419 420 421 422 423 424 425
        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

426 427 428
            ainv = _sqrt(2.0 * alpha - 1.0)
            bbb = alpha - LOG4
            ccc = alpha + ainv
Tim Peters's avatar
Tim Peters committed
429

430
            while True:
431
                u1 = random()
432 433 434
                if not 1e-7 < u1 < .9999999:
                    continue
                u2 = 1.0 - random()
435 436 437 438 439
                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):
440
                    return x * beta
441 442 443

        elif alpha == 1.0:
            # expovariate(1)
Tim Peters's avatar
Tim Peters committed
444
            u = random()
445 446
            while u <= 1e-7:
                u = random()
447
            return -_log(u) * beta
448 449 450 451 452

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

            # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

453
            while True:
454 455 456 457 458 459 460 461 462 463 464 465
                u = random()
                b = (_e + alpha)/_e
                p = b*u
                if p <= 1.0:
                    x = pow(p, 1.0/alpha)
                else:
                    # p > 1
                    x = -_log((b-p)/alpha)
                u1 = random()
                if not (((p <= 1.0) and (u1 > _exp(-x))) or
                          ((p > 1)  and  (u1 > pow(x, alpha - 1.0)))):
                    break
466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486
            return x * beta


    def stdgamma(self, alpha, ainv, bbb, ccc):
        # This method was (and shall remain) undocumented.
        # This method is deprecated
        # for the following reasons:
        # 1. Returns same as .gammavariate(alpha, 1.0)
        # 2. Requires caller to provide 3 extra arguments
        #    that are functions of alpha anyway
        # 3. Can't be used for alpha < 0.5

        # ainv = sqrt(2 * alpha - 1)
        # bbb = alpha - log(4)
        # ccc = alpha + ainv
        import warnings
        warnings.warn("The stdgamma function is deprecated; "
                      "use gammavariate() instead",
                      DeprecationWarning)
        return self.gammavariate(alpha, 1.0)

487

488

489
## -------------------- Gauss (faster alternative) --------------------
490

491
    def gauss(self, mu, sigma):
492 493 494 495 496 497
        """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.
498

499
        """
500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528

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

530
## -------------------- beta --------------------
531 532 533 534 535 536 537 538 539 540 541 542
## 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.
543

544
    def betavariate(self, alpha, beta):
545 546 547 548
        """Beta distribution.

        Conditions on the parameters are alpha > -1 and beta} > -1.
        Returned values range between 0 and 1.
549

550
        """
551

552 553 554 555 556 557 558
        # 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.))
559

560
## -------------------- Pareto --------------------
561

562
    def paretovariate(self, alpha):
563
        """Pareto distribution.  alpha is the shape parameter."""
564
        # Jain, pg. 495
565

566
        u = 1.0 - self.random()
567
        return 1.0 / pow(u, 1.0/alpha)
568

569
## -------------------- Weibull --------------------
570

571
    def weibullvariate(self, alpha, beta):
572 573 574
        """Weibull distribution.

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

576
        """
577
        # Jain, pg. 499; bug fix courtesy Bill Arms
578

579
        u = 1.0 - self.random()
580
        return alpha * pow(-_log(u), 1.0/beta)
581

582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 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
## -------------------- Wichmann-Hill -------------------

class WichmannHill(Random):

    VERSION = 1     # used by getstate/setstate

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

        None or no argument seeds from current time.

        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:
            # Initialize from current time
            import time
            a = long(time.time() * 256)

        if not isinstance(a, (int, long)):
            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
            t = long(time.time() * 256)
            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)

729
## -------------------- test program --------------------
730

731
def _test_generator(n, funccall):
Tim Peters's avatar
Tim Peters committed
732 733 734
    import time
    print n, 'times', funccall
    code = compile(funccall, funccall, 'eval')
735
    total = 0.0
Tim Peters's avatar
Tim Peters committed
736 737 738 739 740 741
    sqsum = 0.0
    smallest = 1e10
    largest = -1e10
    t0 = time.time()
    for i in range(n):
        x = eval(code)
742
        total += x
Tim Peters's avatar
Tim Peters committed
743 744 745 746 747
        sqsum = sqsum + x*x
        smallest = min(x, smallest)
        largest = max(x, largest)
    t1 = time.time()
    print round(t1-t0, 3), 'sec,',
748
    avg = total/n
749
    stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters's avatar
Tim Peters committed
750 751
    print 'avg %g, stddev %g, min %g, max %g' % \
              (avg, stddev, smallest, largest)
752

753 754

def _test(N=2000):
755 756 757 758 759
    _test_generator(N, 'random()')
    _test_generator(N, 'normalvariate(0.0, 1.0)')
    _test_generator(N, 'lognormvariate(0.0, 1.0)')
    _test_generator(N, 'cunifvariate(0.0, 1.0)')
    _test_generator(N, 'vonmisesvariate(0.0, 1.0)')
760 761
    _test_generator(N, 'gammavariate(0.01, 1.0)')
    _test_generator(N, 'gammavariate(0.1, 1.0)')
Tim Peters's avatar
Tim Peters committed
762
    _test_generator(N, 'gammavariate(0.1, 2.0)')
763 764 765 766 767 768 769 770
    _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)')
771

772
# Create one instance, seeded from current time, and export its methods
773 774 775 776 777
# 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.

778 779 780 781 782 783 784
_inst = Random()
seed = _inst.seed
random = _inst.random
uniform = _inst.uniform
randint = _inst.randint
choice = _inst.choice
randrange = _inst.randrange
785
sample = _inst.sample
786 787 788 789 790 791 792 793 794 795 796 797 798 799
shuffle = _inst.shuffle
normalvariate = _inst.normalvariate
lognormvariate = _inst.lognormvariate
cunifvariate = _inst.cunifvariate
expovariate = _inst.expovariate
vonmisesvariate = _inst.vonmisesvariate
gammavariate = _inst.gammavariate
stdgamma = _inst.stdgamma
gauss = _inst.gauss
betavariate = _inst.betavariate
paretovariate = _inst.paretovariate
weibullvariate = _inst.weibullvariate
getstate = _inst.getstate
setstate = _inst.setstate
800
jumpahead = _inst.jumpahead
801

802
if __name__ == '__main__':
803
    _test()