random.py revision 91e27c253c8bb8b6ae8521f1dbb76de7c66ad8cf
1e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum"""Random variable generators.
2e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum
3d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    integers
4d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    --------
5d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters           uniform within range
6d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
7d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    sequences
8d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    ---------
9d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters           pick random element
10f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger           pick random sample
11d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters           generate random permutation
12d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
13e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum    distributions on the real line:
14e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum    ------------------------------
15d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters           uniform
16e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           normal (Gaussian)
17e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           lognormal
18e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           negative exponential
19e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           gamma
20e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           beta
2140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger           pareto
2240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger           Weibull
23e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum
24e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum    distributions on the circle (angles 0 to 2pi)
25e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum    ---------------------------------------------
26e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           circular uniform
27e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum           von Mises
28e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum
2940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond HettingerGeneral notes on the underlying Mersenne Twister core generator:
3040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
3140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* The period is 2**19937-1.
3240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* It is one of the most extensively tested generators in existence
3340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* Without a direct way to compute N steps forward, the
3440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger  semantics of jumpahead(n) are weakened to simply jump
3540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger  to another distant state and rely on the large period
3640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger  to avoid overlapping sequences.
3740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* The random() method is implemented in C, executes in
3840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger  a single Python step, and is, therefore, threadsafe.
3940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
40e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum"""
41d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum
422f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingerfrom warnings import warn as _warn
432f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingerfrom types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
4491e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettingerfrom math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
45d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
46c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettingerfrom os import urandom as _urandom
47c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettingerfrom binascii import hexlify as _hexlify
48d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
49f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger__all__ = ["Random","seed","random","uniform","randint","choice","sample",
500de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro           "randrange","shuffle","normalvariate","lognormvariate",
51f8a52d38ad784b34a60720cb148180d6eb6de373Raymond Hettinger           "expovariate","vonmisesvariate","gammavariate",
52f8a52d38ad784b34a60720cb148180d6eb6de373Raymond Hettinger           "gauss","betavariate","paretovariate","weibullvariate",
53356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger           "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
5423f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger           "SystemRandom"]
55ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
56d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
57d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi
58d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0)
59d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5)
602f726e9093381572b21edbfc42659ea89fbdf686Raymond HettingerBPF = 53        # Number of bits in a float
617c2a85b2d44851c2442ade579b760f86447bf848Tim PetersRECIP_BPF = 2**-BPF
6233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum
63356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
64d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by
6540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# Adrian Baddeley.  Adapted by Raymond Hettinger for use with
663fa19d7ff89be87139e2864fb9186b424d180a58Raymond Hettinger# the Mersenne Twister  and os.urandom() core generators.
6733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum
68145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerimport _random
6940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
70145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerclass Random(_random.Random):
71c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    """Random number generator base class used by bound module functions.
72c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
73c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    Used to instantiate instances of Random to get generators that don't
74c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    share state.  Especially useful for multi-threaded programs, creating
75c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    a different instance of Random for each thread, and using the jumpahead()
76c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    method to ensure that the generated sequences seen by each thread don't
77c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    overlap.
78c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
79c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    Class Random can also be subclassed if you want to use a different basic
80c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    generator of your own devising: in that case, override the following
81c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    methods:  random(), seed(), getstate(), setstate() and jumpahead().
822f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger    Optionally, implement a getrandombits() method so that randrange()
832f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger    can cover arbitrarily large ranges.
84ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
85c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger    """
8633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum
8740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    VERSION = 2     # used by getstate/setstate
8833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum
89d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def __init__(self, x=None):
90d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """Initialize an instance.
9133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum
92d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        Optional argument x controls seeding, as for Random.seed().
93d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """
9433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum
95d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        self.seed(x)
9640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self.gauss_next = None
97ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
980de88fc4b108751b86443852b6741680d704168fTim Peters    def seed(self, a=None):
990de88fc4b108751b86443852b6741680d704168fTim Peters        """Initialize internal state from hashable object.
100d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
10123f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        None or no argument seeds from current time or from an operating
10223f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        system specific randomness source if available.
1030de88fc4b108751b86443852b6741680d704168fTim Peters
104bcd725fc456faca13f4598f87c0517f917711cdaTim Peters        If a is not None or an int or long, hash(a) is used instead.
105d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """
106d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
1073081d59f920229b26293c7a3ee3f1a9da0da53e9Raymond Hettinger        if a is None:
108c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger            try:
109c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger                a = long(_hexlify(_urandom(16)), 16)
110c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger            except NotImplementedError:
111356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger                import time
112356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger                a = long(time.time() * 256) # use fractional seconds
113356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
114145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger        super(Random, self).seed(a)
11546c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters        self.gauss_next = None
11646c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters
117d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def getstate(self):
118d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """Return internal state; can be passed to setstate() later."""
119145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger        return self.VERSION, super(Random, self).getstate(), self.gauss_next
120d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
121d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def setstate(self, state):
122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """Restore internal state from object returned by getstate()."""
123d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        version = state[0]
12440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if version == 2:
12540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            version, internalstate, self.gauss_next = state
126145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger            super(Random, self).setstate(internalstate)
127d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        else:
128d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError("state with version %s passed to "
129d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                             "Random.setstate() of version %s" %
130d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                             (version, self.VERSION))
131d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
132cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when
133cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator.
134d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
135cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support  -------------------
136d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
137cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters    def __getstate__(self): # for pickle
138cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters        return self.getstate()
139d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
140cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters    def __setstate__(self, state):  # for pickle
141cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters        self.setstate(state)
142cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters
1435f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger    def __reduce__(self):
1445f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger        return self.__class__, (), self.getstate()
1455f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger
146cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods  -------------------
147d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
1482f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger    def randrange(self, start, stop=None, step=1, int=int, default=None,
1492f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                  maxwidth=1L<<BPF):
150d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """Choose a random item from range(start, stop[, step]).
151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
152d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        This fixes the problem with randint() which includes the
153d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        endpoint; in Python this is usually not what you want.
1542f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        Do not supply the 'int', 'default', and 'maxwidth' arguments.
155d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """
156d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
157d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # This code is a bit messy to make it fast for the
1589146f27b7799dab231083f194a14c6157b57549fTim Peters        # common case while still doing adequate error checking.
159d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        istart = int(start)
160d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if istart != start:
161d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError, "non-integer arg 1 for randrange()"
162d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if stop is default:
163d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            if istart > 0:
1642f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                if istart >= maxwidth:
1652f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                    return self._randbelow(istart)
166d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                return int(self.random() * istart)
167d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError, "empty range for randrange()"
1689146f27b7799dab231083f194a14c6157b57549fTim Peters
1699146f27b7799dab231083f194a14c6157b57549fTim Peters        # stop argument supplied.
170d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        istop = int(stop)
171d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if istop != stop:
172d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError, "non-integer stop for randrange()"
1732f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        width = istop - istart
1742f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        if step == 1 and width > 0:
17576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # Note that
1762f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            #     int(istart + self.random()*width)
17776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # instead would be incorrect.  For example, consider istart
17876ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # = -2 and istop = 0.  Then the guts would be in
17976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # -2.0 to 0.0 exclusive on both ends (ignoring that random()
18076ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # might return 0.0), and because int() truncates toward 0, the
18176ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # final result would be -1 or 0 (instead of -2 or -1).
1822f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            #     istart + int(self.random()*width)
18376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # would also be incorrect, for a subtler reason:  the RHS
18476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # can return a long, and then randrange() would also return
18576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # a long, but we're supposed to return an int (for backward
18676ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters            # compatibility).
1872f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger
1882f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            if width >= maxwidth:
18958eb11cf62dd04ccc2c364b62fd51b4265e2e203Tim Peters                return int(istart + self._randbelow(width))
1902f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            return int(istart + int(self.random()*width))
191d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if step == 1:
1922f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)
1939146f27b7799dab231083f194a14c6157b57549fTim Peters
1949146f27b7799dab231083f194a14c6157b57549fTim Peters        # Non-unit step argument supplied.
195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        istep = int(step)
196d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if istep != step:
197d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError, "non-integer step for randrange()"
198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if istep > 0:
199ffdb8bb99c4017152a9dca70669f9d6b9831d454Raymond Hettinger            n = (width + istep - 1) // istep
200d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        elif istep < 0:
201ffdb8bb99c4017152a9dca70669f9d6b9831d454Raymond Hettinger            n = (width + istep + 1) // istep
202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        else:
203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError, "zero step for randrange()"
204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
205d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if n <= 0:
206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            raise ValueError, "empty range for randrange()"
2072f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger
2082f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        if n >= maxwidth:
2092f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            return istart + self._randbelow(n)
210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return istart + istep*int(self.random() * n)
211d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
212d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def randint(self, a, b):
213cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters        """Return random integer in range [a, b], including both end points.
214d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """
215d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
216d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return self.randrange(a, b+1)
217d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
2182f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger    def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF,
2192f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                   _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
2202f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        """Return a random int in the range [0,n)
2212f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger
2222f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        Handles the case where n has more bits than returned
2232f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        by a single call to the underlying generator.
2242f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        """
2252f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger
2262f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        try:
2272f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            getrandbits = self.getrandbits
2282f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        except AttributeError:
2292f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            pass
2302f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        else:
2312f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            # Only call self.getrandbits if the original random() builtin method
2322f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            # has not been overridden or if a new getrandbits() was supplied.
2332f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            # This assures that the two methods correspond.
2342f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
2352f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                k = int(1.00001 + _log(n-1, 2.0))   # 2**k > n-1 > 2**(k-2)
2362f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                r = getrandbits(k)
2372f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                while r >= n:
2382f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                    r = getrandbits(k)
2392f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                return r
2402f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        if n >= _maxwidth:
2412f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger            _warn("Underlying random() generator does not supply \n"
2422f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger                "enough bits to choose from a population range this large")
2432f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger        return int(self.random() * n)
2442f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger
245cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods  -------------------
246cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters
247d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def choice(self, seq):
248d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """Choose a random element from a non-empty sequence."""
2495dae505bbd59641a948c81bea981e7c44d4c2343Raymond Hettinger        return seq[int(self.random() * len(seq))]  # raises IndexError if seq is empty
250d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
251d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def shuffle(self, x, random=None, int=int):
252d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """x, random=random.random -> shuffle list x in place; return None.
253d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
254d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        Optional arg random is a 0-argument function returning a random
255d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        float in [0.0, 1.0); by default, the standard random.random.
256d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
257d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        Note that for even rather small len(x), the total number of
258d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        permutations of x is larger than the period of most random number
259d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        generators; this implies that "most" permutations of a long
260d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        sequence can never be generated.
261d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """
262d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
263d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if random is None:
264d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            random = self.random
26585c20a41dfcec04d161ad7da7260e7b94c62d228Raymond Hettinger        for i in reversed(xrange(1, len(x))):
266cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters            # pick an element in x[:i+1] with which to exchange x[i]
267d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            j = int(random() * (i+1))
268d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            x[i], x[j] = x[j], x[i]
269d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
270fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger    def sample(self, population, k):
271f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger        """Chooses k unique random elements from a population sequence.
272f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger
273c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        Returns a new list containing elements from the population while
274c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        leaving the original population unchanged.  The resulting list is
275c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        in selection order so that all sub-slices will also be valid random
276c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        samples.  This allows raffle winners (the sample) to be partitioned
277c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        into grand prize and second place winners (the subslices).
278f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger
279c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        Members of the population need not be hashable or unique.  If the
280c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        population contains repeats, then each occurrence is a possible
281c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        selection in the sample.
282f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger
283c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        To choose a sample in a range of integers, use xrange as an argument.
284c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        This is especially fast and space efficient for sampling from a
285c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        large population:   sample(xrange(10000000), 60)
286f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger        """
287f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger
288c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        # Sampling without replacement entails tracking either potential
28991e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger        # selections (the pool) in a list or previous selections in a set.
290c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger
2912b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton        # When the number of selections is small compared to the
2922b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton        # population, then tracking selections is efficient, requiring
29391e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger        # only a small set and an occasional reselection.  For
2942b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton        # a larger number of selections, the pool tracking method is
2952b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton        # preferred since the list takes less space than the
29691e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger        # set and it doesn't suffer from frequent reselections.
297c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger
298f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger        n = len(population)
299f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger        if not 0 <= k <= n:
300f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger            raise ValueError, "sample larger than population"
3018b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger        random = self.random
302fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger        _int = int
303c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        result = [None] * k
30491e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger        setsize = 21        # size of a small set minus size of an empty list
30591e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger        if k > 5:
30691e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger              setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
30791e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger        if n <= setsize:    # is an n-length list smaller than a k-length set
308311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger            pool = list(population)
309311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger            for i in xrange(k):         # invariant:  non-selected at [0,n-i)
310fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger                j = _int(random() * (n-i))
311311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger                result[i] = pool[j]
3128b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger                pool[j] = pool[n-i-1]   # move non-selected item into vacancy
313c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger        else:
31466d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger            try:
31566d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger                n > 0 and (population[0], population[n//2], population[n-1])
31691e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger            except (TypeError, KeyError):   # handle non-sequence iterables
31766d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger                population = tuple(population)
31891e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger            selected = set()
31991e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger            selected_add = selected.add
320c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger            for i in xrange(k):
321fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger                j = _int(random() * n)
322311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger                while j in selected:
323fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger                    j = _int(random() * n)
32491e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger                selected_add(j)
32591e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger                result[i] = population[j]
326311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger        return result
327f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger
328cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions  -------------------
329cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters
330cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution -------------------
331d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def uniform(self, a, b):
333d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        """Get a random number in the range [a, b)."""
334d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return a + (b-a) * self.random()
335ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
336cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution --------------------
337ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
338d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def normalvariate(self, mu, sigma):
339c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Normal distribution.
340c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
341c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        mu is the mean, and sigma is the standard deviation.
342ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
343c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # mu = mean, sigma = standard deviation
345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
346d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Uses Kinderman and Monahan method. Reference: Kinderman,
347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # A.J. and Monahan, J.F., "Computer generation of random
348d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # variables using the ratio of uniform deviates", ACM Trans
349d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Math Software, 3, (1977), pp257-260.
350d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
351d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        random = self.random
35242406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger        while 1:
353d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            u1 = random()
35473ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger            u2 = 1.0 - random()
355d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            z = NV_MAGICCONST*(u1-0.5)/u2
356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            zz = z*z/4.0
357d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            if zz <= -_log(u2):
358d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                break
359d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return mu + z*sigma
360ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
361cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution --------------------
362ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def lognormvariate(self, mu, sigma):
364c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Log normal distribution.
365c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
366c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        If you take the natural logarithm of this distribution, you'll get a
367c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        normal distribution with mean mu and standard deviation sigma.
368c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        mu can have any value, and sigma must be greater than zero.
369ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
370c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
371d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return _exp(self.normalvariate(mu, sigma))
372ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
373cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution --------------------
374ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
375d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def expovariate(self, lambd):
376c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Exponential distribution.
377c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
378c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        lambd is 1.0 divided by the desired mean.  (The parameter would be
379c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        called "lambda", but that is a reserved word in Python.)  Returned
380c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        values range from 0 to positive infinity.
381ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
382c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
383d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # lambd: rate lambd = 1/mean
384d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # ('lambda' is a Python reserved word)
385ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
386d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        random = self.random
3870c9886d589ddebf32de0ca3f027a173222ed383aTim Peters        u = random()
388d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        while u <= 1e-7:
389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            u = random()
390d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return -_log(u)/lambd
391ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
392cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution --------------------
393ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def vonmisesvariate(self, mu, kappa):
395c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Circular data distribution.
396ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
397c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        mu is the mean angle, expressed in radians between 0 and 2*pi, and
398c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        kappa is the concentration parameter, which must be greater than or
399c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        equal to zero.  If kappa is equal to zero, this distribution reduces
400c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        to a uniform random angle over the range 0 to 2*pi.
401ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
402c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
403d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # mu:    mean angle (in radians between 0 and 2*pi)
404d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # kappa: concentration parameter kappa (>= 0)
405d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # if kappa = 0 generate uniform random angle
4065810297052003f28788f6790ac799fe8e5373494Guido van Rossum
407d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Based upon an algorithm published in: Fisher, N.I.,
408d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # "Statistical Analysis of Circular Data", Cambridge
409d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # University Press, 1993.
4105810297052003f28788f6790ac799fe8e5373494Guido van Rossum
411d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Thanks to Magnus Kessler for a correction to the
412d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # implementation of step 4.
4135810297052003f28788f6790ac799fe8e5373494Guido van Rossum
414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        random = self.random
415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if kappa <= 1e-6:
416d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            return TWOPI * random()
417ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        r = (1.0 + b * b)/(2.0 * b)
421ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
42242406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger        while 1:
423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            u1 = random()
424ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            z = _cos(_pi * u1)
426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            f = (1.0 + r * z)/(r + z)
427d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            c = kappa * (r - f)
428ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
429d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            u2 = random()
430ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
43142406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger            if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                break
433ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
434d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        u3 = random()
435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if u3 > 0.5:
436d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            theta = (mu % TWOPI) + _acos(f)
437d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        else:
438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            theta = (mu % TWOPI) - _acos(f)
439ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return theta
441ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
442cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution --------------------
443ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
444d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def gammavariate(self, alpha, beta):
445c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Gamma distribution.  Not the gamma function!
446c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
447c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        Conditions on the parameters are alpha > 0 and beta > 0.
448c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
449c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
4508ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters
451b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger        # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
4528ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters
453570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum        # Warning: a few older sources define the gamma distribution in terms
454570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum        # of alpha > -1.0
455570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum        if alpha <= 0.0 or beta <= 0.0:
456570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum            raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
4578ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters
458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        random = self.random
459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if alpha > 1.0:
460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            # Uses R.C.H. Cheng, "The generation of Gamma
462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            # variables with non-integral shape parameters",
463d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            # Applied Statistics, (1977), 26, No. 1, p71-74
464d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
465ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger            ainv = _sqrt(2.0 * alpha - 1.0)
466ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger            bbb = alpha - LOG4
467ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger            ccc = alpha + ainv
4688ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters
46942406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger            while 1:
470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                u1 = random()
47173ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger                if not 1e-7 < u1 < .9999999:
47273ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger                    continue
47373ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger                u2 = 1.0 - random()
474d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                v = _log(u1/(1.0-u1))/ainv
475d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                x = alpha*_exp(v)
476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                z = u1*u1*u2
477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                r = bbb+ccc*v-x
478d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
479b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger                    return x * beta
480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        elif alpha == 1.0:
482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            # expovariate(1)
4830c9886d589ddebf32de0ca3f027a173222ed383aTim Peters            u = random()
484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            while u <= 1e-7:
485d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                u = random()
486b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger            return -_log(u) * beta
487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        else:   # alpha is between 0 and 1 (exclusive)
489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
49242406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger            while 1:
493d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                u = random()
494d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                b = (_e + alpha)/_e
495d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                p = b*u
496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                if p <= 1.0:
49742406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger                    x = p ** (1.0/alpha)
498d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                else:
499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                    x = -_log((b-p)/alpha)
500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                u1 = random()
50142406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger                if p > 1.0:
50242406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger                    if u1 <= x ** (alpha - 1.0):
50342406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger                        break
50442406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger                elif u1 <= _exp(-x):
505d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters                    break
506b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger            return x * beta
507b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger
508cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) --------------------
50995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum
510d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def gauss(self, mu, sigma):
511c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Gaussian distribution.
512c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
513c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        mu is the mean, and sigma is the standard deviation.  This is
514c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        slightly faster than the normalvariate() function.
515c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
516c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        Not thread-safe without a lock around calls.
517ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
518c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # When x and y are two variables from [0, 1), uniformly
521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # distributed, then
522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        #
523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        #    cos(2*pi*x)*sqrt(-2*log(1-y))
524d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        #    sin(2*pi*x)*sqrt(-2*log(1-y))
525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        #
526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # are two *independent* variables with normal distribution
527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # (mu = 0, sigma = 1).
528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # (Lambert Meertens)
529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # (corrected version; bug discovered by Mike Miller, fixed by LM)
530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Multithreading note: When two threads call this function
532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # simultaneously, it is possible that they will receive the
533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # same return value.  The window is very small though.  To
534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # avoid this, you have to use a lock around all calls.  (I
535d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # didn't want to slow this down in the serial case by using a
536d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # lock here.)
537d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
538d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        random = self.random
539d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        z = self.gauss_next
540d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        self.gauss_next = None
541d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        if z is None:
542d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            x2pi = random() * TWOPI
543d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            g2rad = _sqrt(-2.0 * _log(1.0 - random()))
544d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            z = _cos(x2pi) * g2rad
545d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters            self.gauss_next = _sin(x2pi) * g2rad
546d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
547d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return mu + z*sigma
54895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum
549cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta --------------------
55085e2e4742d0a1accecd02058a7907df36308297eTim Peters## See
55185e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
55285e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation:
55385e2e4742d0a1accecd02058a7907df36308297eTim Peters##
55485e2e4742d0a1accecd02058a7907df36308297eTim Peters##    def betavariate(self, alpha, beta):
55585e2e4742d0a1accecd02058a7907df36308297eTim Peters##        # Discrete Event Simulation in C, pp 87-88.
55685e2e4742d0a1accecd02058a7907df36308297eTim Peters##
55785e2e4742d0a1accecd02058a7907df36308297eTim Peters##        y = self.expovariate(alpha)
55885e2e4742d0a1accecd02058a7907df36308297eTim Peters##        z = self.expovariate(1.0/beta)
55985e2e4742d0a1accecd02058a7907df36308297eTim Peters##        return z/(y+z)
56085e2e4742d0a1accecd02058a7907df36308297eTim Peters##
56185e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way.
56295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum
563d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def betavariate(self, alpha, beta):
564c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Beta distribution.
565c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
566c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        Conditions on the parameters are alpha > -1 and beta} > -1.
567c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        Returned values range between 0 and 1.
568ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
569c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
570ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
57185e2e4742d0a1accecd02058a7907df36308297eTim Peters        # This version due to Janne Sinkkonen, and matches all the std
57285e2e4742d0a1accecd02058a7907df36308297eTim Peters        # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
57385e2e4742d0a1accecd02058a7907df36308297eTim Peters        y = self.gammavariate(alpha, 1.)
57485e2e4742d0a1accecd02058a7907df36308297eTim Peters        if y == 0:
57585e2e4742d0a1accecd02058a7907df36308297eTim Peters            return 0.0
57685e2e4742d0a1accecd02058a7907df36308297eTim Peters        else:
57785e2e4742d0a1accecd02058a7907df36308297eTim Peters            return y / (y + self.gammavariate(beta, 1.))
57895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum
579cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto --------------------
580cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum
581d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def paretovariate(self, alpha):
582c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Pareto distribution.  alpha is the shape parameter."""
583d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Jain, pg. 495
584cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum
58573ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger        u = 1.0 - self.random()
586d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return 1.0 / pow(u, 1.0/alpha)
587cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum
588cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull --------------------
589cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum
590d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    def weibullvariate(self, alpha, beta):
591c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """Weibull distribution.
592c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger
593c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        alpha is the scale parameter and beta is the shape parameter.
594ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger
595c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger        """
596d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        # Jain, pg. 499; bug fix courtesy Bill Arms
597cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum
59873ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger        u = 1.0 - self.random()
599d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters        return alpha * pow(-_log(u), 1.0/beta)
6006c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum
60140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger## -------------------- Wichmann-Hill -------------------
60240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
60340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettingerclass WichmannHill(Random):
60440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
60540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    VERSION = 1     # used by getstate/setstate
60640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
60740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def seed(self, a=None):
60840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Initialize internal state from hashable object.
60940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
61023f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        None or no argument seeds from current time or from an operating
61123f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        system specific randomness source if available.
61240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
61340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        If a is not None or an int or long, hash(a) is used instead.
61440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
61540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        If a is an int or long, a is used directly.  Distinct values between
61640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        0 and 27814431486575L inclusive are guaranteed to yield distinct
61740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        internal states (this guarantee is specific to the default
61840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        Wichmann-Hill generator).
61940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """
62040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
62140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if a is None:
622c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger            try:
623c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger                a = long(_hexlify(_urandom(16)), 16)
624c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger            except NotImplementedError:
625356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger                import time
626356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger                a = long(time.time() * 256) # use fractional seconds
62740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
62840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if not isinstance(a, (int, long)):
62940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            a = hash(a)
63040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
63140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a, x = divmod(a, 30268)
63240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a, y = divmod(a, 30306)
63340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a, z = divmod(a, 30322)
63440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self._seed = int(x)+1, int(y)+1, int(z)+1
63540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
63640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self.gauss_next = None
63740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
63840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def random(self):
63940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Get the next random number in the range [0.0, 1.0)."""
64040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
64140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # Wichman-Hill random number generator.
64240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #
64340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # Wichmann, B. A. & Hill, I. D. (1982)
64440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # Algorithm AS 183:
64540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # An efficient and portable pseudo-random number generator
64640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # Applied Statistics 31 (1982) 188-190
64740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #
64840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # see also:
64940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #        Correction to Algorithm AS 183
65040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #        Applied Statistics 33 (1984) 123
65140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #
65240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #        McLeod, A. I. (1985)
65340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #        A remark on Algorithm AS 183
65440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        #        Applied Statistics 34 (1985),198-200
65540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
65640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # This part is thread-unsafe:
65740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # BEGIN CRITICAL SECTION
65840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        x, y, z = self._seed
65940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        x = (171 * x) % 30269
66040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        y = (172 * y) % 30307
66140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        z = (170 * z) % 30323
66240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self._seed = x, y, z
66340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # END CRITICAL SECTION
66440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
66540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # Note:  on a platform using IEEE-754 double arithmetic, this can
66640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # never return 0.0 (asserted by Tim; proof too long for a comment).
66740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
66840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
66940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def getstate(self):
67040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Return internal state; can be passed to setstate() later."""
67140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        return self.VERSION, self._seed, self.gauss_next
67240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
67340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def setstate(self, state):
67440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Restore internal state from object returned by getstate()."""
67540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        version = state[0]
67640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if version == 1:
67740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            version, self._seed, self.gauss_next = state
67840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        else:
67940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            raise ValueError("state with version %s passed to "
68040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger                             "Random.setstate() of version %s" %
68140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger                             (version, self.VERSION))
68240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
68340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def jumpahead(self, n):
68440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Act as if n calls to random() were made, but quickly.
68540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
68640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        n is an int, greater than or equal to 0.
68740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
68840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        Example use:  If you have 2 threads and know that each will
68940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        consume no more than a million random numbers, create two Random
69040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        objects r1 and r2, then do
69140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            r2.setstate(r1.getstate())
69240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            r2.jumpahead(1000000)
69340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        Then r1 and r2 will use guaranteed-disjoint segments of the full
69440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        period.
69540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """
69640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
69740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if not n >= 0:
69840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            raise ValueError("n must be >= 0")
69940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        x, y, z = self._seed
70040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        x = int(x * pow(171, n, 30269)) % 30269
70140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        y = int(y * pow(172, n, 30307)) % 30307
70240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        z = int(z * pow(170, n, 30323)) % 30323
70340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self._seed = x, y, z
70440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
70540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def __whseed(self, x=0, y=0, z=0):
70640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Set the Wichmann-Hill seed from (x, y, z).
70740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
70840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        These must be integers in the range [0, 256).
70940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """
71040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
71140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if not type(x) == type(y) == type(z) == int:
71240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            raise TypeError('seeds must be integers')
71340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
71440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            raise ValueError('seeds must be in range(0, 256)')
71540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if 0 == x == y == z:
71640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            # Initialize from current time
71740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            import time
71840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            t = long(time.time() * 256)
71940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            t = int((t&0xffffff) ^ (t>>24))
72040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            t, x = divmod(t, 256)
72140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            t, y = divmod(t, 256)
72240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            t, z = divmod(t, 256)
72340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        # Zero is a poor seed, so substitute 1
72440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self._seed = (x or 1, y or 1, z or 1)
72540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
72640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self.gauss_next = None
72740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
72840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger    def whseed(self, a=None):
72940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """Seed from hashable object's hash code.
73040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
73140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        None or no argument seeds from current time.  It is not guaranteed
73240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        that objects with distinct hash codes lead to distinct internal
73340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        states.
73440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
73540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        This is obsolete, provided for compatibility with the seed routine
73640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        used prior to Python 2.1.  Use the .seed() method instead.
73740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        """
73840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
73940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        if a is None:
74040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            self.__whseed()
74140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger            return
74240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a = hash(a)
74340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a, x = divmod(a, 256)
74440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a, y = divmod(a, 256)
74540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        a, z = divmod(a, 256)
74640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        x = (x + a) % 256 or 1
74740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        y = (y + a) % 256 or 1
74840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        z = (z + a) % 256 or 1
74940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger        self.__whseed(x, y, z)
75040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
75123f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger## --------------- Operating System Random Source  ------------------
752356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
75323f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettingerclass SystemRandom(Random):
75423f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger    """Alternate random number generator using sources provided
75523f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger    by the operating system (such as /dev/urandom on Unix or
75623f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger    CryptGenRandom on Windows).
757356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
758356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger     Not available on all systems (see os.urandom() for details).
759356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    """
760356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
761356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    def random(self):
762356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        """Get the next random number in the range [0.0, 1.0)."""
7637c2a85b2d44851c2442ade579b760f86447bf848Tim Peters        return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
764356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
765356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    def getrandbits(self, k):
766356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        """getrandbits(k) -> x.  Generates a long int with k random bits."""
767356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        if k <= 0:
768356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger            raise ValueError('number of bits must be greater than zero')
769356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        if k != int(k):
770356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger            raise TypeError('number of bits should be an integer')
771356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        bytes = (k + 7) // 8                    # bits / 8 and rounded up
772356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        x = long(_hexlify(_urandom(bytes)), 16)
773356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        return x >> (bytes * 8 - k)             # trim excess bits
774356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
775356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    def _stub(self, *args, **kwds):
77623f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        "Stub method.  Not used for a system random number generator."
777356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger        return None
778356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    seed = jumpahead = _stub
779356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
780356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    def _notimplemented(self, *args, **kwds):
78123f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        "Method should not be called for a system random number generator."
78223f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger        raise NotImplementedError('System entropy source does not have state.')
783356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger    getstate = setstate = _notimplemented
784356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger
785cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program --------------------
786ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
78762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettingerdef _test_generator(n, func, args):
7880c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    import time
78962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    print n, 'times', func.__name__
790b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger    total = 0.0
7910c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    sqsum = 0.0
7920c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    smallest = 1e10
7930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    largest = -1e10
7940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    t0 = time.time()
7950c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    for i in range(n):
79662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger        x = func(*args)
797b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger        total += x
7980c9886d589ddebf32de0ca3f027a173222ed383aTim Peters        sqsum = sqsum + x*x
7990c9886d589ddebf32de0ca3f027a173222ed383aTim Peters        smallest = min(x, smallest)
8000c9886d589ddebf32de0ca3f027a173222ed383aTim Peters        largest = max(x, largest)
8010c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    t1 = time.time()
8020c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    print round(t1-t0, 3), 'sec,',
803b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger    avg = total/n
804d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    stddev = _sqrt(sqsum/n - avg*avg)
8050c9886d589ddebf32de0ca3f027a173222ed383aTim Peters    print 'avg %g, stddev %g, min %g, max %g' % \
8060c9886d589ddebf32de0ca3f027a173222ed383aTim Peters              (avg, stddev, smallest, largest)
807ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum
808f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger
809f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingerdef _test(N=2000):
81062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, random, ())
81162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, normalvariate, (0.0, 1.0))
81262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, lognormvariate, (0.0, 1.0))
81362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, vonmisesvariate, (0.0, 1.0))
81462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (0.01, 1.0))
81562297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (0.1, 1.0))
81662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (0.1, 2.0))
81762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (0.5, 1.0))
81862297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (0.9, 1.0))
81962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (1.0, 1.0))
82062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (2.0, 1.0))
82162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (20.0, 1.0))
82262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gammavariate, (200.0, 1.0))
82362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, gauss, (0.0, 1.0))
82462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger    _test_generator(N, betavariate, (3.0, 3.0))
825cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters
826715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods
82740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# as module-level functions.  The functions share state across all uses
82840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger#(both in the user's code and in the Python libraries), but that's fine
82940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# for most programs and is easier for the casual user than making them
83040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# instantiate their own Random() instance.
83140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger
832d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random()
833d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed
834d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random
835d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform
836d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint
837d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice
838d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange
839f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingersample = _inst.sample
840d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle
841d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate
842d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate
843d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate
844d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate
845d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate
846d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss
847d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate
848d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate
849d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate
850d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate
851d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate
852d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead
8532f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingergetrandbits = _inst.getrandbits
854d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters
855ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__':
856d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters    _test()
857