random.py revision 94547f7646895e032f8fc145529d9efc3a70760d
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. 320e1159583c06fdf85d7d2dbe8b82e42565b9d166Tim Peters* It is one of the most extensively tested generators in existence. 330e1159583c06fdf85d7d2dbe8b82e42565b9d166Tim Peters* Without a direct way to compute N steps forward, the semantics of 340e1159583c06fdf85d7d2dbe8b82e42565b9d166Tim Peters jumpahead(n) are weakened to simply jump to another distant state and rely 350e1159583c06fdf85d7d2dbe8b82e42565b9d166Tim Peters on the large period to avoid overlapping sequences. 360e1159583c06fdf85d7d2dbe8b82e42565b9d166Tim Peters* The random() method is implemented in C, executes in a single Python step, 370e1159583c06fdf85d7d2dbe8b82e42565b9d166Tim Peters and is, therefore, threadsafe. 3840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 39e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum""" 40d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 412f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingerfrom warnings import warn as _warn 422f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingerfrom types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 4391e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettingerfrom math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 44d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 45c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettingerfrom os import urandom as _urandom 46c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettingerfrom binascii import hexlify as _hexlify 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", 52356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 5323f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger "SystemRandom"] 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 607c2a85b2d44851c2442ade579b760f86447bf848Tim PetersRECIP_BPF = 2**-BPF 6133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 62356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 63d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 6440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# Adrian Baddeley. Adapted by Raymond Hettinger for use with 653fa19d7ff89be87139e2864fb9186b424d180a58Raymond Hettinger# the Mersenne Twister and os.urandom() core generators. 6633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 67145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerimport _random 6840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 69145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerclass Random(_random.Random): 70c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Random number generator base class used by bound module functions. 71c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 72c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Used to instantiate instances of Random to get generators that don't 73c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger share state. Especially useful for multi-threaded programs, creating 74c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger a different instance of Random for each thread, and using the jumpahead() 75c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger method to ensure that the generated sequences seen by each thread don't 76c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger overlap. 77c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 78c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Class Random can also be subclassed if you want to use a different basic 79c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger generator of your own devising: in that case, override the following 80c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger methods: random(), seed(), getstate(), setstate() and jumpahead(). 812f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Optionally, implement a getrandombits() method so that randrange() 822f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger can cover arbitrarily large ranges. 83ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 84c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 8533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 8640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 2 # used by getstate/setstate 8733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 88d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 89d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 9033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 91d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 92d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 9333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 94d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 9540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 96ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 970de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 980de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 99d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 10023f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger None or no argument seeds from current time or from an operating 10123f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger system specific randomness source if available. 1020de88fc4b108751b86443852b6741680d704168fTim Peters 103bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 104d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 105d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1063081d59f920229b26293c7a3ee3f1a9da0da53e9Raymond Hettinger if a is None: 107c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger try: 108c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger a = long(_hexlify(_urandom(16)), 16) 109c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger except NotImplementedError: 110356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger import time 111356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger a = long(time.time() * 256) # use fractional seconds 112356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 113145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).seed(a) 11446c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 11546c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 116d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 117d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 118145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger return self.VERSION, super(Random, self).getstate(), self.gauss_next 119d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 120d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 121d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 12340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 2: 12440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, internalstate, self.gauss_next = state 125145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).setstate(internalstate) 126d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 127d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 128d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 129d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 130d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 131cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 132cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 133d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 134cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 135d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 136cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 137cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 138d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 139cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 140cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 141cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 1425f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger def __reduce__(self): 1435f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger return self.__class__, (), self.getstate() 1445f078ff7f0c6bb5086fae077379fc79729c34d2dRaymond Hettinger 145cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 146d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 1472f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger def randrange(self, start, stop=None, step=1, int=int, default=None, 1482f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger maxwidth=1L<<BPF): 149d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 150d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 151d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 152d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 1532f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Do not supply the 'int', 'default', and 'maxwidth' arguments. 154d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 155d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 156d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 1579146f27b7799dab231083f194a14c6157b57549fTim Peters # common case while still doing adequate error checking. 158d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 159d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 160d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 161d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 162d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 1632f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if istart >= maxwidth: 1642f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return self._randbelow(istart) 165d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 166d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1679146f27b7799dab231083f194a14c6157b57549fTim Peters 1689146f27b7799dab231083f194a14c6157b57549fTim Peters # stop argument supplied. 169d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 170d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 171d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 1722f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger width = istop - istart 1732f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if step == 1 and width > 0: 17476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # Note that 1752f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # int(istart + self.random()*width) 17676ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # instead would be incorrect. For example, consider istart 17776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # = -2 and istop = 0. Then the guts would be in 17876ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # -2.0 to 0.0 exclusive on both ends (ignoring that random() 17976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # might return 0.0), and because int() truncates toward 0, the 18076ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # final result would be -1 or 0 (instead of -2 or -1). 1812f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # istart + int(self.random()*width) 18276ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # would also be incorrect, for a subtler reason: the RHS 18376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # can return a long, and then randrange() would also return 18476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # a long, but we're supposed to return an int (for backward 18576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # compatibility). 1862f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 1872f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if width >= maxwidth: 18858eb11cf62dd04ccc2c364b62fd51b4265e2e203Tim Peters return int(istart + self._randbelow(width)) 1892f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(istart + int(self.random()*width)) 190d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 1912f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 1929146f27b7799dab231083f194a14c6157b57549fTim Peters 1939146f27b7799dab231083f194a14c6157b57549fTim Peters # Non-unit step argument supplied. 194d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 196d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 197d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 198ffdb8bb99c4017152a9dca70669f9d6b9831d454Raymond Hettinger n = (width + istep - 1) // istep 199d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 200ffdb8bb99c4017152a9dca70669f9d6b9831d454Raymond Hettinger n = (width + istep + 1) // istep 201d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 205d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 2062f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2072f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if n >= maxwidth: 20894547f7646895e032f8fc145529d9efc3a70760dRaymond Hettinger return istart + istep*self._randbelow(n) 209d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 211d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 212cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 213d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 214d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 215d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 216d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 2172f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, 2182f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 2192f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger """Return a random int in the range [0,n) 2202f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2212f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger Handles the case where n has more bits than returned 2222f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger by a single call to the underlying generator. 2232f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger """ 2242f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 2252f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger try: 2262f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger getrandbits = self.getrandbits 2272f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger except AttributeError: 2282f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger pass 2292f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger else: 2302f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # Only call self.getrandbits if the original random() builtin method 2312f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # has not been overridden or if a new getrandbits() was supplied. 2322f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger # This assures that the two methods correspond. 2332f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 2342f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 2352f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger r = getrandbits(k) 2362f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger while r >= n: 2372f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger r = getrandbits(k) 2382f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return r 2392f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger if n >= _maxwidth: 2402f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger _warn("Underlying random() generator does not supply \n" 2412f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger "enough bits to choose from a population range this large") 2422f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger return int(self.random() * n) 2432f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettinger 244cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 245cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 246d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 247d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 2485dae505bbd59641a948c81bea981e7c44d4c2343Raymond Hettinger return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 249d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 250d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 251d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 252d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 253d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 254d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 255d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 256d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 257d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 258d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 25985c20a41dfcec04d161ad7da7260e7b94c62d228Raymond Hettinger for i in reversed(xrange(1, len(x))): 260cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 261d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 262d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 263d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 264fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger def sample(self, population, k): 265f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """Chooses k unique random elements from a population sequence. 266f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 267c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Returns a new list containing elements from the population while 268c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger leaving the original population unchanged. The resulting list is 269c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger in selection order so that all sub-slices will also be valid random 270c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger samples. This allows raffle winners (the sample) to be partitioned 271c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger into grand prize and second place winners (the subslices). 272f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 273c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Members of the population need not be hashable or unique. If the 274c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger population contains repeats, then each occurrence is a possible 275c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger selection in the sample. 276f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 277c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger To choose a sample in a range of integers, use xrange as an argument. 278c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger This is especially fast and space efficient for sampling from a 279c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger large population: sample(xrange(10000000), 60) 280f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """ 281f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 282c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX Although the documentation says `population` is "a sequence", 283c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX attempts are made to cater to any iterable with a __len__ 284c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX method. This has had mixed success. Examples from both 285c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX sides: sets work fine, and should become officially supported; 286c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX dicts are much harder, and have failed in various subtle 287c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX ways across attempts. Support for mapping types should probably 288c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX be dropped (and users should pass mapping.keys() or .values() 289c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # XXX explicitly). 290c17976e9833f3093adb1019356737e728a24f7c9Tim Peters 291c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger # Sampling without replacement entails tracking either potential 29291e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger # selections (the pool) in a list or previous selections in a set. 293c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 2942b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # When the number of selections is small compared to the 2952b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # population, then tracking selections is efficient, requiring 29691e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger # only a small set and an occasional reselection. For 2972b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # a larger number of selections, the pool tracking method is 2982b55d35850e3e8e0b28aba7878d3f9122a7907acJeremy Hylton # preferred since the list takes less space than the 29991e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger # set and it doesn't suffer from frequent reselections. 300c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 301f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger n = len(population) 302f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if not 0 <= k <= n: 303f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger raise ValueError, "sample larger than population" 3048b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger random = self.random 305fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger _int = int 306c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result = [None] * k 30791e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger setsize = 21 # size of a small set minus size of an empty list 30891e27c253c8bb8b6ae8521f1dbb76de7c66ad8cfRaymond Hettinger if k > 5: 3099e34c047325651853a95f95e538582a4f6d5b7f6Tim Peters setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 310c17976e9833f3093adb1019356737e728a24f7c9Tim Peters if n <= setsize or hasattr(population, "keys"): 311c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # An n-length list is smaller than a k-length set, or this is a 312c17976e9833f3093adb1019356737e728a24f7c9Tim Peters # mapping type so the other algorithm wouldn't work. 313311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger pool = list(population) 314311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger for i in xrange(k): # invariant: non-selected at [0,n-i) 315fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * (n-i)) 316311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger result[i] = pool[j] 3178b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger pool[j] = pool[n-i-1] # move non-selected item into vacancy 318c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger else: 31966d09f1b3029d9cf975ccf26c437c9fb2605db91Raymond Hettinger try: 3203c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger selected = set() 3213c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger selected_add = selected.add 3223c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger for i in xrange(k): 323fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 3243c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger while j in selected: 3253c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger j = _int(random() * n) 3263c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger selected_add(j) 3273c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger result[i] = population[j] 328c17976e9833f3093adb1019356737e728a24f7c9Tim Peters except (TypeError, KeyError): # handle (at least) sets 3293c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger if isinstance(population, list): 3303c3346daa9bf900080428ed12b6d83aa04f7332bRaymond Hettinger raise 331c17976e9833f3093adb1019356737e728a24f7c9Tim Peters return self.sample(tuple(population), k) 332311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger return result 333f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 334cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 335cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 336cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 337d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 338d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 339d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 340d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 341ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 342cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 343ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 345c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Normal distribution. 346c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 347c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. 348ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 349c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 350d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 351d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 352d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 353d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 354d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 355d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 356d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 357d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 35842406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger while 1: 359d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 36073ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 361d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 362d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 364d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 365d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 366ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 367cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 368ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 369d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 370c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Log normal distribution. 371c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 372c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger If you take the natural logarithm of this distribution, you'll get a 373c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger normal distribution with mean mu and standard deviation sigma. 374c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu can have any value, and sigma must be greater than zero. 375ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 376c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 377d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 378ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 379cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 380ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 381d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 382c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Exponential distribution. 383c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 384c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger lambd is 1.0 divided by the desired mean. (The parameter would be 385c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger called "lambda", but that is a reserved word in Python.) Returned 386c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger values range from 0 to positive infinity. 387ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 388c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 390d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 391ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 392d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3930c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 395d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 396d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 397ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 398cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 399ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 400d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 401c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular data distribution. 402ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 403c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean angle, expressed in radians between 0 and 2*pi, and 404c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger kappa is the concentration parameter, which must be greater than or 405c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger equal to zero. If kappa is equal to zero, this distribution reduces 406c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger to a uniform random angle over the range 0 to 2*pi. 407ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 408c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 409d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 410d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 411d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 4125810297052003f28788f6790ac799fe8e5373494Guido van Rossum 413d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 414d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 415d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 4165810297052003f28788f6790ac799fe8e5373494Guido van Rossum 417d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 4195810297052003f28788f6790ac799fe8e5373494Guido van Rossum 420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 421d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 423ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 424d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 425d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 426d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 427ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 42842406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger while 1: 429d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 430ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 431d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 433d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 434ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 436ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 43742406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c): 438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 439ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 441d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 442d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 443d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 444d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 445ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 446d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 447ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 448cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 449ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 450d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 451c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gamma distribution. Not the gamma function! 452c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 453c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > 0 and beta > 0. 454c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 455c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 4568ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 457b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 4588ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 459570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # Warning: a few older sources define the gamma distribution in terms 460570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # of alpha > -1.0 461570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum if alpha <= 0.0 or beta <= 0.0: 462570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 4638ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 464d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 465d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 466d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 467d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 468d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 469d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 470d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 471ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 472ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 473ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 4748ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 47542406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger while 1: 476d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 47773ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger if not 1e-7 < u1 < .9999999: 47873ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger continue 47973ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 480d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 481d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 482d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 483d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 484d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 485b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 486d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 487d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4890c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 490d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 491d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 492b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return -_log(u) * beta 493d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 494d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 495d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 496d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 49842406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger while 1: 499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 501d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 502d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 50342406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger x = p ** (1.0/alpha) 504d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 505d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 506d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 50742406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger if p > 1.0: 50842406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger if u1 <= x ** (alpha - 1.0): 50942406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger break 51042406e6f27e9a42e91db8706d897e0b478b13a4dRaymond Hettinger elif u1 <= _exp(-x): 511d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 512b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 513b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 514cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 51595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 516d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 517c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gaussian distribution. 518c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 519c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. This is 520c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger slightly faster than the normalvariate() function. 521c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 522c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Not thread-safe without a lock around calls. 523ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 524c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 526d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 527d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 528d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 529d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 530d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 531d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 532d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 533d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 534d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 535d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 536d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 537d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 538d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 539d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 540d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 541d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 542d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 543d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 544d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 545d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 546d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 547d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 548d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 549d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 550d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 551d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 552d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 553d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 55495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 555cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 55685e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 55785e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 55885e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 55985e2e4742d0a1accecd02058a7907df36308297eTim Peters## 56085e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 56185e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 56285e2e4742d0a1accecd02058a7907df36308297eTim Peters## 56385e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 56485e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 56585e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 56685e2e4742d0a1accecd02058a7907df36308297eTim Peters## 56785e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 56895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 569d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 570c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Beta distribution. 571c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 572c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > -1 and beta} > -1. 573c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Returned values range between 0 and 1. 574ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 575c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 576ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 57785e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 57885e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 57985e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 58085e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 58185e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 58285e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 58385e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 58495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 585cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 586cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 587d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 588c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Pareto distribution. alpha is the shape parameter.""" 589d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 590cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 59173ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 592d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 593cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 594cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 595cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 596d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 597c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Weibull distribution. 598c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 599c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger alpha is the scale parameter and beta is the shape parameter. 600ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 601c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 602d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 603cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 60473ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 605d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 6066c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 60740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger## -------------------- Wichmann-Hill ------------------- 60840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettingerclass WichmannHill(Random): 61040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 1 # used by getstate/setstate 61240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def seed(self, a=None): 61440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Initialize internal state from hashable object. 61540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61623f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger None or no argument seeds from current time or from an operating 61723f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger system specific randomness source if available. 61840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is not None or an int or long, hash(a) is used instead. 62040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is an int or long, a is used directly. Distinct values between 62240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 0 and 27814431486575L inclusive are guaranteed to yield distinct 62340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger internal states (this guarantee is specific to the default 62440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Wichmann-Hill generator). 62540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 62640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 62740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 628c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger try: 629c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger a = long(_hexlify(_urandom(16)), 16) 630c1c43cad63a88eae694b174c9a0fe6242dd5972bRaymond Hettinger except NotImplementedError: 631356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger import time 632356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger a = long(time.time() * 256) # use fractional seconds 63340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not isinstance(a, (int, long)): 63540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 63640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 30268) 63840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 30306) 63940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 30322) 64040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = int(x)+1, int(y)+1, int(z)+1 64140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 64340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def random(self): 64540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 64640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichman-Hill random number generator. 64840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 64940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichmann, B. A. & Hill, I. D. (1982) 65040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Algorithm AS 183: 65140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # An efficient and portable pseudo-random number generator 65240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 31 (1982) 188-190 65340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 65440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # see also: 65540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Correction to Algorithm AS 183 65640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 33 (1984) 123 65740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 65840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # McLeod, A. I. (1985) 65940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # A remark on Algorithm AS 183 66040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 34 (1985),198-200 66140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # This part is thread-unsafe: 66340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # BEGIN CRITICAL SECTION 66440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 66540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (171 * x) % 30269 66640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (172 * y) % 30307 66740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (170 * z) % 30323 66840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 66940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # END CRITICAL SECTION 67040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Note: on a platform using IEEE-754 double arithmetic, this can 67240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # never return 0.0 (asserted by Tim; proof too long for a comment). 67340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 67440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def getstate(self): 67640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Return internal state; can be passed to setstate() later.""" 67740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return self.VERSION, self._seed, self.gauss_next 67840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def setstate(self, state): 68040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Restore internal state from object returned by getstate().""" 68140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version = state[0] 68240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 1: 68340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, self._seed, self.gauss_next = state 68440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger else: 68540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("state with version %s passed to " 68640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "Random.setstate() of version %s" % 68740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger (version, self.VERSION)) 68840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def jumpahead(self, n): 69040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Act as if n calls to random() were made, but quickly. 69140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 69240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger n is an int, greater than or equal to 0. 69340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 69440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Example use: If you have 2 threads and know that each will 69540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger consume no more than a million random numbers, create two Random 69640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger objects r1 and r2, then do 69740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.setstate(r1.getstate()) 69840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.jumpahead(1000000) 69940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Then r1 and r2 will use guaranteed-disjoint segments of the full 70040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger period. 70140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 70240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 70340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not n >= 0: 70440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("n must be >= 0") 70540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 70640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = int(x * pow(171, n, 30269)) % 30269 70740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = int(y * pow(172, n, 30307)) % 30307 70840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = int(z * pow(170, n, 30323)) % 30323 70940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 71040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def __whseed(self, x=0, y=0, z=0): 71240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Set the Wichmann-Hill seed from (x, y, z). 71340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger These must be integers in the range [0, 256). 71540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 71640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not type(x) == type(y) == type(z) == int: 71840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise TypeError('seeds must be integers') 71940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 72040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError('seeds must be in range(0, 256)') 72140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if 0 == x == y == z: 72240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 72340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 72440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = long(time.time() * 256) 72540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = int((t&0xffffff) ^ (t>>24)) 72640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, x = divmod(t, 256) 72740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, y = divmod(t, 256) 72840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, z = divmod(t, 256) 72940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Zero is a poor seed, so substitute 1 73040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = (x or 1, y or 1, z or 1) 73140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 73240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 73340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 73440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def whseed(self, a=None): 73540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Seed from hashable object's hash code. 73640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 73740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. It is not guaranteed 73840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger that objects with distinct hash codes lead to distinct internal 73940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger states. 74040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 74140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger This is obsolete, provided for compatibility with the seed routine 74240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger used prior to Python 2.1. Use the .seed() method instead. 74340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 74440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 74540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 74640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed() 74740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return 74840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 74940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 256) 75040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 256) 75140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 256) 75240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (x + a) % 256 or 1 75340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (y + a) % 256 or 1 75440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (z + a) % 256 or 1 75540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed(x, y, z) 75640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 75723f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger## --------------- Operating System Random Source ------------------ 758356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 75923f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettingerclass SystemRandom(Random): 76023f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger """Alternate random number generator using sources provided 76123f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger by the operating system (such as /dev/urandom on Unix or 76223f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger CryptGenRandom on Windows). 763356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 764356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger Not available on all systems (see os.urandom() for details). 765356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """ 766356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 767356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def random(self): 768356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 7697c2a85b2d44851c2442ade579b760f86447bf848Tim Peters return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF 770356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 771356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def getrandbits(self, k): 772356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger """getrandbits(k) -> x. Generates a long int with k random bits.""" 773356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if k <= 0: 774356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise ValueError('number of bits must be greater than zero') 775356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger if k != int(k): 776356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger raise TypeError('number of bits should be an integer') 777356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger bytes = (k + 7) // 8 # bits / 8 and rounded up 778356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger x = long(_hexlify(_urandom(bytes)), 16) 779356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger return x >> (bytes * 8 - k) # trim excess bits 780356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 781356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def _stub(self, *args, **kwds): 78223f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger "Stub method. Not used for a system random number generator." 783356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger return None 784356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger seed = jumpahead = _stub 785356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 786356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger def _notimplemented(self, *args, **kwds): 78723f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger "Method should not be called for a system random number generator." 78823f1241dc6495eb255e1a389aef204a3e35a2632Raymond Hettinger raise NotImplementedError('System entropy source does not have state.') 789356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger getstate = setstate = _notimplemented 790356a4599acd4c835ab88c221bd5da073c9895e83Raymond Hettinger 791cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 792ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 79362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettingerdef _test_generator(n, func, args): 7940c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 79562297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger print n, 'times', func.__name__ 796b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total = 0.0 7970c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 7980c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 7990c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 8000c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 8010c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 80262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger x = func(*args) 803b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total += x 8040c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 8050c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 8060c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 8070c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 8080c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 809b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger avg = total/n 810d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 8110c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 8120c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 813ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 814f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 815f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingerdef _test(N=2000): 81662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, random, ()) 81762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, normalvariate, (0.0, 1.0)) 81862297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, lognormvariate, (0.0, 1.0)) 81962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, vonmisesvariate, (0.0, 1.0)) 82062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.01, 1.0)) 82162297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.1, 1.0)) 82262297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.1, 2.0)) 82362297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.5, 1.0)) 82462297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (0.9, 1.0)) 82562297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (1.0, 1.0)) 82662297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (2.0, 1.0)) 82762297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (20.0, 1.0)) 82862297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gammavariate, (200.0, 1.0)) 82962297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, gauss, (0.0, 1.0)) 83062297132215490e9cb406e1a21f03aff40d421cbRaymond Hettinger _test_generator(N, betavariate, (3.0, 3.0)) 831cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 832715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 83340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# as module-level functions. The functions share state across all uses 83440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger#(both in the user's code and in the Python libraries), but that's fine 83540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# for most programs and is easier for the casual user than making them 83640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# instantiate their own Random() instance. 83740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 838d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 839d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 840d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 841d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 842d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 843d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 844d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 845f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingersample = _inst.sample 846d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 847d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 848d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 849d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 850d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 851d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 852d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 853d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 854d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 855d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 856d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 857d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 858d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 8592f726e9093381572b21edbfc42659ea89fbdf686Raymond Hettingergetrandbits = _inst.getrandbits 860d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 861ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 862d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 863