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

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

    sequences
    ---------
           pick random element
           generate random permutation

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

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

Translated from anonymously contributed C/C++ source.

28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
Multi-threading note:  the random number generator used here is not thread-
safe; it is possible that two calls return the same random value.  However,
you can instantiate a different instance of Random() in each thread to get
generators that don't share state, then use .setstate() and .jumpahead() to
move the generators to disjoint segments of the full period.  For example,

def create_generators(num, delta, firstseed=None):
    ""\"Return list of num distinct generators.
    Each generator has its own unique segment of delta elements from
    Random.random()'s full period.
    Seed the first generator with optional arg firstseed (default is
    None, to seed from current time).
    ""\"

    from random import Random
    g = Random(firstseed)
    result = [g]
    for i in range(num - 1):
        laststate = g.getstate()
        g = Random()
        g.setstate(laststate)
        g.jumpahead(delta)
        result.append(g)
    return result

gens = create_generators(10, 1000000)

That creates 10 distinct generators, which can be passed out to 10 distinct
threads.  The generators don't share state so can be called safely in
parallel.  So long as no thread calls its g.random() more than a million
times (the second argument to create_generators), the sequences seen by
each thread will not overlap.

The period of the underlying Wichmann-Hill generator is 6,953,607,871,644,
and that limits how far this technique can be pushed.

Just for fun, note that since we know the period, .jumpahead() can also be
used to "move backward in time":

>>> g = Random(42)  # arbitrary
>>> g.random()
69
0.25420336316883324
70 71
>>> g.jumpahead(6953607871644L - 1) # move *back* one
>>> g.random()
72
0.25420336316883324
73
"""
74
# XXX The docstring sucks.
75

76 77
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
78
from math import floor as _floor
79

80 81 82 83 84
__all__ = ["Random","seed","random","uniform","randint","choice",
           "randrange","shuffle","normalvariate","lognormvariate",
           "cunifvariate","expovariate","vonmisesvariate","gammavariate",
           "stdgamma","gauss","betavariate","paretovariate","weibullvariate",
           "getstate","setstate","jumpahead","whseed"]
Tim Peters's avatar
Tim Peters committed
85

86
def _verify(name, computed, expected):
87 88 89 90
    if abs(computed - expected) > 1e-7:
        raise ValueError(
            "computed value for %s deviates too much "
            "(computed %g, expected %g)" % (name, computed, expected))
91

92
NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
93
_verify('NV_MAGICCONST', NV_MAGICCONST, 1.71552776992141)
94

95
TWOPI = 2.0*_pi
96
_verify('TWOPI', TWOPI, 6.28318530718)
Tim Peters's avatar
Tim Peters committed
97

98
LOG4 = _log(4.0)
99
_verify('LOG4', LOG4, 1.38629436111989)
Tim Peters's avatar
Tim Peters committed
100

101
SG_MAGICCONST = 1.0 + _log(4.5)
102
_verify('SG_MAGICCONST', SG_MAGICCONST, 2.50407739677627)
103

104
del _verify
105

106 107
# Translated by Guido van Rossum from C source provided by
# Adrian Baddeley.
108

109
class Random:
110 111 112 113 114 115 116 117 118 119 120
    """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().
121

122
    """
123

124
    VERSION = 1     # used by getstate/setstate
125

126 127
    def __init__(self, x=None):
        """Initialize an instance.
128

129 130
        Optional argument x controls seeding, as for Random.seed().
        """
131

132
        self.seed(x)
133

134 135
## -------------------- core generator -------------------

136
    # Specific to Wichmann-Hill generator.  Subclasses wishing to use a
137
    # different core generator should override the seed(), random(),
138
    # getstate(), setstate() and jumpahead() methods.
139

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

143 144
        None or no argument seeds from current time.

Tim Peters's avatar
Tim Peters committed
145
        If a is not None or an int or long, hash(a) is used instead.
146 147 148 149 150

        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).
151 152
        """

153
        if a is None:
154 155
            # Initialize from current time
            import time
156 157 158 159 160 161 162 163 164
            a = long(time.time() * 256)

        if type(a) not in (type(3), type(3L)):
            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
165

166 167
        self.gauss_next = None

168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
    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

199 200 201 202 203 204 205 206 207 208 209 210 211 212
    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))

213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
    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

235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
    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) == type(0):
            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)

256 257
        self.gauss_next = None

258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
    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)

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

284
## -------------------- pickle support  -------------------
285

286 287
    def __getstate__(self): # for pickle
        return self.getstate()
288

289 290 291 292
    def __setstate__(self, state):  # for pickle
        self.setstate(state)

## -------------------- integer methods  -------------------
293 294 295 296 297 298 299 300 301 302

    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
303
        # common case while still doing adequate error checking.
304 305 306 307 308 309 310
        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()"
311 312

        # stop argument supplied.
313 314 315
        istop = int(stop)
        if istop != stop:
            raise ValueError, "non-integer stop for randrange()"
316 317 318 319 320 321 322 323 324 325 326
        if step == 1 and istart < istop:
            try:
                return istart + int(self.random()*(istop - istart))
            except OverflowError:
                # This can happen if istop-istart > sys.maxint + 1, and
                # multiplying by random() doesn't reduce it to something
                # <= sys.maxint.  We know that the overall result fits
                # in an int, and can still do it correctly via math.floor().
                # But that adds another function call, so for speed we
                # avoided that whenever possible.
                return int(istart + _floor(self.random()*(istop - istart)))
327 328
        if step == 1:
            raise ValueError, "empty range for randrange()"
329 330

        # Non-unit step argument supplied.
331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
        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):
346
        """Return random integer in range [a, b], including both end points.
347 348 349 350
        """

        return self.randrange(a, b+1)

351 352
## -------------------- sequence methods  -------------------

353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371
    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):
372
            # pick an element in x[:i+1] with which to exchange x[i]
373 374 375
            j = int(random() * (i+1))
            x[i], x[j] = x[j], x[i]

376 377 378
## -------------------- real-valued distributions  -------------------

## -------------------- uniform distribution -------------------
379 380 381 382

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

384
## -------------------- normal distribution --------------------
385

386
    def normalvariate(self, mu, sigma):
387 388 389
        """Normal distribution.

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

391
        """
392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407
        # mu = mean, sigma = standard deviation

        # Uses Kinderman and Monahan method. Reference: Kinderman,
        # A.J. and Monahan, J.F., "Computer generation of random
        # variables using the ratio of uniform deviates", ACM Trans
        # Math Software, 3, (1977), pp257-260.

        random = self.random
        while 1:
            u1 = random()
            u2 = random()
            z = NV_MAGICCONST*(u1-0.5)/u2
            zz = z*z/4.0
            if zz <= -_log(u2):
                break
        return mu + z*sigma
408

409
## -------------------- lognormal distribution --------------------
410

411
    def lognormvariate(self, mu, sigma):
412 413 414 415 416
        """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.
417

418
        """
419
        return _exp(self.normalvariate(mu, sigma))
420

421
## -------------------- circular uniform --------------------
422

423
    def cunifvariate(self, mean, arc):
424 425 426 427 428 429 430 431 432
        """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
433

434
        """
435 436
        # mean: mean angle (in radians between 0 and pi)
        # arc:  range of distribution (in radians between 0 and pi)
437 438 439 440
        import warnings
        warnings.warn("The cunifvariate function is deprecated; Use (mean "
                      "+ arc * (Random.random() - 0.5)) % Math.pi instead",
                      DeprecationWarning)
441

442
        return (mean + arc * (self.random() - 0.5)) % _pi
443

444
## -------------------- exponential distribution --------------------
445

446
    def expovariate(self, lambd):
447 448 449 450 451
        """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.
452

453
        """
454 455
        # lambd: rate lambd = 1/mean
        # ('lambda' is a Python reserved word)
456

457
        random = self.random
Tim Peters's avatar
Tim Peters committed
458
        u = random()
459 460 461
        while u <= 1e-7:
            u = random()
        return -_log(u)/lambd
462

463
## -------------------- von Mises distribution --------------------
464

465
    def vonmisesvariate(self, mu, kappa):
466
        """Circular data distribution.
467

468 469 470 471
        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.
472

473
        """
474 475 476
        # mu:    mean angle (in radians between 0 and 2*pi)
        # kappa: concentration parameter kappa (>= 0)
        # if kappa = 0 generate uniform random angle
477

478 479 480
        # Based upon an algorithm published in: Fisher, N.I.,
        # "Statistical Analysis of Circular Data", Cambridge
        # University Press, 1993.
481

482 483
        # Thanks to Magnus Kessler for a correction to the
        # implementation of step 4.
484

485 486 487
        random = self.random
        if kappa <= 1e-6:
            return TWOPI * random()
488

489 490 491
        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)
492

493 494
        while 1:
            u1 = random()
495

496 497 498
            z = _cos(_pi * u1)
            f = (1.0 + r * z)/(r + z)
            c = kappa * (r - f)
499

500
            u2 = random()
501

502 503
            if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
                break
504

505 506 507 508 509
        u3 = random()
        if u3 > 0.5:
            theta = (mu % TWOPI) + _acos(f)
        else:
            theta = (mu % TWOPI) - _acos(f)
510

511
        return theta
512

513
## -------------------- gamma distribution --------------------
514

515
    def gammavariate(self, alpha, beta):
516 517 518 519 520
        """Gamma distribution.  Not the gamma function!

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

        """
Tim Peters's avatar
Tim Peters committed
521

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

524 525 526 527
        # 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
528

529 530 531 532 533 534 535
        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

536 537 538
            ainv = _sqrt(2.0 * alpha - 1.0)
            bbb = alpha - LOG4
            ccc = alpha + ainv
Tim Peters's avatar
Tim Peters committed
539

540 541 542 543 544 545 546 547
            while 1:
                u1 = random()
                u2 = random()
                v = _log(u1/(1.0-u1))/ainv
                x = alpha*_exp(v)
                z = u1*u1*u2
                r = bbb+ccc*v-x
                if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
548
                    return x * beta
549 550 551

        elif alpha == 1.0:
            # expovariate(1)
Tim Peters's avatar
Tim Peters committed
552
            u = random()
553 554
            while u <= 1e-7:
                u = random()
555
            return -_log(u) * beta
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573

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

            # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

            while 1:
                u = random()
                b = (_e + alpha)/_e
                p = b*u
                if p <= 1.0:
                    x = 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
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594
            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)

595

596

597
## -------------------- Gauss (faster alternative) --------------------
598

599
    def gauss(self, mu, sigma):
600 601 602 603 604 605
        """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.
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

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

638
## -------------------- beta --------------------
639 640 641 642 643 644 645 646 647 648 649 650
## 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.
651

652
    def betavariate(self, alpha, beta):
653 654 655 656
        """Beta distribution.

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

658
        """
659

660 661 662 663 664 665 666
        # 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.))
667

668
## -------------------- Pareto --------------------
669

670
    def paretovariate(self, alpha):
671
        """Pareto distribution.  alpha is the shape parameter."""
672
        # Jain, pg. 495
673

674 675
        u = self.random()
        return 1.0 / pow(u, 1.0/alpha)
676

677
## -------------------- Weibull --------------------
678

679
    def weibullvariate(self, alpha, beta):
680 681 682
        """Weibull distribution.

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

684
        """
685
        # Jain, pg. 499; bug fix courtesy Bill Arms
686

687 688
        u = self.random()
        return alpha * pow(-_log(u), 1.0/beta)
689

690
## -------------------- test program --------------------
691

692
def _test_generator(n, funccall):
Tim Peters's avatar
Tim Peters committed
693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709
    import time
    print n, 'times', funccall
    code = compile(funccall, funccall, 'eval')
    sum = 0.0
    sqsum = 0.0
    smallest = 1e10
    largest = -1e10
    t0 = time.time()
    for i in range(n):
        x = eval(code)
        sum = sum + x
        sqsum = sqsum + x*x
        smallest = min(x, smallest)
        largest = max(x, largest)
    t1 = time.time()
    print round(t1-t0, 3), 'sec,',
    avg = sum/n
710
    stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters's avatar
Tim Peters committed
711 712
    print 'avg %g, stddev %g, min %g, max %g' % \
              (avg, stddev, smallest, largest)
713

714
def _test(N=20000):
715 716 717 718 719 720 721 722 723 724
    print 'TWOPI         =', TWOPI
    print 'LOG4          =', LOG4
    print 'NV_MAGICCONST =', NV_MAGICCONST
    print 'SG_MAGICCONST =', SG_MAGICCONST
    _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, 'expovariate(1.0)')
    _test_generator(N, 'vonmisesvariate(0.0, 1.0)')
725 726
    _test_generator(N, 'gammavariate(0.01, 1.0)')
    _test_generator(N, 'gammavariate(0.1, 1.0)')
Tim Peters's avatar
Tim Peters committed
727
    _test_generator(N, 'gammavariate(0.1, 2.0)')
728 729 730 731 732 733 734 735 736 737 738
    _test_generator(N, 'gammavariate(0.5, 1.0)')
    _test_generator(N, 'gammavariate(0.9, 1.0)')
    _test_generator(N, 'gammavariate(1.0, 1.0)')
    _test_generator(N, 'gammavariate(2.0, 1.0)')
    _test_generator(N, 'gammavariate(20.0, 1.0)')
    _test_generator(N, 'gammavariate(200.0, 1.0)')
    _test_generator(N, 'gauss(0.0, 1.0)')
    _test_generator(N, 'betavariate(3.0, 3.0)')
    _test_generator(N, 'paretovariate(1.0)')
    _test_generator(N, 'weibullvariate(1.0, 1.0)')

739 740 741 742 743 744 745 746 747 748 749 750
    # Test jumpahead.
    s = getstate()
    jumpahead(N)
    r1 = random()
    # now do it the slow way
    setstate(s)
    for i in range(N):
        random()
    r2 = random()
    if r1 != r2:
        raise ValueError("jumpahead test failed " + `(N, r1, r2)`)

751 752 753 754 755
# Create one instance, seeded from current time, and export its methods
# as module-level functions.  The functions are not threadsafe, and state
# is shared 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.
756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776
_inst = Random()
seed = _inst.seed
random = _inst.random
uniform = _inst.uniform
randint = _inst.randint
choice = _inst.choice
randrange = _inst.randrange
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
777
jumpahead = _inst.jumpahead
778
whseed = _inst.whseed
779

780
if __name__ == '__main__':
781
    _test()