random.py revision f8a52d38ad784b34a60720cb148180d6eb6de373
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 42d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import log as _log, exp as _exp, pi as _pi, e as _e 43d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 449146f27b7799dab231083f194a14c6157b57549fTim Petersfrom math import floor as _floor 45d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 46f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger__all__ = ["Random","seed","random","uniform","randint","choice","sample", 470de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "randrange","shuffle","normalvariate","lognormvariate", 48f8a52d38ad784b34a60720cb148180d6eb6de373Raymond Hettinger "expovariate","vonmisesvariate","gammavariate", 49f8a52d38ad784b34a60720cb148180d6eb6de373Raymond Hettinger "gauss","betavariate","paretovariate","weibullvariate", 5040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "getstate","setstate","jumpahead"] 51ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 52d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 53d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 54d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 55d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 5633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 57d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 5840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# Adrian Baddeley. Adapted by Raymond Hettinger for use with 5940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# the Mersenne Twister core generator. 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 61145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerimport _random 6240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerclass Random(_random.Random): 64c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Random number generator base class used by bound module functions. 65c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 66c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Used to instantiate instances of Random to get generators that don't 67c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger share state. Especially useful for multi-threaded programs, creating 68c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger a different instance of Random for each thread, and using the jumpahead() 69c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger method to ensure that the generated sequences seen by each thread don't 70c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger overlap. 71c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 72c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Class Random can also be subclassed if you want to use a different basic 73c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger generator of your own devising: in that case, override the following 74c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger methods: random(), seed(), getstate(), setstate() and jumpahead(). 75ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 76c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 7733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 7840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 2 # used by getstate/setstate 7933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 80d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 81d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 8233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 83d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 84d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 8533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 86d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 8740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 88ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 890de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 900de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 91d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 920de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. 930de88fc4b108751b86443852b6741680d704168fTim Peters 94bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 95d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 96d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 97145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).seed(a) 9846c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 9946c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 100d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 101d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 102145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger return self.VERSION, super(Random, self).getstate(), self.gauss_next 103d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 104d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 105d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 106d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 10740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 2: 10840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, internalstate, self.gauss_next = state 109145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).setstate(internalstate) 110d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 111d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 112d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 113d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 114d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 115cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 116cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 117d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 118cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 119d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 120cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 121cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 123cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 124cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 125cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 1265f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger def __reduce__(self): 1275f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger return self.__class__, (), self.getstate() 1285f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger 129cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 130d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 131d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randrange(self, start, stop=None, step=1, int=int, default=None): 132d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 133d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 134d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 135d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 136d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Do not supply the 'int' and 'default' arguments. 137d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 138d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 139d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 1409146f27b7799dab231083f194a14c6157b57549fTim Peters # common case while still doing adequate error checking. 141d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 142d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 143d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 144d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 145d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 146d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 147d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1489146f27b7799dab231083f194a14c6157b57549fTim Peters 1499146f27b7799dab231083f194a14c6157b57549fTim Peters # stop argument supplied. 150d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 152d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 1539146f27b7799dab231083f194a14c6157b57549fTim Peters if step == 1 and istart < istop: 15476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # Note that 15576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # int(istart + self.random()*(istop - istart)) 15676ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # instead would be incorrect. For example, consider istart 15776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # = -2 and istop = 0. Then the guts would be in 15876ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # -2.0 to 0.0 exclusive on both ends (ignoring that random() 15976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # might return 0.0), and because int() truncates toward 0, the 16076ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # final result would be -1 or 0 (instead of -2 or -1). 16176ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # istart + int(self.random()*(istop - istart)) 16276ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # would also be incorrect, for a subtler reason: the RHS 16376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # can return a long, and then randrange() would also return 16476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # a long, but we're supposed to return an int (for backward 16576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # compatibility). 16676ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters return int(istart + int(self.random()*(istop - istart))) 167d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 168d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1699146f27b7799dab231083f194a14c6157b57549fTim Peters 1709146f27b7799dab231083f194a14c6157b57549fTim Peters # Non-unit step argument supplied. 171d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 172d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 173d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 174d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 175d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep - 1) / istep 176d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 177d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep + 1) / istep 178d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 179d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 180d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 181d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 182d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 183d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 184d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 185d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 186cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 187d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 188d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 189d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 190d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 191cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 192cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 193d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 194d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return seq[int(self.random() * len(seq))] 196d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 197d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 199d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 200d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 201d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 205d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 207d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 208d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 209d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 211d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters for i in xrange(len(x)-1, 0, -1): 212cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 213d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 214d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 215d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 216fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger def sample(self, population, k): 217f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """Chooses k unique random elements from a population sequence. 218f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 219c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Returns a new list containing elements from the population while 220c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger leaving the original population unchanged. The resulting list is 221c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger in selection order so that all sub-slices will also be valid random 222c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger samples. This allows raffle winners (the sample) to be partitioned 223c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger into grand prize and second place winners (the subslices). 224f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 225c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Members of the population need not be hashable or unique. If the 226c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger population contains repeats, then each occurrence is a possible 227c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger selection in the sample. 228f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 229c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger To choose a sample in a range of integers, use xrange as an argument. 230c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger This is especially fast and space efficient for sampling from a 231c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger large population: sample(xrange(10000000), 60) 232f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """ 233f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 234c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger # Sampling without replacement entails tracking either potential 2358b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections (the pool) in a list or previous selections in a 2368b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary. 237c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 2388b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # When the number of selections is small compared to the population, 2398b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # then tracking selections is efficient, requiring only a small 2408b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary and an occasional reselection. For a larger number of 2418b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections, the pool tracking method is preferred since the list takes 2428b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # less space than the dictionary and it doesn't suffer from frequent 2438b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # reselections. 244c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 245f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger n = len(population) 246f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if not 0 <= k <= n: 247f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger raise ValueError, "sample larger than population" 2488b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger random = self.random 249fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger _int = int 250c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result = [None] * k 251f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if n < 6 * k: # if n len list takes less space than a k len dict 252311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger pool = list(population) 253311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger for i in xrange(k): # invariant: non-selected at [0,n-i) 254fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * (n-i)) 255311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger result[i] = pool[j] 2568b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger pool[j] = pool[n-i-1] # move non-selected item into vacancy 257c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger else: 258311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger selected = {} 259c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger for i in xrange(k): 260fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 261311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while j in selected: 262fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 263c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result[i] = selected[j] = population[j] 264311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger return result 265f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 266cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 267cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 268cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 269d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 270d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 271d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 272d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 273ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 274cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 275ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 276d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 277c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Normal distribution. 278c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 279c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. 280ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 281c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 282d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 283d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 284d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 285d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 286d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 287d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 288d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 289d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 290311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 291d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 29273ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 293d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 294d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 295d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 296d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 297d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 298ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 299cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 300ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 301d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 302c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Log normal distribution. 303c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 304c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger If you take the natural logarithm of this distribution, you'll get a 305c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger normal distribution with mean mu and standard deviation sigma. 306c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu can have any value, and sigma must be greater than zero. 307ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 308c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 309d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 310ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 311cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 312ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 313d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 314c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Exponential distribution. 315c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 316c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger lambd is 1.0 divided by the desired mean. (The parameter would be 317c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger called "lambda", but that is a reserved word in Python.) Returned 318c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger values range from 0 to positive infinity. 319ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 320c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 321d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 322d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 323ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 324d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3250c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 326d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 327d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 328d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 329ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 330cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 331ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 332d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 333c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular data distribution. 334ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 335c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean angle, expressed in radians between 0 and 2*pi, and 336c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger kappa is the concentration parameter, which must be greater than or 337c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger equal to zero. If kappa is equal to zero, this distribution reduces 338c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger to a uniform random angle over the range 0 to 2*pi. 339ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 340c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 341d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 342d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 343d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 3445810297052003f28788f6790ac799fe8e5373494Guido van Rossum 345d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 346d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 3485810297052003f28788f6790ac799fe8e5373494Guido van Rossum 349d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 350d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 3515810297052003f28788f6790ac799fe8e5373494Guido van Rossum 352d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 353d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 354d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 355ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 357d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 358d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 359ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 360311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 361d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 362ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 364d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 365d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 366ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 367d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 368ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 369d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 370d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 371ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 372d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 373d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 374d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 375d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 376d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 377ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 378d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 379ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 380cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 381ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 382d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 383c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gamma distribution. Not the gamma function! 384c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 385c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > 0 and beta > 0. 386c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 387c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 3888ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 389b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 3908ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 391570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # Warning: a few older sources define the gamma distribution in terms 392570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # of alpha > -1.0 393570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum if alpha <= 0.0 or beta <= 0.0: 394570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 3958ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 396d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 397d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 398d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 399d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 400d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 401d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 402d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 403ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 404ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 405ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 4068ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 407311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 408d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 40973ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger if not 1e-7 < u1 < .9999999: 41073ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger continue 41173ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 412d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 413d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 416d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 417b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4210c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 423d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 424b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return -_log(u) * beta 425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 427d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 428d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 429d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 430311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 431d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 433d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 434d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 436d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 437d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 439d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 441d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 442d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 443b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 444b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 445cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 44695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 447d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 448c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gaussian distribution. 449c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 450c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. This is 451c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger slightly faster than the normalvariate() function. 452c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 453c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Not thread-safe without a lock around calls. 454ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 455c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 456d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 457d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 463d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 464d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 465d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 466d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 467d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 468d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 469d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 471d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 472d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 473d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 474d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 475d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 477d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 478d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 479d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 483d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 48595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 486cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 48785e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 48885e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 48985e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 49085e2e4742d0a1accecd02058a7907df36308297eTim Peters## 49185e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 49285e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 49385e2e4742d0a1accecd02058a7907df36308297eTim Peters## 49485e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 49585e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 49685e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 49785e2e4742d0a1accecd02058a7907df36308297eTim Peters## 49885e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 49995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 501c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Beta distribution. 502c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 503c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > -1 and beta} > -1. 504c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Returned values range between 0 and 1. 505ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 506c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 507ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 50885e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 50985e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 51085e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 51185e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 51285e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 51385e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 51485e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 51595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 516cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 517cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 519c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Pareto distribution. alpha is the shape parameter.""" 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 521cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 52273ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 524cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 525cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 526cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 528c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Weibull distribution. 529c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 530c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger alpha is the scale parameter and beta is the shape parameter. 531ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 532c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 534cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 53573ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 536d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5376c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 53840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger## -------------------- Wichmann-Hill ------------------- 53940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 54040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettingerclass WichmannHill(Random): 54140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 54240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 1 # used by getstate/setstate 54340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 54440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def seed(self, a=None): 54540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Initialize internal state from hashable object. 54640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 54740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. 54840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 54940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is not None or an int or long, hash(a) is used instead. 55040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 55140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is an int or long, a is used directly. Distinct values between 55240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 0 and 27814431486575L inclusive are guaranteed to yield distinct 55340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger internal states (this guarantee is specific to the default 55440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Wichmann-Hill generator). 55540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 55640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 55740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 55840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 55940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 56040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = long(time.time() * 256) 56140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 56240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not isinstance(a, (int, long)): 56340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 56440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 56540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 30268) 56640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 30306) 56740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 30322) 56840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = int(x)+1, int(y)+1, int(z)+1 56940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 57040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 57140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 57240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def random(self): 57340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 57440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 57540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichman-Hill random number generator. 57640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 57740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichmann, B. A. & Hill, I. D. (1982) 57840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Algorithm AS 183: 57940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # An efficient and portable pseudo-random number generator 58040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 31 (1982) 188-190 58140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 58240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # see also: 58340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Correction to Algorithm AS 183 58440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 33 (1984) 123 58540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 58640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # McLeod, A. I. (1985) 58740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # A remark on Algorithm AS 183 58840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 34 (1985),198-200 58940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # This part is thread-unsafe: 59140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # BEGIN CRITICAL SECTION 59240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 59340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (171 * x) % 30269 59440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (172 * y) % 30307 59540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (170 * z) % 30323 59640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 59740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # END CRITICAL SECTION 59840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Note: on a platform using IEEE-754 double arithmetic, this can 60040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # never return 0.0 (asserted by Tim; proof too long for a comment). 60140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 60240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def getstate(self): 60440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Return internal state; can be passed to setstate() later.""" 60540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return self.VERSION, self._seed, self.gauss_next 60640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def setstate(self, state): 60840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Restore internal state from object returned by getstate().""" 60940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version = state[0] 61040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 1: 61140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, self._seed, self.gauss_next = state 61240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger else: 61340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("state with version %s passed to " 61440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "Random.setstate() of version %s" % 61540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger (version, self.VERSION)) 61640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def jumpahead(self, n): 61840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Act as if n calls to random() were made, but quickly. 61940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger n is an int, greater than or equal to 0. 62140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Example use: If you have 2 threads and know that each will 62340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger consume no more than a million random numbers, create two Random 62440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger objects r1 and r2, then do 62540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.setstate(r1.getstate()) 62640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.jumpahead(1000000) 62740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Then r1 and r2 will use guaranteed-disjoint segments of the full 62840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger period. 62940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 63040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not n >= 0: 63240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("n must be >= 0") 63340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 63440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = int(x * pow(171, n, 30269)) % 30269 63540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = int(y * pow(172, n, 30307)) % 30307 63640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = int(z * pow(170, n, 30323)) % 30323 63740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 63840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def __whseed(self, x=0, y=0, z=0): 64040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Set the Wichmann-Hill seed from (x, y, z). 64140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger These must be integers in the range [0, 256). 64340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 64440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not type(x) == type(y) == type(z) == int: 64640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise TypeError('seeds must be integers') 64740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 64840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError('seeds must be in range(0, 256)') 64940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if 0 == x == y == z: 65040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 65140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 65240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = long(time.time() * 256) 65340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = int((t&0xffffff) ^ (t>>24)) 65440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, x = divmod(t, 256) 65540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, y = divmod(t, 256) 65640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, z = divmod(t, 256) 65740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Zero is a poor seed, so substitute 1 65840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = (x or 1, y or 1, z or 1) 65940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 66140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def whseed(self, a=None): 66340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Seed from hashable object's hash code. 66440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. It is not guaranteed 66640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger that objects with distinct hash codes lead to distinct internal 66740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger states. 66840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger This is obsolete, provided for compatibility with the seed routine 67040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger used prior to Python 2.1. Use the .seed() method instead. 67140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 67240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 67440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed() 67540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return 67640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 67740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 256) 67840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 256) 67940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 256) 68040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (x + a) % 256 or 1 68140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (y + a) % 256 or 1 68240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (z + a) % 256 or 1 68340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed(x, y, z) 68440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 685cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 686ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 687d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test_generator(n, funccall): 6880c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 6890c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print n, 'times', funccall 6900c9886d589ddebf32de0ca3f027a173222ed383aTim Peters code = compile(funccall, funccall, 'eval') 691b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total = 0.0 6920c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 6930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 6940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 6950c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 6960c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 6970c9886d589ddebf32de0ca3f027a173222ed383aTim Peters x = eval(code) 698b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total += x 6990c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 7000c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 7010c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 7020c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 7030c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 704b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger avg = total/n 705d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 7060c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 7070c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 708ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 709f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 710f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingerdef _test(N=2000): 711d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'random()') 712d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'normalvariate(0.0, 1.0)') 713d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'lognormvariate(0.0, 1.0)') 714d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 715b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger _test_generator(N, 'gammavariate(0.01, 1.0)') 716b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger _test_generator(N, 'gammavariate(0.1, 1.0)') 7178ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters _test_generator(N, 'gammavariate(0.1, 2.0)') 718d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.5, 1.0)') 719d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.9, 1.0)') 720d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(1.0, 1.0)') 721d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(2.0, 1.0)') 722d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(20.0, 1.0)') 723d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(200.0, 1.0)') 724d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gauss(0.0, 1.0)') 725d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'betavariate(3.0, 3.0)') 726cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 727715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 72840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# as module-level functions. The functions share state across all uses 72940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger#(both in the user's code and in the Python libraries), but that's fine 73040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# for most programs and is easier for the casual user than making them 73140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# instantiate their own Random() instance. 73240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 733d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 734d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 735d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 736d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 737d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 738d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 739d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 740f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingersample = _inst.sample 741d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 742d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 743d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 744d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 745d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 746d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 747d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 748d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 749d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 750d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 751d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 752d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 753d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 754d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 755ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 756d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 757