random.py revision 2f726e9093381572b21edbfc42659ea89fbdf686
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 44d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import log as _log, exp as _exp, pi as _pi, e as _e 45d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 469146f27b7799dab231083f194a14c6157b57549fTim Petersfrom math import floor as _floor 47d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 48f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger__all__ = ["Random","seed","random","uniform","randint","choice","sample", 490de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "randrange","shuffle","normalvariate","lognormvariate", 50f8a52d38ad784b34a60720cb148180d6eb6de373Raymond Hettinger "expovariate","vonmisesvariate","gammavariate", 51f8a52d38ad784b34a60720cb148180d6eb6de373Raymond Hettinger "gauss","betavariate","paretovariate","weibullvariate", 522f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 532f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger "Random"] 54ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 55d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 56d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 57d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 58d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 592f726e9093381572b21edbfc42659ea89fbdf686Raymond HettingerBPF = 53 # Number of bits in a float 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 61d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 6240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# Adrian Baddeley. Adapted by Raymond Hettinger for use with 6340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# the Mersenne Twister core generator. 6433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 65145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerimport _random 6640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerclass Random(_random.Random): 68c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Random number generator base class used by bound module functions. 69c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 70c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Used to instantiate instances of Random to get generators that don't 71c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger share state. Especially useful for multi-threaded programs, creating 72c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger a different instance of Random for each thread, and using the jumpahead() 73c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger method to ensure that the generated sequences seen by each thread don't 74c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger overlap. 75c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 76c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Class Random can also be subclassed if you want to use a different basic 77c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger generator of your own devising: in that case, override the following 78c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger methods: random(), seed(), getstate(), setstate() and jumpahead(). 792f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Optionally, implement a getrandombits() method so that randrange() 802f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger can cover arbitrarily large ranges. 81ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 82c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 8333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 8440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 2 # used by getstate/setstate 8533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 86d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 87d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 8833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 89d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 90d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 9133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 92d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 9340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 94ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 950de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 960de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 97d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 980de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. 990de88fc4b108751b86443852b6741680d704168fTim Peters 100bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 101d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 102d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1033081d59f920229b26293c7a3ee3f1a9da0da53e9Raymond Hettinger if a is None: 1043081d59f920229b26293c7a3ee3f1a9da0da53e9Raymond Hettinger import time 1053081d59f920229b26293c7a3ee3f1a9da0da53e9Raymond Hettinger a = long(time.time() * 256) # use fractional seconds 106145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).seed(a) 10746c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 10846c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 109d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 110d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 111145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger return self.VERSION, super(Random, self).getstate(), self.gauss_next 112d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 113d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 114d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 115d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 11640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 2: 11740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, internalstate, self.gauss_next = state 118145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).setstate(internalstate) 119d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 120d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 121d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 123d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 124cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 125cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 126d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 127cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 128d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 129cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 130cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 131d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 132cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 133cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 134cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 1355f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger def __reduce__(self): 1365f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger return self.__class__, (), self.getstate() 1375f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger 138cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 139d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1402f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger def randrange(self, start, stop=None, step=1, int=int, default=None, 1412f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger maxwidth=1L<<BPF): 142d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 143d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 144d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 145d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 1462f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Do not supply the 'int', 'default', and 'maxwidth' arguments. 147d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 148d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 149d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 1509146f27b7799dab231083f194a14c6157b57549fTim Peters # common case while still doing adequate error checking. 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 152d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 153d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 154d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 155d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 1562f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if istart >= maxwidth: 1572f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return self._randbelow(istart) 158d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 159d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1609146f27b7799dab231083f194a14c6157b57549fTim Peters 1619146f27b7799dab231083f194a14c6157b57549fTim Peters # stop argument supplied. 162d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 163d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 164d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 1652f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger width = istop - istart 1662f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if step == 1 and width > 0: 16776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # Note that 1682f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # int(istart + self.random()*width) 16976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # instead would be incorrect. For example, consider istart 17076ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # = -2 and istop = 0. Then the guts would be in 17176ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # -2.0 to 0.0 exclusive on both ends (ignoring that random() 17276ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # might return 0.0), and because int() truncates toward 0, the 17376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # final result would be -1 or 0 (instead of -2 or -1). 1742f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # istart + int(self.random()*width) 17576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # would also be incorrect, for a subtler reason: the RHS 17676ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # can return a long, and then randrange() would also return 17776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # a long, but we're supposed to return an int (for backward 17876ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # compatibility). 1792f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 1802f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if width >= maxwidth: 1812f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(istart + self._randbelow(width)) 1822f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(istart + int(self.random()*width)) 183d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 1842f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 1859146f27b7799dab231083f194a14c6157b57549fTim Peters 1869146f27b7799dab231083f194a14c6157b57549fTim Peters # Non-unit step argument supplied. 187d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 188d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 189d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 190d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 1912f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger n = (width + istep - 1) / istep 192d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 1932f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger n = (width + istep + 1) / istep 194d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 196d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 197d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1992f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2002f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if n >= maxwidth: 2012f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return istart + self._randbelow(n) 202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 205cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 207d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 208d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 209d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 2102f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 2112f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 2122f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger """Return a random int in the range [0,n) 2132f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2142f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Handles the case where n has more bits than returned 2152f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger by a single call to the underlying generator. 2162f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger """ 2172f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2182f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger try: 2192f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger getrandbits = self.getrandbits 2202f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger except AttributeError: 2212f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger pass 2222f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger else: 2232f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # Only call self.getrandbits if the original random() builtin method 2242f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # has not been overridden or if a new getrandbits() was supplied. 2252f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # This assures that the two methods correspond. 2262f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 2272f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 2282f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger r = getrandbits(k) 2292f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger while r >= n: 2302f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger r = getrandbits(k) 2312f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return r 2322f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if n >= _maxwidth: 2332f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger _warn("Underlying random() generator does not supply \n" 2342f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger "enough bits to choose from a population range this large") 2352f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(self.random() * n) 2362f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 237cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 238cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 239d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 240d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 241d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return seq[int(self.random() * len(seq))] 242d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 243d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 244d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 245d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 246d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 247d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 248d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 249d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 250d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 251d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 252d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 253d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 254d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 255d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 256d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 257d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters for i in xrange(len(x)-1, 0, -1): 258cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 259d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 260d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 261d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 262fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger def sample(self, population, k): 263f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """Chooses k unique random elements from a population sequence. 264f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 265c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Returns a new list containing elements from the population while 266c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger leaving the original population unchanged. The resulting list is 267c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger in selection order so that all sub-slices will also be valid random 268c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger samples. This allows raffle winners (the sample) to be partitioned 269c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger into grand prize and second place winners (the subslices). 270f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 271c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Members of the population need not be hashable or unique. If the 272c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger population contains repeats, then each occurrence is a possible 273c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger selection in the sample. 274f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 275c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger To choose a sample in a range of integers, use xrange as an argument. 276c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger This is especially fast and space efficient for sampling from a 277c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger large population: sample(xrange(10000000), 60) 278f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """ 279f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 280c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger # Sampling without replacement entails tracking either potential 2818b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections (the pool) in a list or previous selections in a 2828b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary. 283c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 2848b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # When the number of selections is small compared to the population, 2858b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # then tracking selections is efficient, requiring only a small 2868b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary and an occasional reselection. For a larger number of 2878b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections, the pool tracking method is preferred since the list takes 2888b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # less space than the dictionary and it doesn't suffer from frequent 2898b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # reselections. 290c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 291f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger n = len(population) 292f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if not 0 <= k <= n: 293f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger raise ValueError, "sample larger than population" 2948b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger random = self.random 295fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger _int = int 296c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result = [None] * k 297f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if n < 6 * k: # if n len list takes less space than a k len dict 298311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger pool = list(population) 299311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger for i in xrange(k): # invariant: non-selected at [0,n-i) 300fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * (n-i)) 301311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger result[i] = pool[j] 3028b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger pool[j] = pool[n-i-1] # move non-selected item into vacancy 303c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger else: 30466d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger try: 30566d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger n > 0 and (population[0], population[n//2], population[n-1]) 30666d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger except (TypeError, KeyError): # handle sets and dictionaries 30766d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger population = tuple(population) 308311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger selected = {} 309c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger for i in xrange(k): 310fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 311311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while j in selected: 312fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 313c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result[i] = selected[j] = population[j] 314311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger return result 315f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 316cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 317cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 318cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 319d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 320d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 321d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 322d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 323ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 324cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 325ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 326d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 327c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Normal distribution. 328c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 329c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. 330ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 331c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 333d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 334d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 335d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 336d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 337d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 338d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 339d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 340311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 341d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 34273ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 343d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 346d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 348ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 349cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 350ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 351d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 352c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Log normal distribution. 353c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 354c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger If you take the natural logarithm of this distribution, you'll get a 355c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger normal distribution with mean mu and standard deviation sigma. 356c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu can have any value, and sigma must be greater than zero. 357ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 358c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 359d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 360ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 361cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 362ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 364c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Exponential distribution. 365c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 366c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger lambd is 1.0 divided by the desired mean. (The parameter would be 367c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger called "lambda", but that is a reserved word in Python.) Returned 368c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger values range from 0 to positive infinity. 369ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 370c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 371d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 372d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 373ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 374d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3750c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 376d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 377d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 378d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 379ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 380cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 381ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 382d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 383c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular data distribution. 384ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 385c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean angle, expressed in radians between 0 and 2*pi, and 386c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger kappa is the concentration parameter, which must be greater than or 387c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger equal to zero. If kappa is equal to zero, this distribution reduces 388c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger to a uniform random angle over the range 0 to 2*pi. 389ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 390c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 391d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 392d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 393d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 3945810297052003f28788f6790ac799fe8e5373494Guido van Rossum 395d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 396d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 397d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 3985810297052003f28788f6790ac799fe8e5373494Guido van Rossum 399d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 400d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 4015810297052003f28788f6790ac799fe8e5373494Guido van Rossum 402d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 403d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 404d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 405ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 406d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 407d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 408d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 409ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 410311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 411d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 412ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 413d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 416ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 417d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 418ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 421ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 424d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 427ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 428d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 429ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 430cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 431ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 433c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gamma distribution. Not the gamma function! 434c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 435c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > 0 and beta > 0. 436c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 437c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 4388ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 439b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 4408ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 441570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # Warning: a few older sources define the gamma distribution in terms 442570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # of alpha > -1.0 443570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum if alpha <= 0.0 or beta <= 0.0: 444570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 4458ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 446d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 447d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 448d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 449d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 450d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 451d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 452d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 453ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 454ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 455ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 4568ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 457311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 45973ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger if not 1e-7 < u1 < .9999999: 46073ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger continue 46173ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 463d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 464d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 465d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 466d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 467b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 468d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 469d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4710c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 472d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 473d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 474b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return -_log(u) * beta 475d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 478d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 479d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 480311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 483d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 485d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 486d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 489d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 492d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 493b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 494b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 495cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 49695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 498c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gaussian distribution. 499c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 500c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. This is 501c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger slightly faster than the normalvariate() function. 502c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 503c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Not thread-safe without a lock around calls. 504ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 505c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 506d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 507d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 508d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 509d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 510d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 511d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 512d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 513d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 514d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 515d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 516d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 517d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 524d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 53595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 536cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 53785e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 53885e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 53985e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 54085e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54185e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 54285e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 54385e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54485e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 54585e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 54685e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 54785e2e4742d0a1accecd02058a7907df36308297eTim Peters## 54885e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 54995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 550d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 551c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Beta distribution. 552c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 553c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > -1 and beta} > -1. 554c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Returned values range between 0 and 1. 555ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 556c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 557ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 55885e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 55985e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 56085e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 56185e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 56285e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 56385e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 56485e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 56595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 566cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 567cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 568d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 569c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Pareto distribution. alpha is the shape parameter.""" 570d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 571cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 57273ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 573d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 574cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 575cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 576cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 577d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 578c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Weibull distribution. 579c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 580c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger alpha is the scale parameter and beta is the shape parameter. 581ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 582c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 583d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 584cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 58573ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 586d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5876c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 58840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger## -------------------- Wichmann-Hill ------------------- 58940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettingerclass WichmannHill(Random): 59140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 1 # used by getstate/setstate 59340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def seed(self, a=None): 59540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Initialize internal state from hashable object. 59640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. 59840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is not None or an int or long, hash(a) is used instead. 60040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is an int or long, a is used directly. Distinct values between 60240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 0 and 27814431486575L inclusive are guaranteed to yield distinct 60340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger internal states (this guarantee is specific to the default 60440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Wichmann-Hill generator). 60540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 60640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 60840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 60940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 61040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = long(time.time() * 256) 61140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not isinstance(a, (int, long)): 61340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 61440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 30268) 61640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 30306) 61740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 30322) 61840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = int(x)+1, int(y)+1, int(z)+1 61940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 62140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def random(self): 62340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 62440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichman-Hill random number generator. 62640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 62740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichmann, B. A. & Hill, I. D. (1982) 62840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Algorithm AS 183: 62940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # An efficient and portable pseudo-random number generator 63040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 31 (1982) 188-190 63140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 63240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # see also: 63340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Correction to Algorithm AS 183 63440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 33 (1984) 123 63540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 63640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # McLeod, A. I. (1985) 63740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # A remark on Algorithm AS 183 63840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 34 (1985),198-200 63940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # This part is thread-unsafe: 64140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # BEGIN CRITICAL SECTION 64240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 64340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (171 * x) % 30269 64440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (172 * y) % 30307 64540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (170 * z) % 30323 64640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 64740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # END CRITICAL SECTION 64840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Note: on a platform using IEEE-754 double arithmetic, this can 65040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # never return 0.0 (asserted by Tim; proof too long for a comment). 65140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 65240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 65340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def getstate(self): 65440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Return internal state; can be passed to setstate() later.""" 65540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return self.VERSION, self._seed, self.gauss_next 65640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 65740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def setstate(self, state): 65840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Restore internal state from object returned by getstate().""" 65940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version = state[0] 66040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 1: 66140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, self._seed, self.gauss_next = state 66240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger else: 66340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("state with version %s passed to " 66440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "Random.setstate() of version %s" % 66540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger (version, self.VERSION)) 66640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def jumpahead(self, n): 66840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Act as if n calls to random() were made, but quickly. 66940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger n is an int, greater than or equal to 0. 67140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Example use: If you have 2 threads and know that each will 67340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger consume no more than a million random numbers, create two Random 67440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger objects r1 and r2, then do 67540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.setstate(r1.getstate()) 67640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.jumpahead(1000000) 67740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Then r1 and r2 will use guaranteed-disjoint segments of the full 67840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger period. 67940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 68040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not n >= 0: 68240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("n must be >= 0") 68340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 68440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = int(x * pow(171, n, 30269)) % 30269 68540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = int(y * pow(172, n, 30307)) % 30307 68640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = int(z * pow(170, n, 30323)) % 30323 68740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 68840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def __whseed(self, x=0, y=0, z=0): 69040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Set the Wichmann-Hill seed from (x, y, z). 69140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 69240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger These must be integers in the range [0, 256). 69340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 69440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 69540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not type(x) == type(y) == type(z) == int: 69640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise TypeError('seeds must be integers') 69740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 69840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError('seeds must be in range(0, 256)') 69940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if 0 == x == y == z: 70040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 70140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 70240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = long(time.time() * 256) 70340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = int((t&0xffffff) ^ (t>>24)) 70440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, x = divmod(t, 256) 70540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, y = divmod(t, 256) 70640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, z = divmod(t, 256) 70740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Zero is a poor seed, so substitute 1 70840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = (x or 1, y or 1, z or 1) 70940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 71140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def whseed(self, a=None): 71340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Seed from hashable object's hash code. 71440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. It is not guaranteed 71640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger that objects with distinct hash codes lead to distinct internal 71740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger states. 71840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger This is obsolete, provided for compatibility with the seed routine 72040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger used prior to Python 2.1. Use the .seed() method instead. 72140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 72240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 72340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 72440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed() 72540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return 72640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 72740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 256) 72840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 256) 72940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 256) 73040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (x + a) % 256 or 1 73140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (y + a) % 256 or 1 73240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (z + a) % 256 or 1 73340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed(x, y, z) 73440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 735cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 736ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 73762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettingerdef _test_generator(n, func, args): 7380c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 73962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger print n, 'times', func.__name__ 740b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total = 0.0 7410c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 7420c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 7430c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 7440c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 7450c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 74662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger x = func(*args) 747b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total += x 7480c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 7490c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 7500c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 7510c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 7520c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 753b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger avg = total/n 754d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 7550c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 7560c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 757ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 758f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 759f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingerdef _test(N=2000): 76062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, random, ()) 76162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, normalvariate, (0.0, 1.0)) 76262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, lognormvariate, (0.0, 1.0)) 76362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, vonmisesvariate, (0.0, 1.0)) 76462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.01, 1.0)) 76562297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.1, 1.0)) 76662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.1, 2.0)) 76762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.5, 1.0)) 76862297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.9, 1.0)) 76962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (1.0, 1.0)) 77062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (2.0, 1.0)) 77162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (20.0, 1.0)) 77262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (200.0, 1.0)) 77362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gauss, (0.0, 1.0)) 77462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, betavariate, (3.0, 3.0)) 775cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 776715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 77740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# as module-level functions. The functions share state across all uses 77840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger#(both in the user's code and in the Python libraries), but that's fine 77940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# for most programs and is easier for the casual user than making them 78040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# instantiate their own Random() instance. 78140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 782d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 783d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 784d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 785d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 786d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 787d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 788d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 789f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingersample = _inst.sample 790d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 791d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 792d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 793d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 794d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 795d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 796d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 797d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 798d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 799d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 800d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 801d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 802d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 8032f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingergetrandbits = _inst.getrandbits 804d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 805ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 806d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 807